In yt 3.0, we want to make it easier to access "raw" simulation data that a code writes directly to disk. The new unit system makes it much simpler to convert back and forth between physical coordinates and the unscaled "raw" coordinate system used internally in the simulation code. In some cases, this conversion involves transforming to comoving coordinates, so that is also covered here.

Code units¶

Let's take a look at a cosmological enzo dataset to play with converting between physical units and code units:

import yt
ds = yt.load('Enzo_64/DD0043/data0043')

/usr/lib64/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters

The conversion factors between Enzo's internal unit system and the physical CGS system are stored in the dataset's unit_registry object. Code units have names like code_length and code_time. Let's take a look at the names of all of the code units, along with their CGS conversion factors for this cosmological enzo dataset:

reg = ds.unit_registry

for un in reg.keys():
    if un.startswith('code_'):
        fmt_tup = (un, reg.lut[un][0], str(reg.lut[un][1]))
        print ("Unit name:      {:<15}\nCGS conversion: {:<15}\nDimensions:     {:<15}\n".format(*fmt_tup))

Unit name:      code_length    
CGS conversion: 5.5551728502644205e+26
Dimensions:     (length)       

Unit name:      code_mass      
CGS conversion: 4.890451593516587e+50
Dimensions:     (mass)         

Unit name:      code_density   
CGS conversion: 2.8527008703353168e-30
Dimensions:     (mass)/(length)**3

Unit name:      code_specific_energy
CGS conversion: 7.372799999999998e+19
Dimensions:     (length)**2/(time)**2

Unit name:      code_time      
CGS conversion: 647867119969847.0
Dimensions:     (time)         

Unit name:      code_magnetic  
CGS conversion: 5.141019792350707e-05
Dimensions:     sqrt((mass))/(sqrt((length))*(time))

Unit name:      code_temperature
CGS conversion: 1.0            
Dimensions:     (temperature)  

Unit name:      code_pressure  
CGS conversion: 2.0973915250863314e-06
Dimensions:     (mass)/((length)*(time)**2)

Unit name:      code_velocity  
CGS conversion: 8586501033.599192
Dimensions:     (length)/(time)

Unit name:      code_metallicity
CGS conversion: 1.0            
Dimensions:     1

fmt_tup

('code_metallicity', 1.0, '1')

Most of the time you will not have to deal with the unit registry. For example, the conversion factors to code units are stored as attributes of the dataset object:

print ("Length unit: ", ds.length_unit)
print ("Time unit: ", ds.time_unit)
print ("Mass unit: ", ds.mass_unit)
print ("Velocity unit: ", ds.velocity_unit)

Length unit:  128.0 Mpccm/h
Time unit:  647867119969847.0 s
Mass unit:  4.890451593516587e+50 g
Velocity unit:  8586501033.599192 cm/s

Conversion factors will be supplied in CGS by default. We can also ask what the conversion factors are in code units.

print ("Length unit: ", ds.length_unit.in_units('code_length'))
print ("Time unit: ", ds.time_unit.in_units('code_time'))
print ("Mass unit: ", ds.mass_unit.in_units('code_mass'))
print ("Velocity unit: ", ds.velocity_unit.in_units('code_velocity'))

Length unit:  1.0 code_length
Time unit:  1.0 code_time
Mass unit:  1.0 code_mass
Velocity unit:  1.0 code_velocity

as expected, all the conversion factors are unity in code units.

We can also play with unit conversions on ds.domain_width. First, we see for enzo how code length units are defined relative to the domain width:

ds.domain_width

YTArray([1., 1., 1.]) code_length

ds.domain_width.in_cgs()

YTArray([5.55517285e+26, 5.55517285e+26, 5.55517285e+26]) cm

ds.domain_width.in_units('Mpccm/h')

YTArray([128., 128., 128.]) Mpccm/h

Comoving units¶

This last example uses a cosmological unit. In english, I asked for the domain width in comoving megaparsecs, scaled as if the hubble constant were 100 km/s/Mpc. Although $h$ isn't really a unit, yt treats it as one for the purposes of the unit system.

As an aside, Darren Croton's research note on the history, use, and interpretation of $h$ as it appears in the astronomical literature is pretty much required reading for anyone who has to deal with factors of $h$ every now and then.

In yt, comoving length unit symbols are named following the pattern (length symbol)cm, i.e. pccm for comoving parsec or mcm for a comoving meter. A comoving length unit is different from the normal length unit by a factor of $(1+z)$:

z = ds.current_redshift
 
print (ds.quan(1, 'Mpc')/ds.quan(1, 'Mpccm'))
print (1+z)

1.0013930880640793 dimensionless
1.0013930880640796

As we saw before, $h$ is treated like any other unit symbol. It has dimensionless units, just like a scalar:

print (ds.quan(1, 'Mpc')/ds.quan(1, 'Mpc/h'))
print (ds.hubble_constant)

0.71 dimensionless
0.71

These units can be used in readily used in plots and anywhere a length unit is appropriate in yt.

slc = yt.SlicePlot(ds, 0, 'density', width=(128, 'Mpccm/h'))
slc.set_figure_size(6)

The unit registry¶

When you create a YTArray without referring to a unit registry, yt uses the default unit registry, which does not include code units or comoving units.

from yt import YTQuantity

a = YTQuantity(3, 'cm')

print (a.units.registry.keys())

dict_keys(['g', 's', 'K', 'radian', 'dyne', 'erg', 'esu', 'gauss', 'degC', 'statA', 'statV', 'statohm', 'm', 'J', 'W', 'Hz', 'N', 'C', 'A', 'T', 'Pa', 'V', 'ohm', 'ft', 'mile', 'degF', 'R', 'lbf', 'lbm', 'atm', 'h', 'dimensionless', 'min', 'hr', 'day', 'd', 'yr', 'c', 'Msun', 'msun', 'Rsun', 'rsun', 'R_sun', 'r_sun', 'Lsun', 'Tsun', 'Zsun', 'Mjup', 'Mearth', 'AU', 'au', 'ly', 'pc', 'degree', 'arcmin', 'arcsec', 'mas', 'hourangle', 'steradian', 'lat', 'lon', 'eV', 'amu', 'angstrom', 'Jy', 'counts', 'photons', 'me', 'mp', 'mol', 'Sv', 'rayleigh', 'solMass', 'solRad', 'solLum', 'dyn', 'sr', 'rad', 'deg', 'Fr', 'G', 'Angstrom', 'statC', 'm_pl', 'l_pl', 't_pl', 'T_pl', 'q_pl', 'E_pl', 'm_geom', 'l_geom', 't_geom', 'R_earth', 'r_earth', 'R_jup', 'r_jup', 'fm', 'pm', 'nm', 'um', 'mm', 'cm', 'km', 'Mm', 'kpc', 'Mpc', 'Gpc', 'mg', 'kg', 'fs', 'ps', 'ns', 'ms', 'kyr', 'Myr', 'Gyr', 'keV', 'MeV', 'GeV', 'uG'])

When a dataset is loaded, yt infers conversion factors from the internal simulation unit system to the CGS unit system. These conversion factors are stored in a unit_registry along with conversion factors to the other known unit symbols. For the cosmological Enzo dataset we loaded earlier, we can see there are a number of additional unit symbols not defined in the default unit lookup table:

print (sorted([k for k in ds.unit_registry.keys() if k not in a.units.registry.keys()]))

['AUcm', 'Mpccm', 'Ppc', 'Tpc', 'a', 'aucm', 'code_density', 'code_length', 'code_magnetic', 'code_mass', 'code_metallicity', 'code_pressure', 'code_specific_energy', 'code_temperature', 'code_time', 'code_velocity', 'mcm', 'pccm', 'unitary']

Since code units do not appear in the default unit symbol lookup table, one must explicitly refer to a unit registry when creating a YTArray to be able to convert to the unit system of a simulation.

To make this as clean as possible, there are array and quantity-creating convenience functions attached to the Dataset object:

ds.arr()
ds.quan()

These functions make it straightforward to create arrays and quantities that can be converted to code units or comoving units. For example:

a = ds.quan(3, 'code_length')

print (a)
print (a.in_cgs())
print (a.in_units('Mpccm/h'))

3.0 code_length
1.666551855079326e+27 cm
384.0 Mpccm/h

b = ds.arr([3, 4, 5], 'Mpccm/h')
print (b)
print (b.in_cgs())

[3. 4. 5.] Mpccm/h
[1.30199364e+25 1.73599152e+25 2.16998939e+25] cm

Overriding Code Unit Definitions¶

On occasion, you might have a dataset for a supported frontend that does not have the conversions to code units accessible (for example, Athena data) or you may want to change them outright. yt provides a mechanism so that one may provide their own code unit definitions to load, which override the default rules for a given frontend for defining code units. This is provided through the units_override dictionary. We'll use an example of an Athena dataset. First, a call to load without units_override:

ds1 = yt.load("MHDSloshing/virgo_low_res.0054.vtk")
print (ds1.length_unit)
print (ds1.mass_unit)
print (ds1.time_unit)
sp1 = ds1.sphere("c",(0.1,"unitary"))
print (sp1["density"])

1.0 cm
1.0 g
1.0 s
[0.05134981 0.05134912 0.05109047 ... 0.14608461 0.14489453 0.14385277] g/cm**3

This is a galaxy cluster dataset, so it is not likely that the units of density are correct. We happen to know that the unit definitions are different, so we can override the units:

units_override = {"length_unit":(1.0,"Mpc"),
                  "time_unit":(1.0,"Myr"),
                  "mass_unit":(1.0e14,"Msun")}

units_override can take the following keys:

length_unit
time_unit
mass_unit
magnetic_unit
temperature_unit

and the associated values can be (value, unit) tuples, YTQuantities, or floats (in the latter case they are assumed to have the corresponding cgs unit).

ds2 = yt.load("MHDSloshing/virgo_low_res.0054.vtk", units_override=units_override)
print (ds2.length_unit)
print (ds2.mass_unit)
print (ds2.time_unit)
sp2 = ds2.sphere("c",(0.1,"unitary"))
print (sp2["density"])

1.0 Mpc
100000000000000.0 Msun
1.0 Myr
[3.47531683e-28 3.47527018e-28 3.45776515e-28 ... 9.88689766e-28
 9.80635384e-28 9.73584863e-28] g/cm**3

This option should be used very carefully, and only if you know that the dataset does not provide units or that the unit definitions generated are incorrect for some reason.

Comoving units and code units¶

Code units¶

Comoving units¶

The unit registry¶

Overriding Code Unit Definitions¶

Units for Cosmological Datasets¶