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

BUG: RAMSES star_age and particle_birth_time are swapped #3861

Closed
lconaboy opened this issue Mar 25, 2022 · 3 comments · Fixed by #3864
Closed

BUG: RAMSES star_age and particle_birth_time are swapped #3861

lconaboy opened this issue Mar 25, 2022 · 3 comments · Fixed by #3864

Comments

@lconaboy
Copy link
Contributor

lconaboy commented Mar 25, 2022

Bug report

Bug summary

Looking at the yt docs (https://yt-project.org/doc/examining/loading_data.html#ramses-data), ('star', particle_birth_time') should give me the formation time relative to the big bang and ('star', 'star_age') should give me the number of years a star particle's been around. Instead, it seems the other way around -- I'm getting formation times of 0 Myr and ages equal to the output's current time (see example below)., where I'd expect the smallest birth time to be ~171 Myr and the smallest age to be ~0 Myr.

The ages from ramses/utils/f90/getstarlist.f90 match particle_birth_time as currently calculated in yt, mostly to within +/- 1 Myr, with some tail. Indeed, the return (tsim - t) * t_scale line in convert_ramses_ages is how ages are calculated in getstarlist.f90 (see line 336 in https://bitbucket.org/rteyssie/ramses/src/master/utils/f90/getstarlist.f90).

So, in yt particle_birth_time = formation_time = time_simu - t should actually be star_age , and the function star_age returns data.ds.current_time - formation_time = time_tot + time_simu - time_simu + t = time_tot + t = particle_birth_time. These time_simu and time_tot parameters are calculated from the friedman function from ramses, which is reproduced in yt. t is the birth time of the star. I think the fix then is something like: convert_rames_ages in frontends/ramses/io.py instead of returning

(tsim - t) * t_scale

should return

data.ds.current_time - ((tsim - t) * t_scale)

or similar.

Code for reproduction
(Sorry, haven't had time to test with the yt data files, can do if you need me to)

import yt
import numpy as np
from yt.data_objects.particle_filters import add_particle_filter

print('yt version:', yt.__version__)

bbox = [[0.4, 0.4, 0.4], [0.6, 0.6, 0.6]]
ds = yt.load('./output_00031', bbox=bbox)
c = ds.arr([50., 50., 50.], 'Mpccm/h')
r = ds.arr(5., 'Mpccm/h')

# Get a sphere
s = ds.sphere(c, r)

# Print current time
print(f'ds.current_time: {ds.current_time.in_units("Myr"):.6f}')
print(f'ds.cosmology.t_from_z(ds.current_redshift): {ds.cosmology.t_from_z(ds.current_redshift).in_units("Myr"):.6f}')

# Check they all have the right family value
assert((s['star', 'particle_family'] == 2).all())

# Try native particle_birth_time from yt-determined stars
t = s['star', 'particle_birth_time'].in_units('Myr')
a = s['star', 'star_age'].in_units('Myr')
print('-- using star family')
print(f'(birth) min: {t.min():.3f}, min nonzero: {t[t>0.].min():.3f}, max: {t.max():.3f}')
print(f'(age)   min: {a.min():.3f}, min nonzero: {a[a>0.].min():.3f}, max: {a.max():.3f}')
print(f'(b+a)   allclose(a+t, ds.current_time): {np.allclose(np.array(a+t), np.array(ds.current_time.in_units("Myr")))}')


# Try explicitly filtering particles a la docs
def filtered_stars(pfilter, data):
    filter = data['io', 'conformal_birth_time'] != 0
    return filter

add_particle_filter('filtered_stars',
                    function=filtered_stars,
                    filtered_type='io',
                    requires=['conformal_birth_time'] )
ds.add_particle_filter('filtered_stars')

# Get a sphere
s = ds.sphere(c, r)

t = s['filtered_stars', 'particle_birth_time'].in_units('Myr')
a = s['filtered_stars', 'star_age'].in_units('Myr')
print('-- using particle filter')
print(f'(birth) min: {t.min():.3f}, min nonzero: {t[t>0.].min():.3f}, max: {t.max():.3f}')
print(f'(age)   min: {a.min():.3f}, min nonzero: {a[a>0.].min():.3f}, max: {a.max():.3f}')
print(f'(b+a)   allclose(a+t, ds.current_time): {np.allclose(np.array(a+t), np.array(ds.current_time.in_units("Myr")))}')

Actual outcome

yt version: 4.0.2
yt : [INFO     ] 2022-03-24 17:21:00,117 Parameters: current_time              = 4.223387956490154
yt : [INFO     ] 2022-03-24 17:21:00,118 Parameters: domain_dimensions         = [256 256 256]
yt : [INFO     ] 2022-03-24 17:21:00,119 Parameters: domain_left_edge          = [0. 0. 0.]
yt : [INFO     ] 2022-03-24 17:21:00,120 Parameters: domain_right_edge         = [1. 1. 1.]
yt : [INFO     ] 2022-03-24 17:21:00,120 Parameters: cosmological_simulation   = 1
yt : [INFO     ] 2022-03-24 17:21:00,120 Parameters: current_redshift          = 11.690083646079087
yt : [INFO     ] 2022-03-24 17:21:00,120 Parameters: omega_lambda              = 0.685999989509583
yt : [INFO     ] 2022-03-24 17:21:00,120 Parameters: omega_matter              = 0.314000010490417
yt : [INFO     ] 2022-03-24 17:21:00,120 Parameters: omega_radiation           = 0.0
yt : [INFO     ] 2022-03-24 17:21:00,121 Parameters: hubble_constant           = 0.673000030517578
yt : [INFO     ] 2022-03-24 17:22:03,834 Adding particle_type: DM
yt : [INFO     ] 2022-03-24 17:22:03,849 Adding particle_type: star
yt : [INFO     ] 2022-03-24 17:22:03,864 Adding particle_type: cloud
yt : [INFO     ] 2022-03-24 17:22:03,880 Adding particle_type: dust
yt : [INFO     ] 2022-03-24 17:22:03,896 Adding particle_type: star_tracer
yt : [INFO     ] 2022-03-24 17:22:03,912 Adding particle_type: cloud_tracer
yt : [INFO     ] 2022-03-24 17:22:03,928 Adding particle_type: dust_tracer
yt : [INFO     ] 2022-03-24 17:22:03,944 Adding particle_type: gas_tracer
ds.current_time: 381.033802 Myr
ds.cosmology.t_from_z(ds.current_redshift): 382.298844 Myr
-- using star family
(birth) min: 0.000 Myr, min nonzero: 0.002 Myr, max: 210.223 Myr
(age)   min: 170.811 Myr, min nonzero: 170.811 Myr, max: 381.034 Myr
(b+a)   allclose(a+t, ds.current_time): True
-- using particle filter
(birth) min: 0.000 Myr, min nonzero: 0.002 Myr, max: 210.223 Myr
(age)   min: 170.811 Myr, min nonzero: 170.811 Myr, max: 381.034 Myr
(b+a)   allclose(a+t, ds.current_time): True

Expected outcome
Would expect birth times >~ 171 Myr (from looking at the simulation output logs) and stellar ages from 0 (youngest) to 381 - 171 = 210 Myr (current time - oldest).

Version Information

  • Operating System: Ubuntu 18.04
  • Python Version: 3.8.5
  • yt version: 4.0.2
  • Other Libraries (if applicable): numpy

python from anaconda distribution, yt installed through pip.

@welcome
Copy link

welcome bot commented Mar 25, 2022

Hi, and welcome to yt! Thanks for opening your first issue. We have an issue template that helps us to gather relevant information to help diagnosing and fixing the issue.

@matthewturk
Copy link
Member

Thanks for opening this, @lconaboy ! @cphyc , do you think you can take this on?

@cphyc
Copy link
Member

cphyc commented Mar 28, 2022

Hi @lconaboy , could you tell me if #3864 is indeed working as you expect?

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

Successfully merging a pull request may close this issue.

3 participants