Analysis

Power Spectrum Measurement

The code can measure the power spectrum P(k) of a periodic simulation. It uses grid interlacing (see Hockney and Eastwood [1], section 7-8 “Interlacing”) and high order mass assignment (see Sefusatti et al. [9]).

Parameters

The following parameters control P(k) measurement:

nGridPk (default 0)integer

This size of the 3D grid for measuring the power spectrum. This feature is disabled if the grid size is zero which is the default. A reasonable value is to set this to the same size as the IC grid, plus or minus a factor of two.

iPkInterval (default 1)integer

number of timesteps between pk outputs

iPkOrder (default 3)IPKORDER

Selects the mass assignment order. The models should be imported with:

from PKDGRAV import ASSIGNMENT_ORDER

Possible values are: - ASSIGNMENT_ORDER.NGP: Nearest Grid Point - ASSIGNMENT_ORDER.CIC: Cloud in Cell - ASSIGNMENT_ORDER.TSC: Triangular Shaped Cloud - ASSIGNMENT_ORDER.PCS: Piecewise Cubic Spline

The recommended method is PCS. See Sefusatti et al. [9] for details.

bPkInterlace (default True)Boolean

If enabled, the power spectrum measurement tool creates two mass assignment grids. The second contains the mass with the particles offset by half a grid cell. The results of the two mass assignment grids are combined to reduce grid aliasing effects. See Hockney and Eastwood [1], section 7-8 “Interlacing” for details..

Output

The code outputs a pk file for each selected step with four columns:

k

The \(k\) bin in \(h\text{Mpc}^{-1}\). If dBoxSize has been specified (often the case when using the build-in IC generator), \(k\) will be in physical units (scale by dBoxSize). Otherwise you need to divide by your box length.

P(k)

The measured P(k) in \(h^{-1}\text{Mpc}^{3}\). Again, if dBoxSize was not specified then this will need to be corrected by multiplying by \(\texttt{dBoxSize}^3\).

Number of Samples

This is the number of Delta(k) cells that went into the P(k) measurement for each bin.

Linear P(k)

This is the same as the P(k) column, except it includes the linear part (if present), for example if you are using the linear neutrino treatment.

Light Cone

The code can output light cone data as either healpix maps, or raw particle data (or both).

Parameters

bLightCone (default False)Boolean

Enable or disable light cone output.

dRedshiftLCP (default 0.0)float

starting redshift to output light cone

nSideHealpix (default 0)integer

Sets the “nSide” parameter for healpix, and if non-zero enables the output of the healpix maps.

bLightConeParticles (default False)Boolean

output light cone particles

bBowtie (default False)Boolean

output +++ and — octants of the cone; a bowtie

sqdegLCP (default 50.0)float

square degrees of lightcone (if (<0 || >41253) then all sky)

hLCP (default [0.749, 0.454, 1.0])

lightcone direction vector

Output

Healpix

The healpix output is a binary format consisting of 32-bit integer counts of grouped and ungrouped particles (grouped is only valid if group finding is enabled) as well as the potential. The following script (found in tools/hpb2fits.py) will convert the file to fits format.

import numpy as np
import pandas as pd
from astropy.io import fits
from astropy.table import Table
from sys import argv,exit

hpb_type = np.dtype([('grouped', '=i4'),('ungrouped', '=i4'), ('potential', '=f4')])

bPotential = True

NSIDE = 64
N = 12 * NSIDE**2

hdr = fits.Header()
hdr['EXTEND'] = 'T'
primary_hdu = fits.PrimaryHDU(header=hdr)
primary_hdu.writeto('test.fits')


# Use BinTableHDU as a template
hdr = fits.BinTableHDU(Table(names=['SIGNAL'],dtype=['=f4' if bPotential else '=i4']),name='BINTABLE').header
hdr['ORDERING'] = ("RING","Pixel ordering scheme, either RING or NESTED")
hdr['INDXSCHM'] = ("IMPLICIT","Pixel indexing scheme (IMPLICIT or EXPLICIT)")
hdr['NSIDE']    = (NSIDE,"Resolution parameter for HEALPIX")
hdr['COORDSYS'] = ("C","Pixelisation coordinate system")
hdr['PIXTYPE']  = ("HEALPIX", "HEALPIX Pixelisation")
hdr['NAXIS']    = 2
hdr['NAXIS2']   = N
hdr['NAXIS1']   = 1
hdr['BITPIX']   = -32 if bPotential else 32

hdu = fits.StreamingHDU('test.fits',hdr)
for file in argv[1:]:
    with open(file,'rb') as hpb:
        data = np.fromfile(hpb,dtype=hpb_type)
        data = pd.DataFrame(data,columns=data.dtype.names)
        if bPotential:
            hdu.write(data['potential'].values)
        else:
            data = (data.grouped+data.ungrouped).to_frame()
            data.columns=['signal']
            hdu.write(data['signal'].values)
hdu.close()

Particles

The particle output is also a binary format. Note that files may be empty if the assigned processor does not output any particles during the step. This is normal. Each output particle is 40 bytes as follows.

Offset	Type	Field
0	64-bit integer	particle ID
8	3 x float	position (x y z)
20	3 x float	velocity (x y z)
32	float	potential
36	4 bytes	padding

Friend of Friends Group Finder

The code can output FoF groups at each step. The following parameters control the FoF Group Finder:

bFindGroups (default False)Boolean

The Friends-of-friends group finder is enabled when this is true.

dTau (default 0.0)float

The linking length \(\tau\). This is normally set to one fifth of the mean particle separation. For example, if you use the built-in IC generator you could set this parameter with:

dTau = 0.2 / nGrid

nMinMembers (default 10)integer

After groups have been identified, and that have fewer than this many particles are removed.

dEnvironment0 (default disabled)float

first radius for density environment about a group

dEnvironment1 (default disabled)float

second radius for density environment about a group

The output is a binary format with the following structure.

Offset	Type	Field
0	float[3]	Position (x y z) of deepest potential
12	float	Value of deepest potential
16	float[3]	Shrinking Sphere Center (x y z) (see Power et al. [8])
28	float[3]	Center of Mass Position (x y z)
40	float[3]	Center of Mass Velocity
52	float[3]	Angular momentum
64	float[6]	Moment of Inertia
88	float	Velocity Dispersion \(\sigma\)
92	float	\(r_{\text{max}}\) (maximum distance from center of mass)
96	float	Mass
100	float	Mass enclosed in dEnvironment0
104	float	Mass enclosed in dEnvironment1
108	float	Half mass radius
112	int	Number of black holes
116	int	Number of stars
120	int	Number of gas particles
124	int	Number of dark matter particles
128	uint64_t	Global group ID
136

General Analysis Hooks

Aside from the legacy analysis functions, it is possible to add custom analysis hooks. You can use PKDGRAV.add_analysis() to add any object instance that is callable by Python. The callable object will be passed the following named parameters.

Name	Description
msr	the PKDGRAV module
step	current integer step number
time	simulation time
a	the expansion factor (for cosmological simulations)
theta	currently calculated theta

For example:

from PKDGRAV import add_analysis, ASSIGNMENT_ORDER

class MassGrid:
    grid = 0
    order = ASSIGNMENT_ORDER.PCS
    def __init__(self,grid,order=ASSIGNMENT_ORDER.PCS):
        self.grid = grid
        self.order = order
    def __call__(self,msr,step,time,**kwargs):
        print('calculating density grid')
        msr.grid_create(self.grid)
        msr.assign_mass(order=self.order)
        msr.grid_write('output.{:05d}'.format(step))
        msr.grid_delete()
    def ephemeral(self,msr,**kwargs):
        return msr.grid_ephemeral(self.grid)

This class can be used to output a density grid at each step. To enable it, use:

add_analysis(MassGrid(nGrid))