Tutorial
In this tutorial, we demonstrate the usage of the PhaseSpaceDTFE package to estimate the density and velocity fields from a GADGET-4 simulation.
We start by importing the relevant libraries and loading the data:
using JLD2, Plots, HDF5, ProgressMeter, PhaseSpaceDTFE
## set up simulation box
Ni = 64
L = 100.
sim_box = SimBox(L, Ni) ## need this for estimator creation
depth = 5 # depth of estimator search tree
## load data
function load_data(file)
fid = h5open(file, "r")
pos = convert.(Float64, read(fid["PartType1"]["Coordinates"]))
vel = convert.(Float64, read(fid["PartType1"]["Velocities"]))
ids = read(fid["PartType1"]["ParticleIDs"])
time = read_attribute(fid["Header"], "Time")
close(fid)
ordering = sortperm(ids)
return (copy(pos[:,ordering]'), copy(vel[:,ordering]'), time)
end
function load_mass(file)
f = h5open(file, "r")
read_attribute(f["Header"], "MassTable")[2] # particle type 1
end
m = load_mass("../../test/data/snapshot_000.hdf5")
(coords_q, _, _) = load_data("../../test/data/snapshot_000.hdf5")
(coords_x, vels, _) = load_data("../../test/data/snapshot_002.hdf5")([95.94856262207031 98.76901245117188 99.14542388916016; 99.86626434326172 97.77499389648438 97.80353546142578; … ; 94.38851928710938 97.52845764160156 0.3155221939086914; 94.68167114257812 97.48580169677734 99.97648620605469], [-292.1794128417969 -76.47138977050781 28.625730514526367; -60.02801513671875 -113.5673599243164 -126.30517578125; … ; -395.4625244140625 163.2576141357422 155.2444610595703; -418.26104736328125 123.40419006347656 328.3358459472656], 0.9999999999999997)The particle coordinates and velocities are Float64 matrices of size (N, 3). The particle mass m is a single Float64 or a matrix of size (N, 3) for individual particle masses.
Prequel: Delaunay Tesselation Field Estimator
Before going though through PS-DTFE method, we demonstrate the traditional DTFE method by calling the PS-DTFE code only on the final (Eulerian) particle positions. For details, see examples below.
## construct estimator
ps_dtfe = PS_DTFE_periodic(coords_x, coords_x, vels, m, depth, sim_box)
## This is equivalent to the DTFE_periodic routine
# ps_dtfe = DTFE_periodic(coords_x, vels, m, depth, sim_box)
## evaluate density field
Range = 0:2.0:100.
density_field = [PhaseSpaceDTFE.density([L/2., y, z], ps_dtfe) for y in Range, z in Range]
heatmap(Range, Range, log10.(density_field), aspect_ratio=:equal, xlims=(0, L), ylims=(0, L), c=:grays, xlabel="[Mpc]", ylabel="[Mpc]")
Phase-Space Delaunay Tessellation Field Estimator — basic implementation
We now demonstrate the use of the PS-DTFE method with the basic implementation that is suitable to simulations up to size 128^3 particles.
The first step is the construction of the estimator object from the initial (Lagrangian) and final (Eulerian) positions, coords_q and coords_x. This is done only once as a pre-computation step.
## construct estimator
ps_dtfe = PS_DTFE_periodic(coords_q, coords_x, vels, m, depth, sim_box)
## if want to ignore velocities
# ps_dtfe = PS_DTFE_periodic(coords_q, coords_x, m, depth, box)
## for further use without pre-computation, consider saving the estimator to file
# save("ps_dtfe.jld2, "ps-dtfe", ps_dtfe)
nothingThe argument depth specifies the simplex search tree depth in the estimator. Higher tree depths result in faster field evaluations, but require longer construction times. We recommend to start with depth=5 and increase this if required for high-resolution density fields.
The construction time should be of order 1-2 minutes for a 64^3 simulation at depth=7, or a 128^3 simulation at depth=5 on a modern computer.
We now evaluate the density field with the density() function:
# evaluate density field
density_field = [PhaseSpaceDTFE.density([L/2., y, z], ps_dtfe) for y in Range, z in Range]
heatmap(Range, Range, log10.(density_field), aspect_ratio=:equal, xlims=(0, L), ylims=(0, L), c=:grays, xlabel="[Mpc]", ylabel="[Mpc]")
The corresponding number of streams is evaluated with the numberOfStreams() function:
nstreams_field = [numberOfStreams([L/2., y, z], ps_dtfe) for y in Range, z in Range]
heatmap(Range, Range, nstreams_field, aspect_ratio=:equal, xlims=(0, L), ylims=(0, L), clim=(1, 7), xlabel="[Mpc]", ylabel="[Mpc]")
Similarly, the velocity field is evaluated with the velocity()-function:
velocity_field = [velocity([L/2., y, z], ps_dtfe) for y in Range, z in Range]51×51 Matrix{Matrix{Float64}}:
[42.8752 -74.5904 18.2795] … [42.8752 -74.5904 18.2795]
[27.7579 -59.2485 -5.62945] [27.7579 -59.2485 -5.62945]
[23.8106 -157.228 -72.9788] [23.8106 -157.228 -72.9788]
[12.0603 -122.481 -83.4588] [12.0603 -122.481 -83.4588]
[-15.8904 -94.7159 -69.8274] [-13.0526 -92.4946 -75.5898]
[-17.0934 -138.316 -53.3808] … [-17.0934 -138.316 -53.3808]
[6.51849 -107.328 -53.3064] [6.51849 -107.328 -53.3064]
[16.0707 -53.784 -48.2457] [16.0707 -53.784 -48.2457]
[17.5343 5.09039 -55.7492] [24.5043 0.472117 -72.4606]
[21.5663 61.8256 -63.4809] [28.5363 57.2074 -80.1923]
⋮ ⋱ ⋮
[245.46 -306.869 77.1011] [233.843 -313.44 70.9182]
[296.209 -240.779 131.363] [158.601 -410.851 143.251]
[158.889 -345.196 51.9522] [173.018 -389.524 10.5611]
[-3400.49 -1613.24 -1424.04; 151.918 -281.171 -36.3125; 479.391 -222.284 -23.6498] … [196.684 -254.992 84.1738; -559.954 -331.328 204.745; -1327.65 -672.563 -500.751]
[82.1659 -179.572 74.1128] [236.482 -211.028 10.5177]
[95.9255 -156.069 66.1979] [106.702 -153.601 49.6197]
[64.8472 -144.103 62.0064] [61.398 -142.112 54.7716]
[39.2382 -106.284 42.5978] [28.8622 -126.336 44.6443]
[42.8752 -74.5904 18.2795] … [42.8752 -74.5904 18.2795]In multistream regions, the velocity()-function returns the velocities of the individual streams (or NaN if single_stream=true is set in the function). To obtain the stream-mass weighted summation of the velocities, call the velocitySum()-function, which reduces to the velocity()-function in single-stream regions.
velocity_field = [velocitySum([L/2., y, z], ps_dtfe) for y in Range, z in Range]The Phase-Space Delaunay Tessellation Field Estimator — subbox implementation
We now demonstrate the use of the PS-DTFE method for simulations with more than 128^3 particles.
It is not feasible to directly apply the basic PS-DTFE implementation to high-resolution simulations, as the construction of the estimator's simplex search tree would require immense working memory (> 100 GB for 256^3 particles). To circumvent this, the subbox routine internally divides the simulation box into smaller subboxes, constructs an estimator for each of these and writes the estimator to file. The user constructs the ps_dtfe_sb object holding the subbox references as follows:
## construct estimators with velocities
ps_dtfe_sb = ps_dtfe_subbox(coords_q, coords_x, vels, m, depth, sim_box; N_target=32)
## construct estimator without velocities
# ps_dtfe_sb = ps_dtfe_subbox(coords_q, coords_x, m, depth, sim_box; N_target=32)
## it is recommended to save the estimator object (holding the subbox references) for further use
save("ps_dtfe_sb.jld2", "ps-dtfe-sb", ps_dtfe_sb)
ps_dtfe_sb = load("ps_dtfe_sb.jld2")["ps-dtfe-sb"]
nothingCalculating 8 subbox estimators with 32^3 particles each (total simulation 64^3 particles).
subbox i, j, k = [0, 0, 0]
subbox i, j, k = [1, 0, 0]
Progress: 25%|██████████▎ | ETA: 0:00:14subbox i, j, k = [0, 1, 0]
Progress: 38%|███████████████▍ | ETA: 0:00:13subbox i, j, k = [1, 1, 0]
Progress: 50%|████████████████████▌ | ETA: 0:00:10subbox i, j, k = [0, 0, 1]
Progress: 62%|█████████████████████████▋ | ETA: 0:00:07subbox i, j, k = [1, 0, 1]
Progress: 75%|██████████████████████████████▊ | ETA: 0:00:05subbox i, j, k = [0, 1, 1]
Progress: 88%|███████████████████████████████████▉ | ETA: 0:00:02subbox i, j, k = [1, 1, 1]
Progress: 100%|█████████████████████████████████████████| Time: 0:00:18The keyword argument N_target specifies the particle number (N_target^3) of the subboxes. We recommend to use the default value N_target=128.
For a 256^3 simulation with 8 subboxes of size N_target=128 at depth=5-7, the construction time should be of order 10-30 minutes. The estimator objects will require about 20-50 GB of storage space, which can be deleted after the field evaluations (see below).
For internal efficiency, the density field is evaluated by directly passing on the list of coordinates to the density_subbox()-function:
coords_arr = [[L/2., y, z] for y in Range, z in Range]
density_field = density_subbox(coords_arr, ps_dtfe_sb)
heatmap(Range, Range, log10.(density_field), aspect_ratio=:equal, xlims=(0, L), ylims=(0, L), c=:grays, xlabel="[Mpc]", ylabel="[Mpc]")
The number of streams follows analogously with numberOfStreams_subbox()-function:
nstreams_field = numberOfStreams_subbox(coords_arr, ps_dtfe_sb)
heatmap(Range, Range, nstreams_field, aspect_ratio=:equal, xlims=(0, L), ylims=(0, L), clim=(1, 7), xlabel="[Mpc]", ylabel="[Mpc]") Finally, the velocities are evaluated with the velocity()- or velocitySum()-function:
velocitySum_field = velocitySum_subbox(coords_arr, ps_dtfe_sb)We clear the temporary files here. The user might wish to consider storing the estimators for further use.
rm("ps_dtfe", recursive=true)
rm("ps_dtfe_sb.jld2")