Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
V-Nathir opened this issue Oct 28, 2024 · 1 comment
Open
Labels

Comments

@V-Nathir
Copy link
Contributor

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

  • Operating System: Linux
  • Python Version: 3.11.5
  • yt version: 4.4.0
  • Other Libraries (if applicable):
@cphyc cphyc added the bug label Oct 28, 2024
@cphyc
Copy link
Member

cphyc commented Oct 28, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants