I haven’t looked at your code in detail yet, but this is complicated enough that it would help me to understand what’s going on if you can make your script locally runnable by me so I can poke at it with a debugger.
That means making your script use one of the data files we have for this purpose at
yt-project.org/data or by uploading one of your outputs to a publicly visible location:
$ yt upload some_file.tar.gz
You’ll need to share the link the above command prints out.
If you’re uncomfortable sharing your data publicly we can try off-list.
Nathan
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
_______________________________________________
yt-users mailing list -- yt-users@python.org
To unsubscribe send an email to yt-users-leave@python.org