User Guide

This example goes through the main usages of gizio. To follow along, make sure gizio is installed, and the sample data is downloaded.

import gizio

Snapshot

First load the snapshot

snap = gizio.load('data/FIRE_M12i_ref11/snapshot_600.hdf5')

Inspection

Once a snapshot is loaded, several of its general peoperties are extracted as attributes, and could be inspected like the following:

snap.header

{'time': unyt_quantity(13.79874688, 'Gyr'),
 'n_file': 1,
 'mass_tab': unyt_array([0., 0., 0., 0., 0., 0.], 'code_mass'),
 'f_sfr': 1,
 'f_cool': 1,
 'f_fb': 1,
 'f_age': 1,
 'f_met': 11,
 'n_part': array([ 753678, 1104128, 2567905,       0,  361239,       0], dtype=uint32),
 'n_part_pf': [array([ 753678, 1104128, 2567905,       0,  361239,       0], dtype=int32)],
 'box_size': unyt_quantity(60000., 'code_length'),
 'Om0': 0.272,
 'OmL': 0.728,
 'h': 0.702,
 'z': 0.0,
 'cosmological': True}

snap.shape

OrderedDict([('PartType0', 753678),
             ('PartType1', 1104128),
             ('PartType2', 2567905),
             ('PartType3', 0),
             ('PartType4', 361239),
             ('PartType5', 0)])

snap.cosmology

LambdaCDM(H0=70.2 km / (Mpc s), Om0=0.272, Ode0=0.728, Tcmb0=0 K, Neff=3.04, m_nu=None, Ob0=None)

Field Access

Fields could be accessed through a dictionary interface. The available keys could be queried by

snap.keys()

[('PartType0', 'Coordinates'),
 ('PartType0', 'Density'),
 ('PartType0', 'ElectronAbundance'),
 ('PartType0', 'InternalEnergy'),
 ('PartType0', 'Masses'),
 ('PartType0', 'Metallicity'),
 ('PartType0', 'NeutralHydrogenAbundance'),
 ('PartType0', 'ParticleIDs'),
 ('PartType0', 'Potential'),
 ('PartType0', 'SmoothingLength'),
 ('PartType0', 'StarFormationRate'),
 ('PartType0', 'Velocities'),
 ('PartType1', 'Coordinates'),
 ('PartType1', 'Masses'),
 ('PartType1', 'ParticleIDs'),
 ('PartType1', 'Potential'),
 ('PartType1', 'Velocities'),
 ('PartType2', 'Coordinates'),
 ('PartType2', 'Masses'),
 ('PartType2', 'ParticleIDs'),
 ('PartType2', 'Potential'),
 ('PartType2', 'Velocities'),
 ('PartType4', 'Coordinates'),
 ('PartType4', 'Masses'),
 ('PartType4', 'Metallicity'),
 ('PartType4', 'ParticleIDs'),
 ('PartType4', 'Potential'),
 ('PartType4', 'StellarFormationTime'),
 ('PartType4', 'Velocities')]

A field could be accessed through

snap['PartType0', 'Coordinates']

unyt_array([[30711.56320229, 33303.32765069, 33743.24254345],
            [30801.09893721, 33396.72993339, 33673.13785209],
            [30721.03179246, 33313.72662534, 33657.53180413],
            ...,
            [28702.32012144, 32308.36857189, 32512.44014079],
            [28695.79401442, 32304.20735574, 32514.23300421],
            [28700.69483794, 32308.30239816, 32518.00387993]], 'code_length')

Note that the field is loaded as a unyt array.

Particle Selector

Particle selector is a more flexible interface to access fields. Default particle selectors according to particle types are created automatically when loading the snapshot, and could be accessed from

snap.pt

{'gas': <gizio.core.ParticleSelector at 0x7f96d7855f28>,
 'hdm': <gizio.core.ParticleSelector at 0x7f9700eddef0>,
 'ldm': <gizio.core.ParticleSelector at 0x7f9700ed5eb8>,
 'star': <gizio.core.ParticleSelector at 0x7f9700ed5d68>,
 'all': <gizio.core.ParticleSelector at 0x7f9700ed5da0>}

Let’s have a look at the gas particle selector:

gas = snap.pt['gas']
gas.keys()

dict_keys(['z', 'Potential', 'p', 'ne', 'v', 'sfr', 'u', 'nh', 'h', 'id', 'rho', 'm', 't'])

Different from the snapshot interface, shorthand keys are used in the particle selector interface, except for unrecognized ones like 'Potential'. So we have the following equality:

(gas['p'] == snap['PartType0', 'Coordinates']).all()

True

The number of selected particles can be gotten from

len(gas)

753678

Some other particle selector properties could be queried like the following:

gas.pmask

OrderedDict([('PartType0',
              array([ True,  True,  True, ...,  True,  True,  True])),
             ('PartType1', None),
             ('PartType2', None),
             ('PartType3', None),
             ('PartType4', None),
             ('PartType5', None)])

gas.shape

OrderedDict([('PartType0', 753678),
             ('PartType1', 0),
             ('PartType2', 0),
             ('PartType3', 0),
             ('PartType4', 0),
             ('PartType5', 0)])

gas.direct_fields()

{'z': 'Metallicity',
 'Potential': 'Potential',
 'p': 'Coordinates',
 'ne': 'ElectronAbundance',
 'v': 'Velocities',
 'sfr': 'StarFormationRate',
 'u': 'InternalEnergy',
 'nh': 'NeutralHydrogenAbundance',
 'h': 'SmoothingLength',
 'id': 'ParticleIDs',
 'rho': 'Density',
 'm': 'Masses'}

Derived Fields

Among the listed keys, one is special: 't' for temperature. There isn’t a ('PartType0', 'Temperature') field in the snapshot. Actually, temperature is computed from internal energy. So we call it a derived field:

gas['t']

unyt_array([ 59962.26953125,  84842.5 ,  82207.421875, ...,
            121806.1171875, 157144.0625, 134721.953125], dtype=float32, units='K')

We could check the code that does the computation (be careful this is not part of the public interface yet):

import inspect
print(inspect.getsource(gas._field_registry['t']))

    @staticmethod
    def compute_temperature(ps):
        """Compute temperature from internal energy.

        Parameters
        ----------
        ps : ParticleSelector
            A particle selector.

        Returns
        -------
        unyt.unyt_array
            The temperature array.

        """
        # See the note following InternalEnergy on this page:
        # http://www.tapir.caltech.edu/~phopkins/Site/GIZMO_files/gizmo_documentation.html#snaps-reading
        ne = ps["ne"]
        u = ps["u"]
        z_he = ps["z"][:, 1]

        from unyt.physical_constants import kb, mp

        gamma = 5 / 3

        y = z_he / (4 * (1 - z_he))
        mu = (1 + 4 * y) / (1 + y + ne)
        temperature = mu * mp * (gamma - 1) * u / kb
        return temperature.to("K")

It’s also pretty easy to register a custom derived field:

def volume(field_system):
    return (field_system['m'] / field_system['rho']).to('kpc**3')

gas.register_field('vol', volume)
gas['vol']

unyt_array([2.0862592e+06, 2.6118410e+06, 1.1494234e+06, ...,
            1.3297848e+03, 1.0290620e+03, 1.3337478e+03], dtype=float32, units='kpc**3')

Boolean Masking

The particle selection could be further refined by boolean masking. For example, to select hot gas, we could do:

hot_gas = gas[gas['t'].to_value('K') > 1e5]
hot_gas['t'].min()

unyt_quantity(100000.79, dtype=float32, units='K')

Set-like Composing

Another way to construct a new particle selector is to compose existing ones. A set-like interface is implemented for particle selectors to compose. For example, to define baryon as gas or star:

star = snap.pt['star']
baryon = gas | star

Only keys from both gas and star are kept in the result:

set(baryon.keys()) == set(gas.keys()) & set(star.keys())

True

Here’s another example to define dark matter as all but baryon:

dm = snap.pt['all'] - baryon

Snapshot Format Spec

Snapshot format are specified by a class. A GIZMO spec is built into gizio. It could be accessed as an attribute:

snap.spec

<gizio.spec.GIZMOSpec at 0x7f96d7855e48>

To customize a spec, subclass SpecBase or GIZMOSpec and modify its class attributes. For example, here we subclass GIZMOSpec:

from gizio.spec import GIZMOSpec

class MySpec(GIZMOSpec):
    pass

MySpec.FIELD_SPEC

[('Coordinates', 'p', 'code_length'),
 ('Velocities', 'v', 'code_velocity'),
 ('ParticleIDs', 'id', ''),
 ('Masses', 'm', 'code_mass'),
 ('InternalEnergy', 'u', 'code_specific_energy'),
 ('Density', 'rho', 'code_mass / code_length**3'),
 ('SmoothingLength', 'h', 'code_length'),
 ('ElectronAbundance', 'ne', ''),
 ('NeutralHydrogenAbundance', 'nh', ''),
 ('StarFormationRate', 'sfr', 'Msun / yr'),
 ('Metallicity', 'z', 'code_metallicity'),
 ('ArtificialViscosity', 'alpha', ''),
 ('MagneticField', 'b', 'code_magnetic_field'),
 ('DivergenceOfMagneticField', 'divb', 'code_magnetic_field / code_length'),
 ('StellarFormationTime', 'sft', ''),
 ('BH_Mass', 'mbh', 'code_mass'),
 ('BH_Mdot', 'mdot', 'code_mass / code_time'),
 ('BH_Mass_AlphaDisk', 'mad', 'code_mass')]

And then modify the field spec so that 'ParticleIDs' is aliased to 'pid' instead of 'id':

MySpec.FIELD_SPEC[2] = ('ParticleIDs', 'pid', '')
MySpec.FIELD_SPEC

[('Coordinates', 'p', 'code_length'),
 ('Velocities', 'v', 'code_velocity'),
 ('ParticleIDs', 'pid', ''),
 ('Masses', 'm', 'code_mass'),
 ('InternalEnergy', 'u', 'code_specific_energy'),
 ('Density', 'rho', 'code_mass / code_length**3'),
 ('SmoothingLength', 'h', 'code_length'),
 ('ElectronAbundance', 'ne', ''),
 ('NeutralHydrogenAbundance', 'nh', ''),
 ('StarFormationRate', 'sfr', 'Msun / yr'),
 ('Metallicity', 'z', 'code_metallicity'),
 ('ArtificialViscosity', 'alpha', ''),
 ('MagneticField', 'b', 'code_magnetic_field'),
 ('DivergenceOfMagneticField', 'divb', 'code_magnetic_field / code_length'),
 ('StellarFormationTime', 'sft', ''),
 ('BH_Mass', 'mbh', 'code_mass'),
 ('BH_Mdot', 'mdot', 'code_mass / code_time'),
 ('BH_Mass_AlphaDisk', 'mad', 'code_mass')]

Load the snapshot with this new spec and we see 'pid' in the keys.

snap = gizio.load('data/FIRE_M12i_ref11/snapshot_600.hdf5', spec=MySpec())
snap.pt['gas'].keys()

dict_keys(['z', 'Potential', 'p', 'ne', 'v', 'sfr', 'u', 'nh', 'h', 'pid', 'rho', 'm', 't'])