Simplify handling of conformal birth times in RAMSES

[cphyc]

RAMSES ships with a utility to convert the particle creation times from conformal units into Gyr. This creates one file named birthXXXX.outYYYYY, where XXXXX is the number of the output and YYYYY is the index of the computation domain.

This PR detects those files and reads the birth times from there rather than computes them from the conformal time.

• New features are documented, with docstrings and narrative docs
• Adds a test for any bugs fixed. Adds tests for new features.

[matthewturk]

In general, I'm inclined to accept this. Do you have an idea what's up with the tests?

[neutrinoceros]

Do you have an idea what's up with the tests?

Logs are lost, let's try again: @yt-fido test this please

[cphyc]

Before we merge this, I'd like to backport the changes I made in pynbody/pynbody#740.

[neutrinoceros]

@cphyc could you select appropriate labels for this ? I don't know wether this is a pure refactor or an enhancement.

[cphyc]

Well, this is both a refactor and an enhancement :)

[cphyc]

Welp, it looks like the following fields appear broken with RAMSES:

• particle_birth_time,
• age.
  and the whole thing is a complete mess. Not sure how to proceed any further. This may be related to BUG: some RAMSES fields return different values depending on order of operations #4880, but I'm unsure. At present, the field star_age works as intended though.

[matthewturk]

@cphyc do you need a hand?

[cphyc]

@yt-fido test this please

[cphyc]

OK, now that #4880 is fixed, I think everything is working as intended now (finally)!

import yt
import numpy as np

@yt.particle_filter(requires=["conformal_birth_time"])
def star(pfilter, data):
    return data[pfilter.filtered_type, "conformal_birth_time"] != 0

@yt.particle_filter(requires=["conformal_birth_time"])
def DM(pfilter, data):
    return data[pfilter.filtered_type, "conformal_birth_time"] == 0

ds = yt.load(
    "output_00080",
    extra_particle_fields=[
        ("particle_birth_time", "d"),
        ("particle_metallicity", "d"),
    ],
)
ds.add_particle_filter("star")
ds.add_particle_filter("DM")
ad = ds.all_data()

print(f"{ds.current_time.to("Gyr")=}")
print(ad["io", "age"].min().to("Gyr"), ad["io", "age"].max().to("Gyr"))
print(ad["io", "particle_birth_time"].min().to("Gyr"), ad["io", "particle_birth_time"].max().to("Gyr"))
print(ad["io", "star_age"].min().to("Gyr"), ad["io", "star_age"].max().to("Gyr"))

This returns:

ds.current_time.to("Gyr")=unyt_quantity(11.92528501, 'Gyr')
0.0 Gyr 11.925285011256845 Gyr
0.0 Gyr 11.925285011256845 Gyr
1.7763568394002505e-15 Gyr 11.225119034912565 Gyr

So all fields now seem to be properly defined (the 1e-15Gyr in the star_age is due to interpolation).

---

Regarding star_age, it looks quite reasonable if e.g. I make a plot of the star formation history of the simulation.

import matplotlib.pyplot as plt
mask = ad["io", "conformal_birth_time"] != 0
tbins = np.arange(0, ds.current_time.to("Gyr").value, 0.1)
Mstar, _ = np.histogram(
    ad["star", "star_age"].to("Gyr").value,
    weights=ad["star", "particle_mass"].to("Msun").value,
    bins=tbins,
)
plt.step(
	tbins[1:],
	Mstar / np.diff(tbins) / 1e9,
)
plt.yscale("log")
plt.xlabel(r"Lookback time [$\mathrm{Gyr}$]")
plt.ylabel(r"SFR [$M_\odot/\mathrm{yr}$]")
plt.savefig("/tmp/SFR.png")