Derived Fields and Profiles

One of the most powerful features in yt is the ability to create derived fields that act and look exactly like fields that exist on disk. This means that they will be generated on demand and can be used anywhere a field that exists on disk would be used. Additionally, you can create them by just writing python functions.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

import yt
from yt import derived_field

Derived Fields

This is an example of the simplest possible way to create a derived field. All derived fields are defined by a function and some metadata; that metadata can include units, LaTeX-friendly names, conversion factors, and so on. Fields can be defined in the way in the next cell. What this does is create a function which accepts two arguments and then provide the units for that field. In this case, our field is dinosaurs and our units are K*cm/s. The function itself can access any fields that are in the simulation, and it does so by requesting data from the object called data.

@derived_field(name="dinosaurs", units="K * cm/s", sampling_type="cell")
def _dinos(field, data):
    return data["gas", "temperature"] * data["gas", "velocity_magnitude"]

One important thing to note is that derived fields must be defined before any datasets are loaded. Let’s load up our data and take a look at some quantities.

ds = yt.load_sample("IsolatedGalaxy")
dd = ds.all_data()
print(list(dd.quantities.keys()))

['WeightedAverageQuantity', 'TotalQuantity', 'TotalMass', 'CenterOfMass', 'BulkVelocity', 'WeightedStandardDeviation', 'AngularMomentumVector', 'Extrema', 'SampleAtMaxFieldValues', 'MaxLocation', 'SampleAtMinFieldValues', 'MinLocation', 'SpinParameter']

One interesting question is, what are the minimum and maximum values of dinosaur production rates in our isolated galaxy? We can do that by examining the extrema quantity – the exact same way that we would for density, temperature, and so on.

print(dd.quantities.extrema(("gas", "dinosaurs")))

[3.31419530e+04 1.18723107e+15] K*cm/s

We can do the same for the average quantities as well.

print(
    dd.quantities.weighted_average_quantity(
        ("gas", "dinosaurs"), weight=("gas", "temperature")
    )
)

16822848849458.07 K*cm/s

A Few Other Quantities

We can ask other quantities of our data, as well. For instance, this sequence of operations will find the most dense point, center a sphere on it, calculate the bulk velocity of that sphere, calculate the baryonic angular momentum vector, and then the density extrema. All of this is done in a memory conservative way: if you have an absolutely enormous dataset, yt will split that dataset into pieces, apply intermediate reductions and then a final reduction to calculate your quantity.

sp = ds.sphere("max", (10.0, "kpc"))
bv = sp.quantities.bulk_velocity()
L = sp.quantities.angular_momentum_vector()
rho_min, rho_max = sp.quantities.extrema(("gas", "density"))
print(bv)
print(L)
print(rho_min, rho_max)

[-892739.08777744 1107086.20621902  620738.20523763] cm/s
[-5.01415153e+25  7.51695813e+26  4.16988689e+28] cm**2/s
2.2147823812377933e-28 g/cm**3 7.734265039240001e-24 g/cm**3

Profiles

yt provides the ability to bin in 1, 2 and 3 dimensions. This means discretizing in one or more dimensions of phase space (density, temperature, etc) and then calculating either the total value of a field in each bin or the average value of a field in each bin.

We do this using the objects Profile1D, Profile2D, and Profile3D. The first two are the most common since they are the easiest to visualize.

This first set of commands manually creates a profile object the sphere we created earlier, binned in 32 bins according to density between rho_min and rho_max, and then takes the density-weighted average of the fields temperature and (previously-defined) dinosaurs. We then plot it in a loglog plot.

prof = yt.Profile1D(
    sp, ("gas", "density"), 32, rho_min, rho_max, True, weight_field=("gas", "mass")
)
prof.add_fields([("gas", "temperature"), ("gas", "dinosaurs")])

fig, ax = plt.subplots()
ax.loglog(np.array(prof.x), np.array(prof["gas", "temperature"]), "-x")
ax.set(
    xlabel="Density $(g/cm^3)$",
    ylabel="Temperature $(K)$",
)

[Text(0.5, 0, 'Density $(g/cm^3)$'), Text(0, 0.5, 'Temperature $(K)$')]

Now we plot the dinosaurs field.

fig, ax = plt.subplots()
ax.loglog(np.array(prof.x), np.array(prof["gas", "dinosaurs"]), "-x")
ax.set(
    xlabel="Density $(g/cm^3)$",
    ylabel="Dinosaurs $(K cm / s)$",
)

[Text(0.5, 0, 'Density $(g/cm^3)$'), Text(0, 0.5, 'Dinosaurs $(K cm / s)$')]

If we want to see the total mass in every bin, we profile the mass field with no weight. Specifying weight=None will simply take the total value in every bin and add that up.

prof = yt.Profile1D(
    sp, ("gas", "density"), 32, rho_min, rho_max, True, weight_field=None
)
prof.add_fields([("gas", "mass")])

fig, ax = plt.subplots()
ax.loglog(np.array(prof.x), np.array(prof["gas", "mass"].in_units("Msun")), "-x")
ax.set(
    xlabel="Density $(g/cm^3)$",
    ylabel=r"Cell mass $(M_\odot)$",
)

[Text(0.5, 0, 'Density $(g/cm^3)$'), Text(0, 0.5, 'Cell mass $(M_\\odot)$')]

In addition to the low-level ProfileND interface, it’s also quite straightforward to quickly create plots of profiles using the ProfilePlot class. Let’s redo the last plot using ProfilePlot

prof = yt.ProfilePlot(sp, ("gas", "density"), ("gas", "mass"), weight_field=None)
prof.set_unit(("gas", "mass"), "Msun")
prof.show()

Field Parameters

Field parameters are a method of passing information to derived fields. For instance, you might pass in information about a vector you want to use as a basis for a coordinate transformation. yt often uses things like bulk_velocity to identify velocities that should be subtracted off. Here we show how that works:

sp_small = ds.sphere("max", (50.0, "kpc"))
bv = sp_small.quantities.bulk_velocity()

sp = ds.sphere("max", (0.1, "Mpc"))
rv1 = sp.quantities.extrema(("gas", "radial_velocity"))

sp.clear_data()
sp.set_field_parameter("bulk_velocity", bv)
rv2 = sp.quantities.extrema(("gas", "radial_velocity"))

print(bv)
print(rv1)
print(rv2)

[ -8934.11243694 198277.64189059 311013.85834147] cm/s
[-69582773.1781608   32330373.38084824] cm/s
[-69649428.21589676  32236297.22862757] cm/s