yt-project / yt

Main yt repository
http://yt-project.org
Other
469 stars 280 forks source link

Critical density does not change units using data.critical_density.to("") #5032

Open V-Nathir opened 3 weeks ago

V-Nathir commented 3 weeks ago

Bug report

After loading an output (RAMSESDataset) and asking for critical density in different units, the value does not change. I computed by hand to see the comparative:

Code for reproduction

data = yt_load(path_to_halo685x12, "output_00079") # output at redshift > 1
print(data.current_redshift)
print(data.critical_density.to("Msun/kpccm**3"))
print(data.critical_density.to("Msun/kpc**3"))

Actual outcome

Redshift 1.127658514889251 unyt_quantity(465.91611928, 'Msun/kpc3') unyt_quantity(465.91611928, 'Msun/kpccm3')

The comovil units have not been change.

Expected outcome

Computed by hand (ref: BRYAN and NORMAN 1998, section 2.1)

from astropy.constants import G as G_constant

OMEGA_0   = data.parameters["omega_m"]
OMEGA_L   = data.parameters["omega_l"]
G_constant = data.quan(G_constant.to("kpc**3/(Msun*s**2)"), "kpc**3/(Msun*s**2)") 
redshift_  = data.current_redshift
Ezsqr  = OMEGA_0*(1+redshift_)**3 +OMEGA_L
h =  data.quan(data.parameters["H0"] /100, "km/s/Mpc")
H = 100 *h *Ezsqr**(1/2)
critical_density = 3/8/np.pi/G_constant.to("kpc**3/(Msun*s**2)") * H**2

print(critical_density, critical_density.to("Msun/kpccm**3"))

OUTPUT: unyt_quantity(465.88541137, 'Msun/kpc3'), unyt_quantity(48.36969333, 'Msun/kpccm3')

Version Information

cphyc commented 3 weeks ago

I can replicate this with the following MWE:

import yt
ds = yt.load_sample("output_00080")
print(ds.current_redshift)
print(ds.critical_density.to("Msun/kpc**3"))
print(ds.critical_density.to("Msun/kpccm**3"))

Outputs (incorrectly):

0.14255728632206321
155.7792235477523 Msun/kpc**3
155.7792235477523 Msun/kpccm**3