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 feos.si 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. You might find this module useful on it’s own in which case you can also install it as a separate python package. For convenient use with FeOs however, it is re-exported as feos.si.

from feos.si import *
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 ``SINumber`` data type.

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

temperature : 300.15 K
data type   : <class 'feos.si.SINumber'>

The representation of a SINumber (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 * GRAM

$1\,\mathrm{kg}$

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

v = 0.07486757864516083 m³

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

# temperature in Celsius
t_c = t / 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 * BAR + 15_000 * PASCAL
print('pressure          :', pressure)

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

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

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

# power
acceleration = 9.81 * METER / 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(METER**2 / SECOND**2)
print('speed             :', speed)

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

# both numpy and methods of SINumbers work
assert(np.sqrt(METER**2 / SECOND**2) == (METER**2 / 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(METER**2 / SECOND**2)
except Exception as e:
    print("ERROR:", e)

ERROR: must be real number, not si_units.SINumber

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.

PASCAL + KELVIN

---------------------------------------------------------------------------
PanicException                            Traceback (most recent call last)
<ipython-input-9-3b8b977148c6> in <module>
----> 1 PASCAL + KELVIN

PanicException: Inconsistent units Pa + K

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

MOL_M3 = MOL / 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 * RGAS

try:
    p1 = ideal_gas_pressure(0.5 * KILO * MOL / METER**3, 350 * KELVIN)
    print('p1 = ', p1)
    p2 = ideal_gas_pressure(0.5 * KILOGRAM / METER**3, 350 * 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

A numpy array of floating point numbers multiplied or divided by a unit has the SIArrayX data type, where the X stands for the dimension of the numpy array.

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

pressures = [100000, 200000] Pa
data type = <class 'feos.si.SIArray1'>

Currently, indexing into an array only works for arrays 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 SIArray objects. Some of these functions, such as linspace and logspace can directly be constructed using SIArray1.linspace and SIArray1.logspace, respectively.

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

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

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

ps_mpa = ps_np / (MEGA*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'>

Sometimes, we write functions that use arrays as input where it is possible that these arrays only contain one element. A numpy array with a single element qualifies for multiplication rules of a scalar and a SINumber which can cause problems.

Consider the following example:

def print_length(array):
    print(f'Array "{array}" contains {len(array)} elements.\n')

n2 = np.array([1.0, 2.0]) * MOL # data type: SIArray1
n1 = np.array([1.0]) * MOL      # data type: SINumber

try:
    print_length(n2)
    print_length(n1) # <- this fails, since n1 has data type SINumber, not SIArray1
except Exception as e:
    print('ERROR:', e)

Array "[1, 2] mol" contains 2 elements.

ERROR: object of type 'si_units.SINumber' has no len()

This example does not work, because np.array([1.0]) * MOL has data type SINumber instead of SIArray1. We can construct a SIArray1 using the base constructor:

# now it works
n1 = SIArray1(1.0*MOL)
print_length(n1)

Array "[1] mol" contains 1 elements.

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/api.html#si-base-units-and-associated-constants>`__.

# The seven constants that inform the base units
print('Hyperfine transition frequency of Cs:', DVCS)
print('Speed of light                      :', CLIGHT)
print('Planck constant                     :', PLANCK)
print('Elementary charge                   :', QE)
print('Boltzmann constant                  :', KB)
print('Avogradro constant                  :', NAV)
print('Luminous efficacy                   :', 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:', G)
print('Ideal gas constant    :', 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 SINumber and SIArrayX (where X stands for the dimnesion). 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.