Units in EXP and pyEXP

Attention

The Units interface described below is available in the EXP devel branch only! You will need to git checkout devel and recompile to use this feature.

Overview

The EXP N-body code assumes the gravitational constant is unity and that mass, position and velocity units can have any value defined by the user. In other words, the EXP N-body code does not enforce any particular unit set consistent with G=1.

pyEXP will also accept mass, position, time, and velocity with any unit set defined by the imported simulation as well as an arbitrary value of the gravitational constant. Explicit units types and the gravitational constant G may be set at construction time as we will describe below in The EXP unit schema. pyEXP requires the user to set units explicitly when coefficient sets are constructed. There is no default. See pyEXP Tutorial for an end-to-end example of coefficient computation.

All unit information is written into the EXP coefficient files as HDF5 attributes. Also see The HDF5 specification for details. You may add units to previously computed coefficient files using the script described in Backward compatibility.

The EXP unit schema

EXP requires one of the following two sequences of 4 unit types:

  1. Mass, Length, Time, gravitational constant (G)

  2. Mass, Length, Velocity, gravitational constant (G)

Other combinations are possible in principle, such as Mass, Length, Velocity and Time. However, this would require an automatic deduction of the gravitational constant. EXP does not current know about physical unit conversion. This feature may be added in the future.

Each separate unit is defined as a tuple which takes the form:

('unit type', 'unit name', <float value>)

The type and name strings are checked against the allowed values as follows:

  • The unit type is one of ‘length’, ‘mass’, ‘time’, ‘velocity’ or ‘G’. The internal parser will accept a variety of aliases for these such as ‘Length’, ‘Len’, ‘len’, ‘l’, ‘L’ for ‘length’; ‘Mass’, ‘m’, ‘M’ for ‘mass’; ‘Time’, ‘t’, ‘T’ for ‘time’, ‘vel’, ‘Vel’, ‘Velocity’, ‘v’, ‘V’ for ‘velocity’; and ‘Grav’, ‘grav’, ‘grav_constant’, ‘Grav_constant’, ‘gravitational_constant’, ‘Gravitational_constant’ for ‘G’.

  • The unit name is one of the usual unit names for each of the unit type. The allowed list is a subset of the standard astropy strings. For example, ‘pc’ or ‘kpc’ for ‘length’. There is also a ‘none’ type which implies no assigned physical units.

  • The floating value defines the number of each unit for the type. The value can be any valid floating-point number.

The allowed types and names may be queried interactively in pyEXP using the getAllowedUnitTypes() and getAllowedUnitNames(type), see The Units interface. For development reference, these allowed strings are defined in expui/UnitValidator.cc in the EXP source tree.

As an example, a Gadget user might define mass units as:

('mass', 'Msun', 1.0e10)

The full units specification is a list of tuples that includes one of the four sequences. In C++, the list is a std::vector of tuples.

As an example, all native EXP simulations have the following unit lists:

Units in pyEXP

The pyEXP.basis classes will use the units specified by the user in the pyEXP.coefs class created from simulation snapshots (see pyEXP Tutorial) or read by the coefficient factory from an existing HDF5 file to produce accelerations using the value of the gravitational constant and length units provided by the user.

The Units interface

The pyEXP user interface includes two member functions for explicitly setting and removing units as part of the Coef class. For setting units, we have:

setUnits(type, unit, value)
setUnits(list)
removeUnits(type)
getAllowedUnitTypes()
getAllowedTypeAliases(type)
getAllowedUnitName(type)

where type and unit are strings and value is a float. The list is a list of tuples of (name, unit, value). The last three members return the list of unit types, the recognized aliases for each type, and the allowed unit names for each unit type, respectively.

For an example, suppose you are making a set of coefficients for a simulation with default Gadget units. Say your coefficients instance is called coefs. The following command will register the unit set:

coefs.setUnits([ ('mass', 'Msun', 1e10), ('length', 'kpc', 1.0),
('velocity', 'km/s', 1.0), ('G', 'mixed', 43007.1) ])

These units will be in the HDF5 that you create using coefs.WriteH5Coefs('filename'). You can query, for example, the allowed ‘mass’ units with the call coefs.getAllowedUnitNames('mass') which returns:

['gram', 'earth_mass', 'solar_mass', 'kilograms', 'kg', 'g', 'Mearth', 'Msun', 'None', 'none']

The allowed mass units for ‘G’ are:

['gram', 'earth_mass', 'solar_mass', 'kilograms', 'kg', 'g', 'Mearth', 'Msun', 'None', 'none']

A quick note: ‘mixed’ is an allowed alias when dealing with gravitational constants that have physical units. You can see all unit types with getAllowedUnitTypes(); this returns ['mass', 'length', 'time', 'velocity', 'G']. You can see the recognized aliases for each type using getAllowedTypeAliases(type). The gravitational constant is

a special case. The recognized unit names for ‘G’ are: ['','mixed', 'none', 'unitless']`.

The C++ UI echos the functions above and adds functions to retrieve units

// Add or replace a unit by specifying its unit tuple
void setUnits(const std::string name, const std::string unit, float value);
// Add or replace a unit(s) by specifying an array of unit tuple
void setUnits(std::vector<std::tuple<std::string, std::string, float>> list);
// Remove a unit tuple by type
void removeUnits(const std::string type);
// Retrieve the currently define unit set
std::vector<std::tuple<std::string, std::string, float>> getUnits();
// Get a list of unit types and their aliases
std::vector<std::string> getAllowedTypes();
// Get a list aliases for each type
std::vector<std::string> getAllowedTypeAliases(const std::string& type);
// Get a list of unit name and their aliases for a given unit type
std::vector<std::string> getAllowedUnits(const std::string& type)

and to interact with HDF files that will only be of interest to developers creating new coefficient classes.

The HDF5 specification

Note

This and the following sections assumes basic familiarity with HDF5 and h5py in particular.

EXP HDF5 coefficient files contain a full set of metadata describing their basis type, basis parameters, geometry, and units as attributes of the root (“/”) group. The snapshots group contains all of the coefficient data and metadata (such as center and rotation information) for each snapshot time in the coefficient set.

The units information is stored in the root group as dataset named “Units”. The dataset is a sequence or list of 4 tuples. Each tuple has three fields: two fixed length strings of sixteen (16) characters and a float value.

For an EXP run, the units specification appears as dataset in the root group with the following hierarchy:

  DATASET "Units" {
   DATATYPE  H5T_COMPOUND {
      H5T_STRING {
         STRSIZE 16;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_UTF8;
         CTYPE H5T_C_S1;
      } "name";
      H5T_STRING {
         STRSIZE 16;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_UTF8;
         CTYPE H5T_C_S1;
      } "unit";
      H5T_IEEE_F32LE "value";
   }
   DATASPACE  SIMPLE { ( 4 ) / ( 4 ) }
   DATA {
   (0): {
         "G",
         "none",
         1
      },
   (1): {
         "length",
         "none",
         1
      },
   (2): {
         "mass",
         "none",
         1
      },
   (3): {
         "time",
         "none",
         1
      }
   }
}

Backward compatibility

A units attribute list could be easily added to existing HDF5 files using h5py using the schema described above. For example, assuming that you already have an open HDF5 coefficient file named f, the units described in the scheme above could be added to f with the following code:

import h5py
import numpy as np

# Define the compound datatype with fixed-length ASCII strings and a float32
dt = np.dtype([
    ('name', 'S16'),   # Fixed-length ASCII string of 16 bytes
    ('unit', 'S16'),   # Fixed-length ASCII string of 16 bytes
    ('value', '<f4')   # Little-endian 32-bit float
])

# Create the dataset data as a numpy array
data = np.array([
    (b'G',       b'none', 1.0),
    (b'length',  b'none', 1.0),
    (b'mass',    b'none', 1.0),
    (b'time',    b'none', 1.0),
], dtype=dt)

# Open the HDF5 coefficient file and add the dataset
with h5py.File('outcoef.disk.runtag', 'a') as f:
    # 'a' mode opens file for read/write, creates if not exist

    # Create dataset 'Units'
    if 'Units' in f:
        # Optional: remove existing dataset if you want to overwrite
        del f['Units']
    dset = f.create_dataset('Units', data=data, dtype=dt)

# File is automatically closed on leaving the 'with' block, flush and
# update saved data on disk

You can update the data array in the example above to match your unit set.