Skip to content

Commit 6e3a24b

Browse files
The9CatThe9Cat
authored
Merge pull request #38 from EXP-code/AddUnits
Add some units documentation
2 parents fd7fd2d + 0d28624 commit 6e3a24b

File tree

3 files changed

+305
-0
lines changed

3 files changed

+305
-0
lines changed

‎index.rst‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -89,6 +89,7 @@ EXP concepts
8989
topics/bfetheory
9090
topics/codeintro
9191
topics/timesseries
92+
topics/units
9293
topics/multistep
9394
topics/centering
9495
topics/yamlconfig
@@ -106,6 +107,9 @@ mathematics used in EXP.
106107
:doc:`topics/timesseries`
107108
Descriptions for methods to analyze BFE time series
108109

110+
:doc:`topics/units`
111+
Physical units in EXP
112+
109113
:doc:`topics/multistep`
110114
A quick review of EXP's ODE solver and multi time step ladder
111115

‎intro/overview.rst‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,11 @@ Here's the code for computing coefficients::
9696
#
9797
coef = basis.createFromArray(data[:,0], data[:,1:4], time=3.0)
9898

99+
# Add your simulation units (we assume a unit-free simulation here)
100+
#
101+
coef.setUnits([('mass', 'none', 1.0), ('length', 'none', 1.0),
102+
('time', 'none', 1.0), ('G', 'none', 1.0)])
103+
99104
# Make an HDF5 file
100105
#
101106
coefs = pyEXP.coefs.Coefs.makecoefs(coef)
@@ -144,6 +149,11 @@ More detailed information on YAML and config parameters is available
144149
in the :ref:`What is YAML?<yamlconfig>` and :ref:`How to visualize the
145150
BFE bases used to make your coefficients<visualizing-bases>` pages.
146151

152+
The unit system in EXP is described in more detail in :ref:`Units in
153+
EXP and pyEXP<units>`. For this example, we assume a unit-free
154+
simulation but it's easy to add any unit system that is convenient for
155+
you.
156+
147157
pyEXP is then ready to make the coefficients from your phase-space
148158
data. This example assumes that the mass and positions of your
149159
particles are in columns 1, 2, 3, 4 of the file and that the positions

‎topics/units.rst‎

Lines changed: 291 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,291 @@
1+
.. role:: python(code)
2+
:language: python
3+
:class: highlight
4+
5+
.. _units:
6+
7+
Units in EXP and pyEXP
8+
======================
9+
10+
Overview
11+
--------
12+
13+
The EXP N-body code assumes the gravitational constant is unity and
14+
that mass, position and velocity units can have any value defined by
15+
the user. In other words, the EXP N-body code does not enforce any
16+
particular unit set consistent with ``G=1``.
17+
18+
pyEXP will also accept mass, position, time, and velocity with any
19+
unit set defined by the imported simulation as well as an arbitrary
20+
value of the gravitational constant. Explicit units types and the
21+
gravitational constant ``G`` may be set at construction time as we
22+
will describe below in :ref:`unit-schema`. pyEXP **requires** the
23+
user to set units explicitly when coefficient sets are
24+
constructed. There is no default. See :ref:`intro-pyEXP-tutorial` for
25+
an end-to-end example of coefficient computation.
26+
27+
All unit information is written into the EXP coefficient files as HDF5
28+
attributes. Also see :ref:`hdf5-unit-attributes` for details. You
29+
may add units to previously computed coefficient files using the
30+
script described in :ref:`units_old`.
31+
32+
.. _unit-schema:
33+
34+
The EXP unit schema
35+
-------------------
36+
37+
EXP requires one of the following two sequences of 4 unit types:
38+
39+
1. Mass, Length, Time, gravitational constant (G)
40+
41+
2. Mass, Length, Velocity, gravitational constant (G)
42+
43+
44+
Other combinations are possible in principle, such as Mass, Length,
45+
Velocity and Time. However, this would require an automatic deduction
46+
of the gravitational constant. EXP does not current know about
47+
physical unit conversion. This feature may be added in the future.
48+
49+
Each separate unit is defined as a ``tuple`` which takes the form::
50+
51+
('unit type', 'unit name', <float value>)
52+
53+
The type and name strings are checked against the allowed values as follows:
54+
55+
- The ``unit type`` is one of 'length', 'mass', 'time', 'velocity' or
56+
'G'. The internal parser will accept a variety of aliases for these
57+
such as 'Length', 'Len', 'len', 'l', 'L' for 'length'; 'Mass', 'm',
58+
'M' for 'mass'; 'Time', 't', 'T' for 'time', 'vel',
59+
'Vel', 'Velocity', 'v', 'V' for 'velocity'; and 'Grav', 'grav',
60+
'grav_constant', 'Grav_constant', 'gravitational_constant',
61+
'Gravitational_constant' for 'G'.
62+
63+
- The ``unit name`` is one of the usual unit names for each of the
64+
``unit type``. The allowed list is a subset of the standard
65+
`astropy` strings. For example, 'pc' or 'kpc' for 'length'. There
66+
is also a 'none' type which implies no assigned physical units.
67+
68+
- The floating value defines the number of each unit for the type.
69+
The value can be any valid floating-point number.
70+
71+
The allowed types and names may be queried interactively in pyEXP
72+
using the :python:`getAllowedUnitTypes()` and
73+
:python:`getAllowedUnitNames(type)`, see :ref:`units-interface`. For
74+
development reference, these allowed strings are defined in
75+
``expui/UnitValidator.cc`` in the EXP source tree.
76+
77+
As an example, a Gadget user might define mass units as::
78+
79+
('mass', 'Msun', 1.0e10)
80+
81+
The full units specification is a list of tuples that includes one of
82+
the four sequences. In C++, the list is a :cpp:any:`std::vector<>` of
83+
tuples.
84+
85+
As an example, all native EXP simulations have the following unit
86+
lists:
87+
88+
.. code-block:: python
89+
[('mass', 'none', 1.0f), ('length', 'none', 1.0f), ('time', 'none', 1.0f), ('G', 'none', 1.0f)]
90+
91+
92+
Units in pyEXP
93+
--------------
94+
95+
The ``pyEXP.basis`` classes will use the units specified by the user
96+
in the ``pyEXP.coefs`` class created from simulation snapshots (see
97+
:ref:`intro-pyEXP-tutorial`) or read by the coefficient factory from
98+
an existing HDF5 file to produce accelerations using the value of the
99+
gravitational constant and length units provided by the user.
100+
101+
102+
.. _units-interface:
103+
104+
The Units interface
105+
-------------------
106+
107+
The `pyEXP` user interface includes two member functions for explicitly
108+
setting and removing units as part of the `Coef` class. For setting
109+
units, we have:
110+
111+
.. code-block:: python
112+
113+
setUnits(type, unit, value)
114+
setUnits(list)
115+
removeUnits(type)
116+
getAllowedUnitTypes()
117+
getAllowedTypeAliases(type)
118+
getAllowedUnitName(type)
119+
120+
where ``type`` and ``unit`` are strings and ``value`` is a float. The
121+
list is a list of tuples of ``(name, unit, value)``. The last three
122+
members return the list of unit types, the recognized aliases for each
123+
type, and the allowed unit names for each unit type, respectively.
124+
125+
For an example, suppose you are making a set of coefficients for a
126+
simulation with default Gadget units. Say your coefficients instance
127+
is called ``coefs``. The following command will register the unit set:
128+
129+
.. code-block:: python
130+
131+
coefs.setUnits([ ('mass', 'Msun', 1e10), ('length', 'kpc', 1.0),
132+
('velocity', 'km/s', 1.0), ('G', 'mixed', 43007.1) ])
133+
134+
These units will be in the HDF5 that you create using
135+
:python:`coefs.WriteH5Coefs('filename')`. You can query, for example,
136+
the allowed 'mass' units with the call
137+
:python:`coefs.getAllowedUnitNames('mass')` which returns:
138+
139+
.. code-block:: python
140+
141+
['gram', 'earth_mass', 'solar_mass', 'kilograms', 'kg', 'g', 'Mearth', 'Msun', 'None', 'none']
142+
143+
The allowed mass units for 'G' are:
144+
145+
.. code-block:: python
146+
147+
['gram', 'earth_mass', 'solar_mass', 'kilograms', 'kg', 'g', 'Mearth', 'Msun', 'None', 'none']
148+
149+
150+
A quick note: 'mixed' is an allowed alias when dealing with
151+
gravitational constants that have physical units. You can see all
152+
unit types with :python:`getAllowedUnitTypes()`; this returns
153+
:python:`['mass', 'length', 'time', 'velocity', 'G']`. You can see
154+
the recognized aliases for each type using
155+
:python:`getAllowedTypeAliases(type)`. The gravitational constant is
156+
a special case. The recognized unit names for 'G' are: :python:`['','mixed', 'none', 'unitless']``.
157+
158+
The C++ UI echos the functions above and adds functions to retrieve
159+
units
160+
161+
.. code-block:: C++
162+
163+
// Add or replace a unit by specifying its unit tuple
164+
void setUnits(const std::string name, const std::string unit, float value);
165+
// Add or replace a unit(s) by specifying an array of unit tuple
166+
void setUnits(std::vector<std::tuple<std::string, std::string, float>> list);
167+
// Remove a unit tuple by type
168+
void removeUnits(const std::string type);
169+
// Retrieve the currently define unit set
170+
std::vector<std::tuple<std::string, std::string, float>> getUnits();
171+
// Get a list of unit types and their aliases
172+
std::vector<std::string> getAllowedTypes();
173+
// Get a list aliases for each type
174+
std::vector<std::string> getAllowedTypeAliases(const std::string& type);
175+
// Get a list of unit name and their aliases for a given unit type
176+
std::vector<std::string> getAllowedUnits(const std::string& type)
177+
178+
and to interact with HDF files that will only be of interest to
179+
developers creating new coefficient classes.
180+
181+
182+
.. _hdf5-unit-attributes:
183+
184+
The HDF5 specification
185+
----------------------
186+
187+
.. note:: This and the following sections assumes basic familiarity
188+
with HDF5 and `h5py` in particular.
189+
190+
EXP HDF5 coefficient files contain a full set of metadata describing
191+
their basis type, basis parameters, geometry, and units as attributes
192+
of the root ("/") group. The ``snapshots`` group contains all of the
193+
coefficient data and metadata (such as center and rotation
194+
information) for each snapshot time in the coefficient set.
195+
196+
The units information is stored in the root group as dataset named
197+
"Units". The dataset is a sequence or list of 4 tuples. Each tuple
198+
has three fields: two fixed length strings of sixteen (16) characters
199+
and a float value.
200+
201+
For an EXP run, the units specification appears as dataset in the root
202+
group with the following hierarchy::
203+
204+
DATASET "Units" {
205+
DATATYPE H5T_COMPOUND {
206+
H5T_STRING {
207+
STRSIZE 16;
208+
STRPAD H5T_STR_NULLTERM;
209+
CSET H5T_CSET_UTF8;
210+
CTYPE H5T_C_S1;
211+
} "name";
212+
H5T_STRING {
213+
STRSIZE 16;
214+
STRPAD H5T_STR_NULLTERM;
215+
CSET H5T_CSET_UTF8;
216+
CTYPE H5T_C_S1;
217+
} "unit";
218+
H5T_IEEE_F32LE "value";
219+
}
220+
DATASPACE SIMPLE { ( 4 ) / ( 4 ) }
221+
DATA {
222+
(0): {
223+
"G",
224+
"none",
225+
1
226+
},
227+
(1): {
228+
"length",
229+
"none",
230+
1
231+
},
232+
(2): {
233+
"mass",
234+
"none",
235+
1
236+
},
237+
(3): {
238+
"time",
239+
"none",
240+
1
241+
}
242+
}
243+
}
244+
245+
246+
.. _units_old:
247+
248+
Backward compatibility
249+
----------------------
250+
251+
A units attribute list could be easily added to existing HDF5 files
252+
using `h5py` using the schema described above. For example, assuming
253+
that you already have an open HDF5 coefficient file named `f`, the
254+
units described in the scheme above could be added to `f` with the
255+
following code:
256+
257+
.. code-block:: python
258+
259+
import h5py
260+
import numpy as np
261+
262+
# Define the compound datatype with fixed-length ASCII strings and a float32
263+
dt = np.dtype([
264+
('name', 'S16'), # Fixed-length ASCII string of 16 bytes
265+
('unit', 'S16'), # Fixed-length ASCII string of 16 bytes
266+
('value', '<f4') # Little-endian 32-bit float
267+
])
268+
269+
# Create the dataset data as a numpy array
270+
data = np.array([
271+
(b'G', b'none', 1.0),
272+
(b'length', b'none', 1.0),
273+
(b'mass', b'none', 1.0),
274+
(b'time', b'none', 1.0),
275+
], dtype=dt)
276+
277+
# Open the HDF5 coefficient file and add the dataset
278+
with h5py.File('outcoef.disk.runtag', 'a') as f:
279+
# 'a' mode opens file for read/write, creates if not exist
280+
281+
# Create dataset 'Units'
282+
if 'Units' in f:
283+
# Optional: remove existing dataset if you want to overwrite
284+
del f['Units']
285+
dset = f.create_dataset('Units', data=data, dtype=dt)
286+
287+
# File is automatically closed on leaving the 'with' block, flush and
288+
# update saved data on disk
289+
290+
You can update the :python:`data` array in the example above to match your
291+
unit set.

0 commit comments

Comments
 (0)