Hi,
I am using YT to process the outputs of a non-standard RAMSES-RT simulation (non-standard in the sense that it uses custom particle/hydro variables), and it works excellently for almost everything. However I am now trying to make plots of the surface star formation rate density, and have hit an obstacle which I can't seem to surmount.
Specifically I want to estimate the SFRD by looking at all star particles that formed recently, which requires correctly processing the particle birth times. This seems to be a tricky business in RAMSES, as depending on the type of the simulation run (cosmological/RT/etc) the particle time units are different. Luckily for this particular case I have separate output files which contain the star particle formation times in the correct physical units (e.g. Gyr). What I would like to do is load these in, override the "star_age" field of the star particles, and then use that to filter out the star particles that formed recently. I have tried this and it seems to work well, up until I want to deposit this filtered particle field. At this point the code fails with the NeedsGridType and YTIllDefinedFilter exceptions.
Here's a minimal working example:
################################################################################
#!/usr/bin/env python3
import yt
from yt.data_objects.particle_filters import add_particle_filter
from yt import YTArray
import numpy as np
from scipy.io import FortranFile
# Load outputs of non-standard RAMSES-RT run
particles = [("particle_birth_time", "d"), ("particle_metallicity", "d"),
("particle_initial_mass", "d")]
fields = ["Density", "x-velocity", "y-velocity", "z-velocity", "nontp",
"Pressure", "Metallicity", "HI", "HII", "HeII", "HeIII",
"ivar_ref"]
ds = yt.load("output_00050/info_00050.txt",
fields=fields, extra_particle_fields=particles)
# Add initial filter to distinguish star particles
def _stars(pfilter, data):
return data[("io", "conformal_birth_time")] != 0 # star particles
add_particle_filter("stars", function=_stars, filtered_type="all",
requires=["conformal_birth_time"])
ds.add_particle_filter("stars")
# Get star formation times in Gyr
# (Based on pynbody's get_tform, read in star formation time from converted
# output)
def convert_times():
ages = np.array([])
base = "output_00050/birth_00050."
ncpu = 1024
nstar = 0
for i in range(ncpu):
fname = base + "out{:05d}".format(i + 1)
with FortranFile(fname) as birth_file:
read_ages = birth_file.read_reals(np.float64)
new = np.where(read_ages > 0)[0] # star particles only
nstar += len(new)
ages = np.append(ages, read_ages[new])
return nstar, ages
nstar, stellar_ages = convert_times()
# Override existing ("stars", "star_age") field with correct values
def _stellar_age(field, data):
return data.ds.current_time - YTArray(stellar_ages, "Gyr")
ds.add_field(("stars", "star_age"), function=_stellar_age,
units="Myr", particle_type=True, force_override=True)
# Filter stars born within past 10 Myr
def _new_stars(pfilter, data):
age = data[(pfilter.filtered_type, "star_age")].in_units("Myr")
mask = np.logical_and(age.in_units("Myr") <= 10, age > 0.0)
return mask
add_particle_filter("new_stars", function=_new_stars, filtered_type="stars",
requires=["star_age"])
ds.add_particle_filter("new_stars")
# Create new deposit field with estimate of recent star formation rate density
def _star_formation(field, data):
stellar_density = data[("deposit", "new_stars_density")]
stellar_age = ds.quan(10, "Myr")
return stellar_density / stellar_age
ds.add_field(("deposit", "sfr"), function=_star_formation,
units="Msun/kpc**3/yr", sampling_type="cell")
# Define a region of interest to project over
centre = ds.arr([0.5, 0.5, 0.5], "code_length")
fwidth = ds.quan(0.1, "code_length")
left_edge = centre - (fwidth/2)*1.1
right_edge = centre + (fwidth/2)*1.1
box = ds.box(left_edge=left_edge, right_edge=right_edge)
proj = ds.proj(field=("deposit", "sfr"), axis="z", weight_field=None,
center=centre, data_source=box)
################################################################################
The short summary of the error output is as follow:
################################################################################
yt.fields.field_exceptions.NeedsGridType: (0, None)
During handling of the above exception, another exception occurred:
yt.utilities.exceptions.YTIllDefinedFilter: Filter '<yt.data_objects.particle_filters.ParticleFilter object at 0x7f04591b2390>' ill-defined. Applied to shape (0, 3) but is shape (14150626,).
################################################################################
I was wondering if I am making a mistake somewhere, or if perhaps I'm simply going about this in the wrong way? The issue seems to revolve around using data[("deposit", "new_stars_density")]; if for example I instead use data[("deposit", "stars_density")] in the _star_formation() function then it all works fine.
Many thanks in advance for any help!
Lewis