Working with SI-units

Goal

Learn how to work with SINumber and SIArray objects which represent physical quantities, i.e. one or more floating point numbers with an associated unit.

The si-units module

Most interfaces in FeOs use dimensioned quantities as input. For example, to define a thermodynamic state at given temperature, pressure and amount of substance, all of these properties have to be multiplied by an apropriate unit before we can call the function that creates the state.

FeOs uses the quantity crate which generates the si_units Python module. To have consistency throughout the ecosystem, we recommend importing the package as import si_units as si

import si_units as si
import numpy as np

The si module contains units according to the Standard Interational System of Units (SI), constants and prefixes. A scalar floating point number multiplied or divided by a unit has the ``SIObject`` data type.

temperature = 300.15 * si.KELVIN
print(f"temperature : {temperature}")
print(f"data type   : {type(temperature)}")

temperature : 300.15 K
data type   : <class 'si_units._core.SIObject'>

The representation of a SIObject (e.g. using print) uses a prefix (e.g. k for kilo) so that the numerical value is convenient. For example 1000 g will be represented as 1 kg. Internally, all values are stored with respect to the basic SI unit, i.e. METER, KILOGRAM, SECOND, AMPERE, MOL, KELVIN, and CANDELA.

1000 * si.GRAM

$1\,\mathrm{kg}$

# volume of an ideal gas
t = 300.15 * si.KELVIN
p = 0.5 * si.BAR
n = 1.5 * si.MOL
v = n * si.RGAS * t / p
print(f"v = {v}")

v = 0.07486757864516083 m³

We can use division to transform a SIObject to a float for the desired unit:

# temperature in Celsius
t_c = t / si.CELSIUS
print(f"t = {t_c:3.2f} °C")
print(f"data type = {type(t_c)}")

t = 27.00 °C
data type = <class 'float'>

Mathematical operations and interaction with numpy

The SINumber object supports some mathematical operations, i.e. we can

• add or subtract two objects, if they have the same unit,

• multiply or divide two objects, which returns a floating point number if they have the same unit,

• raise an object to a power, and

• take the square and cubic root, which only works if the units have the correct exponents to allow for the operation.

# addition
pressure = 2.5 * si.BAR + 15_000 * si.PASCAL
print('pressure          :', pressure)

# subtraction
temperature = 543.15 * si.KELVIN - 230.0 * si.CELSIUS
print('temperature       :', temperature)

# division
velocity = 360_000 * si.METER / si.HOUR
print('velocity          :', velocity)

# division as transformation to target unit
v_cm_minute = velocity / (si.CENTI * si.METER / si.MINUTE) # this is a floating point number
print(f'velocity (cm/min) : {v_cm_minute:8.1f}')

# power
acceleration = 9.81 * si.METER / si.SECOND**2
print('acceleration      :', acceleration)

pressure          : 265 kPa
temperature       : 40 K
velocity          : 100  m/s
velocity (cm/min) : 600000.0
acceleration      : 9.81 m s^-2

For the square and cubic root we can use the numpy functions, numpy.sqrt and numpy.cbrt, respectively, or use the methods provided by si_units.

# square root
speed = np.sqrt(si.METER**2 / si.SECOND**2)
print('speed             :', speed)

# cubic root
box_length = (27_000 * si.ANGSTROM**3).cbrt()
print('length            :', box_length)

# both numpy and methods of SINumbers work
assert(np.sqrt(si.METER**2 / si.SECOND**2) == (si.METER**2 / si.SECOND**2).sqrt())

speed             : 1  m/s
length            : 3 nm

The mathematical functions from the math module do not work because it is written in C and SINumbers are not compatible with the C function arguments. If you are having trouble with errors, consider transforming your properties to floats (by division with the proper unit) and then multiplying the correct unit to the result.

from math import sqrt

try:
    sqrt(si.METER**2 / si.SECOND**2)
except Exception as e:
    print("ERROR:", e)

ERROR: must be real number, not si_units._core.SIObject

If you try to perform arithmetic operations with SINumbers that are not allowed (e.g. adding two properties with different units), an error is raised.

si.PASCAL + si.KELVIN

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Input In [9], in <module>
----> 1 si.PASCAL + si.KELVIN

RuntimeError: Inconsistent units Pa and K

We can enforce correct units in our interfaces using the has_unit method:

MOL_M3 = si.MOL / si.METER**3 # we can create own units for future use

def ideal_gas_pressure(density, temperature):
    if not density.has_unit(MOL_M3):
        raise ValueError("Please provide the molar density, e.g. in units of mol/m³.")
    else:
        return density * temperature * si.RGAS

try:
    p1 = ideal_gas_pressure(0.5 * si.KILO * si.MOL / si.METER**3, 350 * si.KELVIN)
    print('p1 = ', p1)
    p2 = ideal_gas_pressure(0.5 * si.KILOGRAM / si.METER**3, 350 * si.KELVIN)
    print('p2 = ', p2)
except Exception as e:
    print("ERROR:", e)

p1 =  1.4550309581768168 MPa
ERROR: Please provide the molar density, e.g. in units of mol/m³.

Arrays of quantities

The “number” part of the SIObject can be any type that supports the necessary operations. Aside from floats, those are, e.g., numpy arrays or tensors.

ps = np.array([1.0, 2.0]) * si.BAR
print(f"pressures = {ps}")
print(f"data type = {type(ps)}")

pressures = array([100000., 200000.]) Pa
data type = <class 'si_units._core.SIObject'>

The most important functionalities of array or tensor datatypes are also implemented for SIObject. For example, you can index into an array of dimension one.

print(ps[1])

200 kPa

There are several useful numpy methods to create arrays. Most of them can be simply multiplied by units to convert them to SIObject objects. Some of these functions, such as linspace and logspace can directly be constructed using si.linspace and si.logspace, respectively.

ps_np = np.linspace(1, 2, 5) * si.BAR
ps_si = si.linspace(1 * si.BAR, 2 * si.BAR, 5)
print(f"pressures (numpy) = {ps_np}")
print(f"pressures (si)    = {ps_si}")

pressures (numpy) = array([100000., 125000., 150000., 175000., 200000.]) Pa
pressures (si)    = array([100000., 125000., 150000., 175000., 200000.]) Pa

Just like scalars, division by a unit that matches the property stored in an array yields a numpy ndarray.

ps_mpa = ps_np / (si.MEGA*si.PASCAL)
print(f"pressures = {ps_mpa} MPa")
print(f"data type = {type(ps_mpa)}")

pressures = [0.1   0.125 0.15  0.175 0.2  ] MPa
data type = <class 'numpy.ndarray'>

Derived units, constants and prefixes

In conjunction to the base units, the si module also exports derived units (e.g. HOUR = 3600 * SECOND), constants and prefixes (such as KILO or FEMTO). Of course, we could multiply by a floating point number instead of using prefixes (e.g. cm = 1e-2 * METER vs. cm = CENTI * METER) but we think using prefixes make the code more readable.

For a complete overview of exported units constants and prefixes, take a look at the documentation of the ``si_units` Python package <https://itt-ustutt.github.io/quantity/base/>`__.

# The seven constants that inform the base units
print('Hyperfine transition frequency of Cs:', si.DVCS)
print('Speed of light                      :', si.CLIGHT)
print('Planck constant                     :', si.PLANCK)
print('Elementary charge                   :', si.QE)
print('Boltzmann constant                  :', si.KB)
print('Avogradro constant                  :', si.NAV)
print('Luminous efficacy                   :', si.KCD)

Hyperfine transition frequency of Cs: 9.19263177 GHz
Speed of light                      : 2.99792458e8  m/s
Planck constant                     : 6.62607015e-34  Js
Elementary charge                   : 1.602176634e-19 C
Boltzmann constant                  : 1.380649e-23  J/K
Avogradro constant                  : 6.02214076e23 mol^-1
Luminous efficacy                   : 683 lm/W

# Derived constants
print('Gravitational constant:', si.G)
print('Ideal gas constant    :', si.RGAS)

Gravitational constant: 6.6743e-11 m³/kg/s²
Ideal gas constant    : 8.31446261815324  J/mol/K

Summary

In feos, most interfaces use dimensioned units in form of SIObject from the si-units Python library. This enables us to write functions that can check for proper units and thus circumvent unit errors that could occur, e.g. providing mass densities instead of molar densities.

In this example, we learned how to create and use these objects in conjunction with methods provided by numpy.

Concluding remkars

Hopefully you found this example helpful. If you have comments, critique or feedback, please let us know and consider opening an issue on github.