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