Processing a simulation snapshot#

You’ve found the tutorials for learning more about using hydrodynamical zoom-in snapshot in cogsworth simulation! In this first part of two we’ll talk about how to prepare a snapshot and extract the necessary information.

Learning Goals#

By the end of this tutorial you should know:

How to prepare a snapshot for use in cogsworth
How to compute the galactic potential of a snapshot
How to access the initial state of a subset of snapshot particles

Beware - extra dependencies required here!

You’ll need to have installed the extra dependencies of cogsworth to postprocess hydrodynamical simulations. Check out the installation page for more details on how to do this! (You’ll probably just need to run pip install cogsworth[extras])

Note

Note that running this notebook requires a hydrodynamical zoom-in snapshot. The original was run using FIRE-2 m11h, snapshot 600, which can be downloaded from FlatHub.

[1]:

import cogsworth
import pynbody
import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
import gala.potential as gp

[2]:

# this all just makes plots look nice
%config InlineBackend.figure_format = 'retina'

plt.rc('font', family='serif')
plt.rcParams['text.usetex'] = False
fs = 24

# update various fontsizes to match
params = {'figure.figsize': (12, 8),
          'legend.fontsize': fs,
          'axes.labelsize': fs,
          'xtick.labelsize': 0.9 * fs,
          'ytick.labelsize': 0.9 * fs,
          'axes.linewidth': 1.1,
          'xtick.major.size': 7,
          'xtick.minor.size': 4,
          'ytick.major.size': 7,
          'ytick.minor.size': 4}
plt.rcParams.update(params)

The Goal#

Regular Population classes use empirically informed analytical distributions for the star formation history and galactic potential and sample either a fixed total number or mass of systems. This is great for many purposes, but in some cases it may be important to consider formation scenarios that properly account for initial spatial clustering and more a complex galactic potential. This is where the hydro module comes in use.

Using this module you can use hydrodynamical zoom-in simulations as the initial conditions/galactic potential for a cogsworth simulation. In this tutorial we’ll focus on preparing simulation snapshots for use in cogsworth and in the next we’ll talk about building a HydroPopulation.

Preparing a snapshot#

Each snapshot uses an arbitrary coordinate system for their simulation box. However, we would like to focus on a particlar halo in these simulations and convert everything to those galactocentric coordinates.

cogsworth makes this easy for you with cogsworth.hydro.utils.prepare_snapshot(). This provides a wrapper over a couple of pynbody functions - all you need to do is specifiy the snapshot file, or folder in which to find and any parameters you want to pass to the halo finder.

Let’s try this out on FIRE’s m11h!

Note

In order to process FIRE snapshots with pynbody (which is what cogsworth uses under the hood), you need to create a .pynbodyrc file in your home directory with the following content:

[gadget-units]
pos: kpc a h^-1

[gadgethdf-type-mapping]
dm: PartType1

[3]:

snap = cogsworth.hydro.utils.prepare_snapshot("../../../data/snapshot_600.hdf5")

cogsworth warning: Looks like you're loading a snapshot that doesn't specify its units, I'm going to infer them but make sure the outputted units looks right!
cogsworth warning: I couldn't find a halo catalogue, so I'll use pynbody's built-in halo finder to centre the snapshot

This function has just:

centred the snapshot on the primary halo
rotated the coordinates to put the halo edge-on and then face-on
converted everything to physical units

Note the warnings that show up as well - FIRE-2 simulations don’t store units in the output files and thus pynbody has to infer the units based on some assumed defaults. Additionally, I didn’t download an associated halo catalogue and instead use a shrinking sphere method to find the main halo center.

This function has returned a pynbody.snapshot.SimSnap instance and you can do all of the usual pynbody things with this.

Let’s check that our function did what we wanted by plotting an image of the gas!

[4]:

pynbody.plot.image(snap.gas, width=40, resolution=1000);

/home/tom/miniconda3/envs/cogsworth/lib/python3.10/site-packages/pynbody/snapshot/gadgethdf.py:394: UserWarning: Unable to infer units from HDF attributes
  warnings.warn("Unable to infer units from HDF attributes")

../../_images/tutorials_hydro_prep_9_1.png

Tip

If you want to learn more about using pynbody you should go and check out their awesome documentation

Computing the potential#

Now that the snapshot is ready to go we can work out a representation of the galactic potential. Using a self-consistent field method, one can determine the galactic potential resulting from all of the particles in the simulation. By default cogsworth constructs this potential using all of the particles in a snapshot as follows

[7]:

pot = gp.load("../../../cogsworth/hydro/m11h.yml")

[ ]:

pot = cogsworth.hydro.potential.get_snapshot_potential(snap)

[19]:

fig, ax = plt.subplots(figsize=(10, 10))

# define a basic grid to use
grid = np.linspace(-5, 5, 100) * u.kpc

# plot using the grid for x, y, but fixing z=0 kpc
pot.plot_contours(grid=(grid, grid, 0), ax=ax, levels=15);

ax.set(xlabel=r'$x \, [\rm kpc]$', ylabel=r'$y \, [\rm kpc]$')

plt.show()

../../_images/tutorials_hydro_prep_14_0.png

You can also change the number of components to use for the SCF function, for example we can decrease nmax and lmax to create a lower resolution potential.

[8]:

low_res_pot = cogsworth.hydro.potential.get_snapshot_potential(snap, nmax=2, lmax=1)

You can similarly find the potential of a subset of particles, say only every 1000th dark matter particle for instance

[9]:

weird_pot = cogsworth.hydro.potential.get_snapshot_potential(snap[::1000],
                                                             components=[{"attr": "dm",
                                                                          "label": "dark_matter",
                                                                          "r_s": 3}])

Rewinding particles#

When initialising a population from a snapshot, you’ll want to use the initial state of each particle rather than its present day state (which you’ll find the snapshot). Therefore you can use cogsworth.hydro.rewind.rewind_to_formation() to integrate the orbits of particles back to their starting location (using the potential we just calculated) and access their formation masses.

Tip

Doing this for every particle would be very expensive, so in most cases you’ll likely want to operate on a subset of your (star) particles

Let’s try rewinding just the first 3 particles of the simulation

[17]:

snap.s[:3]["iord"]

[17]:

SimArray([10846552,  5726876,  6569263], dtype=uint32)

[20]:

init_star_particles = cogsworth.hydro.rewind.rewind_to_formation(snap.s[:3], pot=pot,
                                                                 dt=-1 * u.Myr, processes=1)

100%|█████████████████████████████████████████████| 3/3 [00:09<00:00,  3.19s/it]

[21]:

init_star_particles

[21]:

	mass	Z	t_form	x	y	z	v_x	v_y	v_z
id
10846552	5054.304478	0.006751	3.874034	-39.968784	37.217411	-100.197458	26.202575	-12.450640	36.281095
5726876	8013.486763	0.014378	13.735638	-6.546165	-4.008135	0.889492	46.652619	-79.332788	39.956209
6569263	9579.370813	0.017819	13.735705	-6.535423	-4.008934	0.891565	43.619342	-80.168116	37.953939

These star particles can then be used to sample a population (we’ll do that in another tutorial).

In reality you’ll likely want to focus on a less contrived subset of your star particles - perhaps focusing on only the stars formed recently like this?

[22]:

recent_stars = snap.s[(snap.properties["time"] - snap.s["tform"]).in_units("Myr") < 50]
recent_stars

[22]:

<SimSnap "../../../data/snapshot_600.hdf5::star:indexed" len=3316>

Wrap-up#

And that’s all for this first tutorial of the hydrodynamical zoom-in series, nice job! Press your right arrow key to head over to part two where we use what we’ve done here to create a full blown cogsworth population!

Note

This tutorial was generated from a Jupyter notebook that can be found here.