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