Derived fields in yt for RAMSES-RTZ

Overview

Make derived fields in yt for simulations run with the RAMSES-RTZ code. Currently useful for simulations such as the MEGATRON suite, but can in principle be used for any RAMSES-RTZ simulation with the same thermochemistry and stellar SEDs.

Spectral derived fields

• Emission line module for collisional and recombination lines in the optical and UV (requires atomic data generated by the chianti package, see How to use below)
• Nebular and stellar continuum spectra (requires data on glamdring/infinity, see How to use below))
• Corrections for unresolved stromgren spheres (only emission lines for now, will be extended to the nebular continuum in the future)
• INCOMING: Ability to generate IFU spectra using bins of cells. This requires finished functionality of yt-project/yt#5227.

Star derived fields

• Correct star particle ages
• Filters for Pop. II and Pop. III stars
• Stellar spectra

Chemistry derived fields

• Electron number density
• Mean molecular weight
• Metallicity

Cooling derived fields

• H2 cooling reconstruction
• Full metal cooling

Radiative transfer derived fields

• Energy density of the radiation field for each frequency bin

Loading cutouts

Cutouts generated for the Megatron simulations can be read using yt_derived_fields.cutouts.loader.load_cutout as follows:

from yt_derived_fields.cutouts.loader import load_cutout
import unyt as u
import yt

ds = load_cutout("MEGATRON_CP_NEW/halo_cutouts/output_00094/halo_3378_gas.bin", boxsize=50 * u.Mpc)

# Due to the way data are loaded, we'll have *many* nans (in regions where no data were saved)
# Let's filter them out
ad = ds.all_data().exclude_nan(("gas", "density"))

p = yt.ProjectionPlot(ds, "x", ("gas", "density"), data_source=ad, center=ds.domain_center, width=(20, "kpc"))
p.save()

How to use

1. Install yt: Since the changes to YT can be very specific to the MEGATRON simulation, it is recommended to use a branch of yt that has been modified for this purpose. The branch megatron_rtFix on Anatole Storck's yt fork in GitHub contains the necessary changes. Here is how to clone and install it for your python environment.

   git clone https://github.com/AnatoleStorck/yt.git
   cd yt
   git switch megatron_rtFix
   pip install -e .

2. Clone and install this repository and copy the TOML configuration specific to MEGATRON for global yt configuration. If you don't want to make the configuration global, you can also copy the yt.toml file to your analysis directory.

   git clone https://github.com/AnatoleStorck/yt_derived_fields.git
   cd yt_derived_fields
   pip install -e .

   # This will make the configuration global. If ~/.config/yt/ doesn't exist, you may need to create it first.
   cp yt_derived_fields/megatron_derived_fields/yt.toml ~/.config/yt/.

3. Generate atomic data for observational derived fields: The emission line module requires atomic data generated by the chianti package. On glamdring/infinity, chianti databases are already set up. If you want to generate your own chianti database, follow the instructions on the Chianti database website. You can then run the spectral_utils/generate_atomic_data.py script to generate the atomic data for the emission line module (if you have made a custom chianti database, make sure to set the XUVTOP environment variable to point to it before running the script).

   # This will create the collisional and recombination line data for the emission line module.
   # Only Optical and UV lines are currently included, but more (such as infrared lines) can be added in the future.
   python yt_derived_fields/spectral_utils/generate_atomic_grids.py

   The script currently mostly generates optical/UV lines, but more lines can be added in the future (specifically more infrared lines).

4. Use the derived field functions: Generate new fields to be used with various yt functions. Chemistry fields include the electron number density, mean molecular weight, molecular hydrogen number density, and others. Cooling includes two H2 cooling fields, full metal cooling, and others. Star derived fields include a proper way to derive the star particle ages and to get spectra.

   import yt
   import pandas as pd
   
   import yt_derived_fields.generate_fields as gen_fields
   
   ds = yt.load(<path_to_your_favorite_megatron_output>)
   hc = pd.read_csv(<path_to_your_halo_catalog_in_pandas_format>)
   
   # Generate all chemistry, cooling, radiative, stellar fields. Along with all emission lines and nebular continuum fields
   # This will create a lot of fields, but nothing is actually calculated until you access the field for the first time.
   gen_fields.create_derived_fields(
      ds,
      simple_ne=False,              # Whether to derive ne using only H and He (no metals).
      H2_cooling="moseley",         # Make sure you have the right cooling model for your simulation.
      pop3_stars=False,             # All CCM runs have pop3 stars, all CP runs do not.
      stromgren_correction=True,    # Replace the emissivity of cells with unresolved stromgren spheres.
      parallel=True,                # Parallel processing for generating the nebular continuum and stellar spectra.
   )
   
   # Create a sphere around a halo from the halo catalog
   # (This could also be any other sized region you want, see the yt documentation for more details on how to create data objects)
   some_halo = hc[hc["id"] == 42]
   halo_sphere = ds.sphere(([some_halo.x,
                             some_halo.y,
                             some_halo.z], "Mpccm/h"),
                            (some_halo.r200b, "kpccm/h"))
   
   # An example field from each module
   ne = halo_sphere["gas", "electron_number_density"]                   # the electron number density in each cell
   cool_H2 = halo_sphere["gas", "cooling_H2"]                           # The H2 cooling rate in each cell
   O3_5007 = halo_sphere["O3-5007_luminosity"]                          # The O3-5007 line luminosity in every cell
   LW_photon_density = halo_sphere["rt", "photon_energy_density_LW"]    # The energy density of the Lyman-Werner radiation field in each cell
   pop2_spectra = halo_sphere["pop2", "spectra"]                        # The spectra of every pop2 star particle in the sphere
   nebular_continuum = halo_sphere["gas", "nebc_total"]                 # The nebular continuum emission in every cell