Hello Enzo developers,
My name is Bruno Villasenor, I'm a graduate student at UCSC working with Brant Robertson and I've been extending Cholla to run cosmological simulations.
Part of my validation process has been to compare with results from other codes, running simulations from the same initial conditions using similar integration schemes. I am writing because we have encountered some differences in small (but resolved) scales in dark matter-only cosmological simulations performed with Cholla and ENZO. While these differences are <~20% in P(k), the differences between results from Cholla and Nyx for the same test are <~0.1%. I suspect my ENZO tests are run incorrectly in some way, and Brant and I were hoping you might share some suggestions or tests that may help us resolve the discrepancy.
I’ve included in this email results from 256^3 cosmological simulations of a 50Mpc/h box, using a uniform 256^3 cells grid for the potential calculation. I have refinement disabled for ENZO.
For Cholla I am computing gravity using the following methods:
- Cloud-In-Cell interpolation of the particle masses onto a uniform grid.
- FFT based Poisson solver using the same Greens function kernel as Equation 61 from Bryan et al. 2014.
- Cloud-In-Cell interpolation for the acceleration at the particle positions from the gravitational field.
- Drift-Kick-Drift scheme for the particle update, same as Equation 62 from Bryan et al. 2014.
The simulations start from exactly the same initial conditions, evolved for the same cosmology, and have used the same ParticleCourantSafetyNumber equal to 0.1.
These are the cosmological parameters I used for the test:
OmegaMatterNow = 0.3111
OmegaLambdaNow = 0.6889
HubbleConstantNow = 0.6766 // in 100 km/s/Mpc
Attached you can find figures for the following comparisons:
-density_dm_comparison: projections ( 25Mpc/h projected) of the dark matter density field from each code and their fractional difference.
-power_dm_comparison: density power spectrum for several redshifts and the fractional difference between code results.
There are significant differences in the power spectrum at small scales. Relative to the ENZO results, I get about 10-25% more power in the Cholla simulation, and the differences in the density projection go up to 50%. Performing the same comparison between Cholla and Nyx, the differences are much smaller, less than 0.1% in the power spectrum and less than 0.5% in the density projection.
For each code (Cholla, Nyx, and ENZO), the density field is computed from the particle positions using the same CIC method for each simulation, and the density power spectrum is computed from the density field in the same manner for each simulation.
There is probably something wrong with the way I am running the Enzo simulations, attached you can find the ENZO parameter file I used (parameter_file_dm_only.txt ). The ENZO version I used was v2.5.
Could you please let me know if there is anything wrong with the parameters I used for the Enzo simulation? I would also appreciate any insights you may have on how these differences may arise.
Also, if you have any reference dark matter only ENZO simulations that I could compare to and that you are willing to share, that would be great.
I am happy to answer any question you might have or share any other information.
I appreciate any suggestions and thank you for your time.
Thanks for the answer.
Yes, at the beginning of the simulation the timestep is limited by expansion ( max_delta(a)/a ). For my tests I used a value of 0.015 for both Enzo and Cholla and for Nyx I'm not completely sure which parameter controls this condition, these are the timestep parameters I used for Nyx:
# TIME STEP CONTROL
particles.cfl = 0.1 # 'cfl' for particles
nyx.cfl = 0.1 # cfl number for hyperbolic system
nyx.init_shrink = 1.0 # scale back initial timestep
nyx.change_max = 1.1 # factor by which timestep can change
nyx.dt_cutoff = 5.e-20 # level 0 timestep below which we halt
From the Nyx user's guide:
nyx.init_shrink : factor by which to shrink the initial time step. (default 1.0)
nyx.change_max: factor by which the time step can grow in subsequent steps.(default 1.1)
I will repeat the Enzo simulation using max_delta(a)/a = 0.01 and see if there is a significant change in the results, but I doubt this change will make a significant difference.
Thanks for the suggestion and I appreciate any other ideas.