Tutorial

Cosmological zoom-in simulations

1) Generate low-resolution initial conditions.

Download MUSIC, a code to generate nested grid initial conditions for zoom-in cosmological simulations.

https://bitbucket.org/ohahn/music

As a first step, create dark-matter only low-resolution initial conditions. After compiling MUSIC, edit a parameter file "ics.conf" as follows. Perform the code using the command "./MUSIC ics.conf".

[setup]

boxlength = 50

zstart = 125

levelmin = 8

levelmin_TF = 8

levelmax = 8

padding = 10

overlap = 5

ref_center = 0.5, 0.5, 0.5

ref_extent = 0.2, 0.2, 0.2

align_top = no

baryons = no

use_2LPT = no

use_LLA = no

periodic_TF = yes

[cosmology]

Omega_m = 0.264745

Omega_L = 0.735255

w0 = -1.0

wa = 0.0

Omega_b = 0.0455

H0 = 71

sigma_8 = 0.811

nspec = 0.961

transfer = eisenstein

[random]

seed[7] = 32455

seed[8] = 23456

seed[9] = 34567

seed[10] = 45678

seed[11] = 56789

seed[12] = 67890

seed[13] = 34234

[output]

##Gadget-2 (type=1: high-res particles, type=5: rest)

format = gadget2

filename = ics_gadget_MW_dm3.dat

[poisson]

fft_fine = yes

accuracy = 1e-5

pre_smooth = 3

post_smooth = 3

smoother = gs

laplace_order = 6

grad_order = 6

2) Run a dark-matter only simulation down to z=0 using the created initial conditions.

Download GADGET to run dark-matter only simulations with the initial conditions generated from Step 1.

Be careful when changing parameters, in particular for the following values:

BoxSize 50 % Mpc/h

SofteningHalo 0.0078 % L / N /25, L=50, N=2^8

SofteningHaloMaxPhys 0.002 % (1024/N) * (L/100) x 1e-3

3) Make MW-like halo catalogs using a halo finder, rockstar.

3-1) parameter file

https://bitbucket.org/gfcstanford/rockstar

Download and compile rockstar, then edit a configuration file, "quickstart.cfg". What you should change in quickstart.cfg is "PARTICLE_MASS" and "FORCE_RES".

Example (quickstart.cfg):

FILE_FORMAT = "ASCII"

PARTICLE_MASS = 5.47409e+08 # must specify (in Msun/h) for ART or ASCII

h0 = 0.71

Ol = 0.735255

Om = 0.264745

MIN_HALO_OUTPUT_SIZE = 100

OUTPUT_FORMAT = BOTH

DELETE_BINARY_OUTPUT_AFTER_FINISHED = 1

FORCE_RES = 0.002 #Force resolution of simulation, in Mpc/h

3-2) Make an input file for rockstar

You can read the final output data (Gadget file at z=0) from Step 2 and make an ascii file. The file format is assumed to be:

# X Y Z VX VY VZ ID

Make your own script. Or, here is an example IDL script.

pro example

step = 303

dir = "/rsgrps/myjeon/output_MW/snapshot_"+string(step,format='(i3.3)')+"/"

file = dir+'snap_'+string(step,format='(i3.3)')+'.hdf5'

file_id = H5F_OPEN(file)

pc = 3.085678d18

mpc = pc * 1d6

pos = READ_GADGET_HDF5(file, 1, 'Coordinates')

vel_dm = READ_GADGET_HDF5(file, 1, 'Velocity')

dmmass = read_gadget_hdf5(file, 1, 'Mass', /floatunits)

ID_dm = read_gadget_hdf5(file, 1, 'ParticleIDs')

pos = pos/mpc*0.71

vel_dm = vel_dm/1d5

print, min(dmmass)*0.71

openw, 1, 'rockstar_MW_303.dat'

printf, 1, "#a =", atime

for i=0L, n_elements(dmmass)-1 do begin

printf, 1, pos(0,i), pos(1,i), pos(2,i), vel_dm(0,i), vel_dm(1,i), vel_dm(2,i), ID_dm(i), format='(6(e18.8), i15)'

endfor

close, 1

stop

end

3-3) Run Rockstar

Then, you can simply run the code using the following command.

"./rockstar -c quickstart.cfg rockstar_MW_303.dat"

3-4) Output and create MW-like halo catalogs

Once it is done, the output date will be stored in "halos_0.0.ascii". Read the file and generate MW-like halo catalogs with a virial mass between 7x10^11 M_sun and 10^12 M_sun. Here is an IDL script to generate the catalogs.

pro generate_MW_catalogs

readcol, 'data_ascii', id, num_p, mvir, mbound_vir, rvir, vmax, rvmax, vrms, x, y, z, vx, vy, vz, Count=n, /silent

mvir = mvir / 0.71

mpc = 3.085678d18 * 1d6

massive = where(mvir gt 7d11 and mvir lt 1d12, nm)

x = x*mpc/0.71 & y = y*mpc/0.71 & z = z*mpc/0.71

rvir = rvir*mpc/1d3/0.71 ; in kpc

file = 'center_MW'

openw, 1, file

for inow=0, nm-1 do begin

cx = x(massive(inow))

cy = y(massive(inow))

cz = z(massive(inow))

mmvir = mvir(massive(inow))

rrvir = rvir(massive(inow))

num_p_now = num_p(massive(inow))

printf, 1, cx, cy, cz, mmvir, rrvir, format='(5(e18.8),2x)'

print, cx, cy, cz, mmvir, rrvir, num_p_now, format='(6(e18.8),2x)'

endfor

close, 1

stop

end

4) Search a target halo and find the central position of the target halo.

Prepare an input file (center_MW: Step 3-4) that includes the position of halos using rockstar. Identify all particles within 5 R_vir (R_vir is a virial radius of a halo at z=0) and then trace these particles back in the first snapshot (snapshot_000.hdf5). The information what you need is the center and the extent of the zoom-in regions. You can make your own script. Here is an IDL script.

pro example2

; positions of halos found from rockstar.

file = 'center_MW'

readcol,file,Format='D, D, D, D, D', fgx, fgy, fgz, fgmvir, fgrvir, Count=nfg, /SILENT

step = 303

dir = "/rsgrps/myjeon/output_MW/snapshot_"+string(step,format='(i3.3)')+"/"

file = dir+'snap_'+string(step,format='(i3.3)')+'.hdf5'

file_id = H5F_OPEN(file)

group_id = h5g_open(file_id,"/Header")

redshift = read_attribute_hdf5(group_id, "Redshift")

atime = 1./(1.+redshift)

pos = READ_GADGET_HDF5(file, 1, 'Coordinates')

pc = 3.085678d18

mpc = pc * 1d6

vel_dm = READ_GADGET_HDF5(file, 1, 'Velocity')

dmmass = read_gadget_hdf5(file, 1, 'Mass', /floatunits)

ID_dm = read_gadget_hdf5(file, 1, 'ParticleIDs')

; read the initial position at z=125 from snapshot_000.hdf5

for i=0, nfg-1 do begin

cx = fgx(i)

cy = fgy(i)

cz = fgz(i)

mvir = fgmvir(i)

rvir = fgrvir(i)*5.0

dx = cx - pos(0,*)

dy = cy - pos(1,*)

dz = cz - pos(2,*)

dr = sqrt(dx*dx + dy*dy + dz*dz)

valid = where(dr le rvir, nvalid)

print, i, 'nvalid', nvalid

ID_save = ID_dm(valid)

step = 0

dir = "/rsgrps/myjeon/output_MW/snapshot_"+string(step,format='(i3.3)')+"/"

file = dir+'snap_'+string(step,format='(i3.3)')+'.hdf5'

pos1 = READ_GADGET_HDF5(file, 1, 'Coordinates', /noconv)

ID_dm1 = read_gadget_hdf5(file, 1, 'ParticleIDs')

ID_now = indgen(nvalid)

xnow = fltarr(nvalid) & ynow = fltarr(nvalid) & znow = fltarr(nvalid)

match, ID_save, ID_dm1, suba, subb

xnow = reform(pos1(0,subb))

ynow = reform(pos1(1,subb))

znow = reform(pos1(2,subb))

ID_now = reform(ID_dm1(subb))

xext = max(xnow) - min(xnow)

yext = max(ynow) - min(ynow)

zext = max(znow) - min(znow)

xcen = xext/2.0+min(xnow)

ycen = yext/2.0+min(ynow)

zcen = zext/2.0+min(znow)

box = 50.0

if xext/box le 0.1 and yext/box le 0.1 and zext/box le 0.1 then begin

print, xcen/box, ycen/box, zcen/box

print, xext/box, yext/box, zext/box

endif

endfor

stop

end

5) Re-generate zoom-in initial conditions with the position of the target halo.

Re-run MUSIC with a new parameter file. Change "ref_center" and "ref_extent" using the values from the previous step. Increase the resolution by resetting "levelmax" above 8, which is the previous value in step 1. You should "NOT" change the seed numbers to generate random numbers. Use the same seeds from step 1.

[setup]

boxlength = 50

zstart = 125

levelmin = 8

levelmin_TF = 10

levelmax = 13

padding = 10

overlap = 5

ref_center = xcen/box, ycen/box, zcen/box

ref_extent = xext/box, yext/box, zext/box

align_top = no

baryons = no

use_2LPT = no

use_LLA = no

periodic_TF = yes

[cosmology]

Omega_m = 0.264745

Omega_L = 0.735255

w0 = -1.0

wa = 0.0

Omega_b = 0.0455

H0 = 71

sigma_8 = 0.811

nspec = 0.961

transfer = eisenstein

[random]

seed[7] = 32455

seed[8] = 23456

seed[9] = 34567

seed[10] = 45678

seed[11] = 56789

seed[12] = 67890

seed[13] = 34234

[output]

##Gadget-2 (type=1: high-res particles, type=5: rest)

format = gadget2

filename = ics_gadget_MW_dm.dat

[poisson]

fft_fine = yes

accuracy = 1e-5

pre_smooth = 3

post_smooth = 3

smoother = gs

laplace_order = 6

grad_order = 6