diff --git a/index.rst b/index.rst index 5f395b7..848fe8e 100644 --- a/index.rst +++ b/index.rst @@ -89,6 +89,7 @@ EXP concepts topics/bfetheory topics/codeintro topics/timesseries + topics/units topics/multistep topics/centering topics/yamlconfig @@ -106,6 +107,9 @@ mathematics used in EXP. :doc:`topics/timesseries` Descriptions for methods to analyze BFE time series +:doc:`topics/units` + Physical units in EXP + :doc:`topics/multistep` A quick review of EXP's ODE solver and multi time step ladder diff --git a/intro/overview.rst b/intro/overview.rst index 688e5e5..1f48132 100644 --- a/intro/overview.rst +++ b/intro/overview.rst @@ -96,6 +96,11 @@ Here's the code for computing coefficients:: # coef = basis.createFromArray(data[:,0], data[:,1:4], time=3.0) + # Add your simulation units (we assume a unit-free simulation here) + # + coef.setUnits([('mass', 'none', 1.0), ('length', 'none', 1.0), + ('time', 'none', 1.0), ('G', 'none', 1.0)]) + # Make an HDF5 file # coefs = pyEXP.coefs.Coefs.makecoefs(coef) @@ -144,6 +149,11 @@ More detailed information on YAML and config parameters is available in the :ref:`What is YAML?` and :ref:`How to visualize the BFE bases used to make your coefficients` pages. +The unit system in EXP is described in more detail in :ref:`Units in +EXP and pyEXP`. For this example, we assume a unit-free +simulation but it's easy to add any unit system that is convenient for +you. + pyEXP is then ready to make the coefficients from your phase-space data. This example assumes that the mass and positions of your particles are in columns 1, 2, 3, 4 of the file and that the positions diff --git a/topics/units.rst b/topics/units.rst new file mode 100644 index 0000000..a610ddb --- /dev/null +++ b/topics/units.rst @@ -0,0 +1,291 @@ +.. role:: python(code) + :language: python + :class: highlight + +.. _units: + +Units in EXP and pyEXP +====================== + +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 :ref:`unit-schema`. pyEXP **requires** the +user to set units explicitly when coefficient sets are +constructed. There is no default. See :ref:`intro-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 :ref:`hdf5-unit-attributes` for details. You +may add units to previously computed coefficient files using the +script described in :ref:`units_old`. + +.. _unit-schema: + +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', ) + +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 :python:`getAllowedUnitTypes()` and +:python:`getAllowedUnitNames(type)`, see :ref:`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 :cpp:any:`std::vector<>` of +tuples. + +As an example, all native EXP simulations have the following unit +lists: + + .. code-block:: python + [('mass', 'none', 1.0f), ('length', 'none', 1.0f), ('time', 'none', 1.0f), ('G', 'none', 1.0f)] + + +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 +:ref:`intro-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. + + +.. _units-interface: + +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: + + .. code-block:: python + + 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: + + .. code-block:: python + + 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 +:python:`coefs.WriteH5Coefs('filename')`. You can query, for example, +the allowed 'mass' units with the call +:python:`coefs.getAllowedUnitNames('mass')` which returns: + + .. code-block:: python + + ['gram', 'earth_mass', 'solar_mass', 'kilograms', 'kg', 'g', 'Mearth', 'Msun', 'None', 'none'] + +The allowed mass units for 'G' are: + + .. code-block:: python + + ['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 :python:`getAllowedUnitTypes()`; this returns +:python:`['mass', 'length', 'time', 'velocity', 'G']`. You can see +the recognized aliases for each type using +:python:`getAllowedTypeAliases(type)`. The gravitational constant is + a special case. The recognized unit names for 'G' are: :python:`['','mixed', 'none', 'unitless']``. + +The C++ UI echos the functions above and adds functions to retrieve +units + + .. code-block:: C++ + + // 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> list); + // Remove a unit tuple by type + void removeUnits(const std::string type); + // Retrieve the currently define unit set + std::vector> getUnits(); + // Get a list of unit types and their aliases + std::vector getAllowedTypes(); + // Get a list aliases for each type + std::vector getAllowedTypeAliases(const std::string& type); + // Get a list of unit name and their aliases for a given unit type + std::vector getAllowedUnits(const std::string& type) + +and to interact with HDF files that will only be of interest to +developers creating new coefficient classes. + + +.. _hdf5-unit-attributes: + +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 + } + } + } + + +.. _units_old: + +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: + + .. code-block:: python + + 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', '