Cosmological Hydrodynamical Codes width buildin or external Shock detection comparison

The goal of the Chewbacca project is to investigate how much the prediction of the shock front distributions (in strength and volume filling) by cosmological, hydrodynamic Al simulations depend on the underlying numerical methods. Therefore we will investigate 3 different cosmological hydro codes ( TVD-code, Gadget-2 and ENZO). We want to understand how much they differ already in their predictions of the low density, low temperature gas, as well as how the different ability of resolving low/high density regions and the structure within them influences the shock distribution. We also want to test internal versus different external shock detection methods (e.g. based on Density/Temperature versus velocity jumps).
Name explanation Chewbacca.
Further webpages related to this Project:

Initial conditions

To follow a resonably resolved cosmological box with a massive cluster we setup a box of size 100 Mpc/h with a intentionally large value of Sig8: There is a set of initial conditions for 64³, 128³, 256³ and 512³ dark matter particles, inclusive some description on the data format available at http://canopus.cnu.ac.kr/shocks/.

Here you can find some theoretical predictions of (cumulative) massfunctions for the standard LCDM model and the model addapte here. First column contains log10(M) in units of Msol/h, second column has the prediction of cumulative mass function log10(N) using Press & Schechter in units of 1/(Mpc/h)³ and third column the same using Sheth & Tormen.

Here is an idl program to create Gadget initial conditions from the above dark matter particle positions, setup_ic.pro.

The ICs for the gadget runs and the simulation data can be found on clx at CINECA under /clx/userhpe/hpedeue7/Data/ShockComp/

Performance

Gadget-2(§) (dm only run)
Run Tcpu [h] Nproc Ntstep Mdm [Msol/h] Rsoft [kpc/h]
64³ 1.2 4 1153 28.6e10 31.0
128³ 19.1 8 1985 3.57e10 15.5
256³ 258.1 16 3035 0.45e10 7.75
512³ 6544.49 64 4943 5.6e8 3.775
Gadget-2(§)
Run Tcpu [h] Nproc Ntstep Mdm/Mgas [Msol/h] Rsoft [kpc/h]
64³ 5.2 4 3316 24.0e10 / 4.55e10 31.0
128³ 144.4 8 9646 3.0e10 / 0.57e10 15.75
256³ 1893.53 16 25896 3.76e9 / 7.11e8 7.875
Gadget-2(§) (box like runs)
Run Tcpu [h] Nproc Ntstep Mdm/Mgas [Msol/h] Rsoft [kpc/h]
64³ 4.9 4 2321 24.0e10 / 4.55e10 31.0
128³ 109.8 8 6346 3.0e10 / 0.57e10 15.75
256³ 1484.57 32 17205 3.76e9 / 7.11e8 7.875
ENZO(+)
Run Tcpu [h] Nproc Ntstep Mdm [Msol/h] Grid [kpc/h]
64³ 0.52 16 241 1562.5
128³ 7.7 16 304 781.25
TVD(*)
Run Tcpu [h] Nproc Ntstep Mdm [Msol/h] Grid [kpc/h]
(512-256)³ ~200 2 468 195.312

Checks on DM particles

Particle positions

Shown are the dark matter particle positions focusing on a 30Mpc/h section of the simulation. Black are the particles from the 256³ run, blue the one from the 128³ run and red the ones from the 64³ run.
Gadget-2 Gadget-2 (grid like IC) TVD

Shown are the dark matter particle positions focusing on a 30Mpc/h section of the simulation. Black are the particles from the TVD run, blue are the from the Gadget-2 and red are from the Gadget-2 run with the grid like IC. From left to right are the 64³, 128³ and 256³ runs. Additional a 15Mpc/h zoom for the 256³ run is shown. comments ...
64³: 128³: 256³: 256³ (zoom):

Massfunctions

Shown are the cumulative mass functions for the different runs. They are always compared to the prediction from using Sheth & Tormen for the cosmology addapted (thick yellow line). Additinal the thin yellow lines mark the plus/minus sqrt(N) to get an proxy for the statistical errors.
Results are based on spherical overdensity halo finder (using the overdensity from a spherical collaps model) are shown. To apply the halofinder to the TVD simulations, the results are converted into an artificial gadget format using ryu2snap.pro, then exactly the same post processing is used.
To compare the simulations, once has to assign a formal resolution to it. Although the massfunctions seems to be resoved down to the smalest haloes, build from less than 20 particles, Power et al. 2003 showed that to converge on resolving the structure within haloes, once need much more particles. From Fig. 14 once can read that 500-1000 particles will be needed to resolve Rvir at the density contrast of interest for clusters (1e4-1e5). Therefore i define the resolution of the gadget runs as 500 times the mass of the dm particles. For the TVD runs i assume that a radial profile needs 5 cells to be resolved. I then use the predicted Rvir vs. Mvir relation from the highest dm run to assign a massresolution according to the Mvir corresponding to Rvir = 5 cells. Although these values look quite artificial and conservative in the mass functions or the Rvir vs. Mvir plots (See small vertical lines), they seem to predict the convergence of the baryon fractions across all runs very well. Additional to the mass-functions also the Mvir vs. Rvir relations are shown for the gadget runs. The lower points are the fraction of the measured Rvir and the best fit to the highest resolution run.
Validating the ICs (using Gadget)
Gadget-2 (dm only run) Gadget-2 (SPH run, mass is the total = dark matter + baryonic mass)
Massfunctions (individual codes)
Gadget-2 (SPH run, grid like IC) TVD (Results are based on converting the grid to artificial SPH particles, ryu2snap.pro) ENZO(Results are based on converting the grid to artificial SPH particles, enzo2snap.pro)
Gadget-2 vs. TVD Gadget-2 vs. ENZO TVD vs. ENZO

Checks on gas distribution

Baryon fraction

Shown are the baryon fraction for the different runs. The vertical lines are indicating the mass resolution as defined previously. Interesting to note that they very well represent the point at which the results seems to be converged in respect to the higher resolution runs.
In agreement with the definition of resolution we used before, the TVD 512-256 runs seems to be bracketed by the Gadget 64/128 runs. It has to be noticed that the baryon fraction between the TVD and gadget simulations seems to be consistent, exept the two most massive systems, where the gadget runs predict significantly smaller baryon fractions than the TVD, which confirms known trends. However the TVD runs also predict baryon fractions significantly smaller than unity, which is is supprosingly consistent with the gadget results and and in contrast with other findings in literature, where grid codes tend to predict baryon fractions larger than unity.
Baryon fraction
Gadget-2 Gadget-2 (box like ICs) TVD ENZO
Gadget vs.TVD Gadget vs.ENZO TVD vs.ENZO

Density and Temperature distribution

Shown are the filling factors of temperature and density in respect to the total volume/mass of the simulation. First all simulations along with their resolution are shown, then the simulations are comapred.
Gadget-2
rho Volume rho Mass T Volume T Mass
Gadget-2 (box like ICs)
rho Volume rho Mass T Volume T Mass
TVD
rho Volume rho Mass (box like ICs) T Volume T Mass
ENZO
rho Volume rho Mass (box like ICs) T Volume T Mass
Comparison
rho Volume rho Mass (box like ICs) T Volume T Mass

2D projected Maps

Shown are maps of projections through the hole 100 Mpc/h box for different resolutions. From left to right 64³,128³,256³ and 512³. A small idl program to produce such maps from the TVD data is ryu2map.pro. The idl program to produce the maps from the Enzo data is enzo2map.pro.
The fits files with the maps can be downloaded from the Maps subdirectory.
Projected gas density
gas b0
64³ 128³ 256³ 64³ 128³ 256³
TVD ENZO
(128-64)³ (256-128)³ (512-256)³ 128³ 256³ 512³
Projected gas temperature (mass weighted)
gas b0
64³ 128³ 256³ 64³ 128³ 256³
TVD ENZO
(128-64)³ (256-128)³ (512-256)³ 128³ 256³ 512³

Checks on clusters

Gadget:

Listed are the properties of the two most massive clusters for the runs with the different resolutions.
Cluster A
Run X [Mpc/h] Y [Mpc/h] Z [Mpc/h] Mvir [Msol/h] Rvir [kpc/h]
dm 64³ 60.6652 75.5600 88.17896 2.611e15 2890.30
128³ 60.6522 75.8406 87.8016 2.695e15 2921.52
256³ 60.5536 75.8819 87.7350 2.747e15 2940.53
gas 64³ 60.6398 75.6028 87.9461 2.587e15 2881.69
128³ 60.6273 75.7304 87.7092 2.699e15 2922.69
b0 64³ 61.6577 75.6777 88.9890 1.127e15 2290.93
128³ 61.5856 75.6988 88.8988 1.3375e15 2312.70
256³ 61.4941 75.6521 88.8680 1.359e15 2325.42
Cluster B
Run X [Mpc/h] Y [Mpc/h] Z [Mpc/h] Mvir [Msol/h] Rvir [kpc/h]
dm 64³ 92.4939 47.1792 48.0768 1.687e15 2498.80
128³ 92.6371 47.1574 48.0683 1.586e15 2448.27
256³ 92.6508 47.1263 48.0983 1.606e15 2458.80
gas 64³ 92.6904 47.2706 48.0845 1.631e15 2471.17
128³ 92.5842 47.1623 48.0751 1.605e15 2457.75
b0 64³ 92.7272 47.6547 47.9963 1.546e15 2427.32
128³ 92.8785 47.7715 48.0926 1.623e15 2466.80
256³ 93.1952 47.8979 48.0630 1.639e15 2475.01

Shown are projected temperature maps (mass weighted) of a 8Mpc large region centered on the two main clusters for different resolutions. From left to right 64³,128³,256³ and 512³.
Cluster A
gas b0
64³ 128³ 256³ 64³ 128³ 256³
Cluster B
gas b0
64³ 128³ 256³ 64³ 128³ 256³

Shown are threedimensional profiles of the cluster A and B for the different runs (colorcoded). The gray lines mark the resolution limits used before. Theirby the mass resolution of the SPH runs is translated into a spacial resultion using the Mvir vesus Rvir relation, analogous as we did for the grid codes in other directions. The dashed lines on the botom mark the radius enclosing 500 gas particles for the SPH runs. The extension .b0 labels the runs with grid like initial conditions.
Gadget
Cluster A Cluster B
density temperature entropy density temperature entropy
TVD
Cluster A Cluster B
density temperature entropy density temperature entropy
ENZO
Cluster A Cluster B
density temperature entropy density temperature entropy
Gadget vs. TVD
Cluster A Cluster B
density temperature entropy density temperature entropy
Gadget vs. ENZO
Cluster A Cluster B
density temperature entropy density temperature entropy
TVD vs. ENZO
Cluster A Cluster B
density temperature entropy density temperature entropy

Checks on entropy profile

Shown are the entropry profoiles, where for the gadget runs the most dens parts in each radial bin are excluded. It looks like that the differences theerby gets smaller (specially when focussing on the entropy plateau at 2-4 Rvir. Also the profile with increasing resolution also evolves towards the results of the grid codes.
Cluster B
Volume fraction for different treshholds Entropy profiele using different tresholds Convergence test for low treshold Comparison using low treshold

Checks on Machnumbers

Gadget:

Shown are the differential machnumber statistics of the full box without deviding into internal and external shocks (like upper pannels of fig.7 in Pfrommer et al. 2006).
Differential machnumber distributions
Gadget-2 Gadget-2 (box like ICs)