Comoving units and code units

Notebook

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:

In [1]:
import yt
ds = yt.load('Enzo_64/DD0043/data0043')

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:

In [2]:
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_magnetic  
CGS conversion: 0.005133186931548846
Dimensions:     sqrt((mass))/(sqrt((length))*(time))

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

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

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

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

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

Unit name:      code_metallicity
CGS conversion: 1.0            
Dimensions:     1              

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

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

Unit name:      code_specific_energy
CGS conversion: 7.375871999999998e+19
Dimensions:     (length)**2/(time)**2
In [3]:
fmt_tup
Out[3]:
('code_specific_energy', 7.375871999999998e+19, '(length)**2/(time)**2')

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:

In [4]:
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:  648009786344140.2 s
Mass unit:  4.891307900838622e+50 g
Velocity unit:  8588289701.681004 cm/s

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

In [5]:
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:

In [6]:
ds.domain_width
Out[6]:
YTArray([ 1.,  1.,  1.]) code_length
In [7]:
ds.domain_width.in_cgs()
Out[7]:
YTArray([  5.55517285e+26,   5.55517285e+26,   5.55517285e+26]) cm
In [8]:
ds.domain_width.in_units('Mpccm/h')
Out[8]:
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)$:

In [9]:
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:

In [10]:
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.

In [11]:
slc = yt.SlicePlot(ds, 0, 'density', width=(128, 'Mpccm/h'))
slc.set_figure_size(6)
Out[11]:

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.

In [12]:
from yt import YTQuantity

a = YTQuantity(3, 'cm')

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

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

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:

In [13]:
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:

In [14]:
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
In [15]:
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:

In [16]:
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:

In [17]:
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).

In [18]:
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.

(3)_Comoving_units_and_code_units.ipynb; 3)_Comoving_units_and_code_units_evaluated.ipynb; 3)_Comoving_units_and_code_units.py)

Units for Cosmological Datasets

yt has additional capabilities to handle the comoving coordinate system used internally in cosmological simulations. Simulations that use comoving coordinates, all length units have three other counterparts corresponding to comoving units, scaled comoving units, and scaled proper units. In all cases ‘scaled’ units refer to scaling by the reduced Hubble parameter - i.e. the length unit is what it would be in a universe where Hubble’s parameter is 100 km/s/Mpc.

To access these different units, yt has a common naming system. Scaled units are denoted by dividing by the scaled Hubble parameter h (which is in itself a unit). Comoving units are denoted by appending cm to the end of the unit name.

Using the parsec as an example,

pc
Proper parsecs, \(\rm{pc}\).
pccm
Comoving parsecs, \(\rm{pc}/(1+z)\).
pccm/h
Comoving parsecs normalized by the scaled hubble constant, \(\rm{pc}/h/(1+z)\).
pc/h
Proper parsecs, normalized by the scaled hubble constant, \(\rm{pc}/h\).
© Copyright 2013-2018, the yt Project. Last updated on Feb 16, 2018. Created using Sphinx 1.5.2.