Correct m200 mean halo masses #1
Conversation
Cleaned up code to calculate cosmological properties, now more general. Also added extra cosmological fields to config file.
Fixed bug when calculated matter density. Updated overdensity thresholds when calculating spherical overdensity masses.
Missing call to calculation of background density when reading gadget/hdf/ramses. Swift interface may need update.
| //if opt.virlevel<0, then use virial overdensity based on Bryan and Norman 1997 virialization level is given by | ||
| if (opt.virlevel<0) opt.virlevel=opt.virBN98; | ||
|
|
||
| //normally Hubbleflow=lvscale*Hubble but we only care about peculiar velocities | ||
| //ignore hubble flow | ||
| Hubbleflow=0.; | ||
| cout<<"Cosmology (h,Omega_m,Omega_cdm,Omega_b,Omega_L) = ("<< opt.h<<","<<opt.Omega_m<<","<<opt.Omega_cdm<<","<<opt.Omega_b<<","<<opt.Omega_Lambda<<")"<<endl; | ||
| PrintCosmology(opt); |
There was a problem hiding this comment.
The same calculations done here are done in all data input routines, including the swift interface: Omega-k, critical density, background density, and VirBN98. All these are also done in CalcCosmoParams, which is defined in substructureproperties.cxx but neither declared nor used. On top of that there is more common code to all the data input routines: the calculation of the adjustment factor from the comove option, the adjustment of the virlevel option when it's negative, and the printing of the cosmological parameters.
Taking this into account, I'd take this opportunity to move these three additional operations into CalcCosmoParams (which then wouldn't take an extra Double_t argument, as it would be calculated internally), and then add a call to CalcCosmoParams from each data input routine replacing the current common contents.
As a result, the calls to GetHubble inside each data input routine are not necessary anymore, and the Hubble variable (and probably the aadjust variable) can be removed.
There was a problem hiding this comment.
Just to clarify: the code displayed above doesn't give the full context for the comment, which is better seen in the "Files changed" tab of the pull request.
| Double_t m200mval=log(opt.rhobg*200.0); | ||
| Double_t mBN98val=log(opt.virBN98*opt.rhobg); | ||
| Double_t mBN98val=log(opt.virBN98*opt.rhocrit); |
There was a problem hiding this comment.
Is this intended? I don't know the physical meaning behind rhocrit and rhobg, but they are clearly different values which seemed to have been computed differently in different parts of the code. I see that this pull is clearly trying to fix that, but this particular change (and that in line 2216 below) looks more like yet a different problem being solved, so I just wanted to make sure it was effectively intended.
There was a problem hiding this comment.
Thanks for the thorough review! I agree with moving the duplicated code into a separate function. I guess we'll have to discuss how far we're prepared to diverge from Pascal's repository.
rhobg is the mean mass density of the universe and rhocrit is the critical density above which the universe would collapse in on itself. There is also a density parameter omega which is defined as omega=rhobg/rhocrit. All three quantities vary with time. The variable omega_m in the code is the density parameter computed at the present day (redshift z=0 or expansion factor a=1) but the code was using it to compute rhocrit from rhobg at earlier times, which was wrong. The changes at lines 410 and 414 are part of the fix for this - instead of computing rhocrit inline when we need it, we compute it once and store it.
The change to mBN98val does look like a separate bug fix. I believe the new formula is correct (and the output agrees with the EAGLE code in my tests) but I'll check the BN98 paper to make sure.
There was a problem hiding this comment.
The quantity opt.virBN98 is delta_c in formula 6 in Bryan and Norman (1998). The text states that it should be multiplied by the critical density to get the mean density within the virial radius, so the changes to mBN98val are fixing a bug.
| @@ -2213,10 +2213,10 @@ void GetInclusiveMasses(Options &opt, const Int_t nbodies, Particle *Part, Int_t | |||
| Double_t change=MAXVALUE,tol=1e-2; | |||
| Int_t ii,icmv,numinvir,num200c,num200m; | |||
| Double_t virval=log(opt.virlevel*opt.rhobg); | |||
| Double_t mBN98val=log(opt.virBN98*opt.rhobg); | |||
| Double_t m200val=log(opt.rhobg/opt.Omega_m*200.0); | |||
| Double_t mBN98val=log(opt.virBN98*opt.rhocrit); | |||
| @@ -237,6 +237,7 @@ void SetVelociraptorSimulationState(cosmoinfo c, siminfo s) | |||
| libvelociraptorOpt.ellxscale*=libvelociraptorOpt.a; | |||
| libvelociraptorOpt.uinfo.eps*=libvelociraptorOpt.a; | |||
|
|
|||
| /* | |||
It might be best not to merge this immediately but I'm making a pull request to document this fix.
Here I've tried to cherry pick the changes from Pascal's expandedproperties branch which are needed to fix the incorrect m200 halo masses which Velociraptor was generating. The problem was that the density threshold used for the spherical mass estimates had an incorrect dependence on redshift so it diverged from the correct value as redshift increased. The calculation of the threshold was duplicated over the various read routines and the SWIFT interface. In the fixed version there's a single function to do the calculation.
So far I've checked that this generates the same m200 mean and critical mass functions as the EAGLE code at low and high redshift in the 25Mpc DMONLY box. I haven't tried running through the swift interface yet.