Reading Data using YT

Install yt

https://yt-project.org/#getyt

I recommend installing yt using conda. I also installed Python and managed my Python packages and versions using conda. This will guarantee that all the packages (such as yt, numpy, scipy etc) work together.

Loading Athena++ output

In Python, first import yt:

import yt

Then simply use yt to load your file (hdf5 definitely work, I think vtk works too):

filename="output.athdf"

ds = yt.load(filename)

This will load the generic data without units specified. I usually set the units already when I load the data. For example, using:

unit_base={"length_unit": (1.0,"pc"), "time_unit": (1.0,"s*pc/km"), "mass_unit": (2.38858753789e-24,"g/cm**3*pc**3")}

ds = yt.load(filename, units_override=unit_base)

Reading data into NumPy arrays

ds stores the metadata of the output data set. For example, you can use

print(ds.field_list)

to see all the fields that the data set contains.

To read a field, for example, the density field (in Athena++ the name is "rho"), you need to take a slice of the data first:

grid = ds.covering_grid(level=0, left_edge=left_edge, dims=dims)

If you want to read the whole grid, you can set:

left_edge = ds.domain_left_edge

dims = ds.domain_dimensions

Now, finally, you have the density field,

rho = grid["rho"] # rho is a YTArray with units attached

However, the variable "rho" here has a unit attached, usually in code unit as a default. If you want to manipulate a NumPy array without units, you should do

rho = grid["rho"].in_units("g/cm**3").to_ndarray() # rho is a NumPy array

Add derived fields

Sometimes, you want to add some field derived from the default output fields. For example, you often work with both the physical density in units of g/cm**3, as well as the density that is expressed in the unit of the number density of hydrogen atoms. You can add a derived field called "nH" when you read the data. First, you need to define a function:

from yt.utilities.physical_constants import mh

def _ndensity(field, data): return data["rho"]/(1.4271*mh)

Then you can add the field:

ds.add_field(("nH"), function=_ndensity, units='cm**(-3)',display_name=r'$n_{\rm H}$', sampling_type="cell")