You can interact with this notebook online: Launch notebook

Converting Arepo snapshots to be usable by TARDIS

Note

This notebook cannot be run interactively!

What is Arepo?

Arepo is a massively parallel gravity and magnetohydrodynamics code for astrophysics, designed for problems of large dynamic range. It employs a finite-volume approach to discretize the equations of hydrodynamics on a moving Voronoi mesh, and a tree-particle-mesh method for gravitational interactions. Arepo is originally optimized for cosmological simulations of structure formation, but has also been used in many other applications in astrophysics.

[WSP20]

This parser is intended for loading Arepo output files (‘snapshots’), extracting the relevant (line-of-sight dependent) data and exporting it to csvy files, which can in turn be used in TARDIS models (see CSVY model).

Note

This parser has been developed for the (not publicly available) development version of Arepo, not the public version. Although it should also work with snapshots from the public version, this has not been tested. If you run into trouble loading the snapshot using the built-in functions, try providing the data manually.

import numpy as np
import matplotlib.pyplot as plt
from tardis.io.parsers import arepo
import json
%matplotlib inline

Loading the simulation data

As a first step, the relevant data has to be loaded from an Arepo snapshot, i.e. an output file of Arepo. In case you have the arepo-snap-util package installed, you can use the built in wrapper (as described below) to load the relevant data. In case you do not have this package installed or want to load the snapshot in a different way, you can manually provide the relevant data and continue with the next step.

If you’re using the built-in tool, you will need to provide the path to the snapshot file, a list with the elements you want to include in your TARDIS model, as well as the species file with which the Arepo simulation was run. (The species file should contain one header line, followed by two columns, where the first contains the names of all the species and is used to find the indices of the individual species within the snapshot. The second column is not used by the loader.)

In case you have the arepo-snap-util package installed, you can load the data directly from a snapshot:

snapshot = arepo.ArepoSnapshot(
    "arepo_snapshot.hdf5", ["ni56", "si28"], "species55.txt"
)
pos, vel, rho, mass, xnuc, time = snapshot.get_grids()

This will load the necessary from the snapshot. See the API documentation for more options on how to load snapshots.

Note

The built in functionality only loads Arepo particle species 0, i.e. gas particles. For most (if not all) purposes connected to TARDIS, gas particles should be the only species of interest. In case you need to load a different particle species, please manually load the data.

This will fail with an error if you do not have the arepo-snap-util package installed. Since this is not a default dependency of TARDIS, the data can also be loaded manually. (This manual load can effectively be used to load all kinds of models unrelated to Arepo, as long as the data comes in the correct format.)

The data could be loaded from a previously generated json file, where all necessary input from the snapshot has been dumped in. (The ``json`` file is only an example and more for illustrative purposes. As long as the data is provided in the correct format (see below) it is up to the user how it is saved/loaded.)

Code: (click to expand)

### Dumping the data:
data = {
    "pos" : pos.tolist(),
    "vel" : vel.tolist(),
    "rho" : rho.tolist(),
    "mass" : mass.tolist(),
    "xnuc": [xnuc[x].tolist() for x in list(xnuc.keys())],
    "time": time,
}
json_string = json.dumps(data)
with open('arepo_snapshot.json', 'w') as outfile:
    json.dump(json_string, outfile)

### Loading the data:
with open('arepo_snapshot.json') as json_file:
    data = json.loads(json.load(json_file))

# The following lines only parse the .json file. You might not need this depending on how you saved
# the snapshot data
pos, vel, rho, mass, nucs, time = data["pos"], data["vel"], data["rho"], data["mass"], data["xnuc"], data["time"]
pos = np.array(pos)
vel = np.array(vel)
rho = np.array(rho)
mass = np.array(mass)

# The nuclear data should be in a dict where each element has its own entry (with the key being the element name)
xnuc = {
    "ni56" : np.array(nucs[0]),
    "si28" : np.array(nucs[1]),
}

print("Position data shape: ", pos.shape)
print("Velocity data shape: ", vel.shape)
print("Density data shape: ", rho.shape)
print("Mass data shape: ", mass.shape)
print("Nuclear data shape (per element): ", xnuc["ni56"].shape)

Position data shape:  (3, 1218054)
Velocity data shape:  (3, 1218054)
Density data shape:  (1218054,)
Mass data shape:  (1218054,)
Nuclear data shape (per element):  (1218054,)

In case you want to load the snapshot data itself with your own tools, you will need to provide the following data: - Position in the center of mass frame (“pos”) -> np.array - Velocity (“vel”) -> np.array - Density (“rho”) -> np.array - Mass (“mass”) -> np.array - Nuclear fraction (“xnuc”) of each element you want to include -> dict containing np.array

- Time of the snapshot (“time”) -> float

The data should have the same shape as the one provided by the built-in tool (see the given shapes above for an example).

Extracting a profile and converting it to a csvy file

Now You can create a TARDIS model. There are two possibilities on how to extract the profiles from the snapshot:

• Cone profile: This extracts the data within a specified cone

• Full profile: This averages over the whole simulation

profile = arepo.ConeProfile(pos, vel, rho, mass, xnuc, time)

This loads the data (in this example for a cone profile), which can then be cut to the ranges which you want to include in your TARDIS model. The syntax for the other profiles is similar:

• arepo.FullProfile(<args>)

Next you can create the profiles according to the model option you selected. A diagnostic plot will be shown per default, but this behaviour can be turned off with the option show_plot=False. The plot will always show both the positive and negative axis.

Note

The keyword opening_angle=40 is only needed for the cone profile. The other modes do not accept this keyword! The angle itself is the opening angle of the full cone and NOT the angle between the central x-axis and the cone!

profile.create_profile(opening_angle=40)
profile.rebin(100)
profile.plot_profile()
plt.show()

Plotting the data

Using profile.plot_profile(save=<filename>, dpi=<dpi>) you can create a plot of the exported profile. Due to the large number of cells usually used in Arepo simulations, it is recommended to plot the data only after rebinning the profiles or exporting them to a csvy file (which automatically rebins the data).

Manually rebinning the data

Using profile.rebin(<nshells>), you can manually rebin the data , using the mean value of the data in each bin. Rebinning only works when decreasing the number of bins. In case you need to increase resolution, you will need to create a new profile.

Selecting a specific region

In many cases you only want a very specific region from the snapshot, e.g. cutting out the dense, optically thick regions. This can be achieved using the keywords inner_radius and outer_radius.

profile.create_profile(opening_angle=40, inner_radius=1e11, outer_radius=2e11)
profile.rebin(100)
profile.plot_profile()
plt.show()

Once you have created a profile of the desired region, you can export the profile to a csvy file, which in turn can be used in a TARDIS model. Here you have to specify how many shells you want to export. The profiles are rebinned automatically using Scipys binned_statistic function, using the mean value of the data in each bin.

profile.export(20, "snapshot_converted_to_tardis.csvy", overwrite=True)

'snapshot_converted_to_tardis.csvy'

During the export, the yaml header is automatically written and includes the time of the screenshot as the time for both the nuclear data as well as the density profile. Per default, the positive axis will be exported. the negative axis can be exported with direction="neg".

All abundences will normalised such that roughly sum to 1, but slight deviations are expected to occur.

Note

The export function will not overwrite existing files, but rather add an incrementing number to the filename.

Using the parser as a command line tool

You can also use the tardis.io.arepo package as a command line utility, e.g. if you are running batch jobs and want to automatically convert your snapshots from within you job-script.

To export the same profile as in the example above you can run:

python ./<location_of_arepo_parser>/arepo.py snapshot.hdf5 snapshot_converted.csvy -o 40 -n 20 --inner_radius 1e11 --outer_radius 2e11 -e ni56 si28 --plot plot.png

This will also save diagnostic plots of both the raw and rebinned profiles. For more information on how to use the command line tool run:

python ./<location_of_arepo_parser>/arepo.py --help

Note

The command line tool does only work with snapshot files, not with e.g. json files. It is in any case not advised to use json files as an intermediate step, as the become huge for higher resolutions.