DM only simulations

[Bruno Villasenor Alvarez]

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.

Best regards,

Bruno

[Wise, John H]

Hi Bruno,

Thanks for informing us about these differences between Nyx, Cholla, and
Enzo. I think we all would like to figure out what's causing these
small-scale differences.

I looked at your parameter file, and nearly all of the parameters were
fine. Because all of the methods are the same and this is a unigrid
simulation, my best guess would be that the timestepping is different.
Both Nyx and Enzo restrict the timesteps by particle velocities,
expansion, and acceleration with the same prescription. I saw that you
set the MaxExpansionRate to be 0.015, where both in Nyx and Enzo, the
default is 0.01. Usually the timesteps are restricted in this way early
in the simulation, so any differences might be magnified later on as
halos collapse. What was this factor set to in Nyx & Cholla?

Thanks,
John

> -*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.
>
>
> Best regards,
>
> Bruno
>

> --
> You received this message because you are subscribed to the Google
> Groups "enzo-dev" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to enzo-dev+u...@googlegroups.com
> <mailto:enzo-dev+u...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/enzo-dev/14677bc4-9ae4-43e6-9f40-d90e680eee63%40googlegroups.com
> <https://groups.google.com/d/msgid/enzo-dev/14677bc4-9ae4-43e6-9f40-d90e680eee63%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
John Wise
Associate Professor of Physics
Center for Relativistic Astrophysics, Georgia Tech
http://cosmo.gatech.edu

[Bruno Villasenor Alvarez]

Hi John, 

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.

Best,

Bruno

> an email to enzo...@googlegroups.com
> <mailto:enzo...@googlegroups.com>.

[Bruno Villasenor Alvarez]

Hi,

To follow up on the previous idea of changing MaxExpansionRate, I repeated the Enzo and Cholla simulations using max_delta(a)/a=0.01, I didn't repeat the Nyx simulation. Attached is the power spectrum comparison for the new simulations (power_dm_comparison_expansion0.01.png ), the differences between Enzo and Cholla are still <~25% and the differences between Nyx and Cholla are still <~0.1%.

If there are any other suggestions or ideas please let me know, I'm happy to answer any other questions or do other tests.

Thanks,

Bruno

[Wise, John H]

Hi Bruno,

I'm just forwarding you this email from Mike Norman that was lost
(somehow). Apparently, Enzo uses kick-drift-kick instead of DKD, as
described in the paper. You can look at the source code
(UpdateParticlePositions.C) to see how it works. I'm not that familiar
with that part myself.

Thanks,
John

------------------------------------------------------------------------
*From:* Norman, Michael <mlno...@sdsc.edu>
*Sent:* Thursday, August 29, 2019 5:22 PM
*To:* enzo...@googlegroups.com
*Cc:* James Bordner
*Subject:* RE: [enzo-dev] DM only simulations

Thanks for the detailed note Bruno. It is important to understand this
discrepancy. I seem to recall looking into Enzo’s source and finding
that despite what it says in the method paper, it is actually doing
Kick-Drift-Kick. Could someone confirm that? I could be remembering
wrong. Reason I think I know this is that we compared Enzo-E and Enzo on
a DM only problem, and we do KDK in Enzo-E (I think!).

James, can you please confirm which update scheme is used in both codes?

Mike

> <https://groups.google.com/d/msgid/enzo-dev/14677bc4-9ae4-43e6-9f40-d90e680eee63%40googlegroups.com?utm_medium=email&utm_source=footer

> <https://groups.google.com/d/msgid/enzo-dev/14677bc4-9ae4-43e6-9f40-d90e680eee63%40googlegroups.com?utm_medium=email&utm_source=footer>>.
>
>
> --
> John Wise
> Associate Professor of Physics
> Center for Relativistic Astrophysics, Georgia Tech
> http://cosmo.gatech.edu
>

> --
> You received this message because you are subscribed to the Google
> Groups "enzo-dev" group.
> To unsubscribe from this group and stop receiving emails from it, send

> an email to enzo-dev+u...@googlegroups.com
> <mailto:enzo-dev+u...@googlegroups.com>.

> To view this discussion on the web visit

> https://groups.google.com/d/msgid/enzo-dev/f5ad1661-d621-46f9-9b6a-d951ef5145b2%40googlegroups.com
> <https://groups.google.com/d/msgid/enzo-dev/f5ad1661-d621-46f9-9b6a-d951ef5145b2%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Bruno Villasenor Alvarez]

Hi John and Mike,

Thanks for the answer. I looked at UpdateParticlesPositions.C and indeed it seems like Enzo is doing KDK instead of DKD.

Nyx also uses KDK for the particles integration and for Cholla I've tried both KDK and DKD, actually for the comparisons I sent in the previous messages Cholla was using KDK, if instead I use DKD in Cholla then the differences between Nyx and Cholla go up to ~1% in large scales and ~3% in small scales, logically the differences between Cholla and Enzo also change by 1-3% when using DKD for Cholla, but the ~20% differences in small scales are still there. Just for reference, I attached the power spectrum comparison when using DKD in Cholla.

I looked deeper in the Enzo source code and I saw that in Grid_UpdateParticleVelocity.C there are 3 possible ways of updating the particle velocities:

#ifdef VELOCITY_METHOD1

        /* i) partially time-centered. */

#ifdef VELOCITY_METHOD2

        /* ii) partially backward. */

#ifdef VELOCITY_METHOD3

        /* iii) Semi-implicit way */

The method used was VELOCITY_METHOD_3 which is the most similar to the way I update the velocities in Cholla and the way described in the Nyx paper. 

If the differences are not due to the particles integration scheme then the other two candidates are the CIC interpolation for density and acceleration and the Poisson solver, in theory, Cholla and Enzo are doing the same thing.

I'm happy to run other tests if there are any suggestions.

Thanks,

Bruno