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
crate which generates the si_units Python module. To have consistency throughout the ecosystem, we recommend importing the package as importsi_unitsassi
[1]: importsi_unitsassiimportnumpyasnpThe 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.
[2]: temperature=300.15*si.KELVINprint(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.
[3]: 1000*si.GRAM[3]: $1\,\mathrm{kg}$
[4]: # volume of an ideal gast=300.15*si.KELVINp=0.5*si.BARn=1.5*si.MOLv=n*si.RGAS*t/pprint(f"v = {v}")v = 0.07486757864516083 m³ We can use division to transform a SIObject to a float for the desired unit:
[5]: # temperature in Celsiust_c=t/si.CELSIUSprint(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.
[6]: # additionpressure=2.5*si.BAR+15_000*si.PASCALprint('pressure :',pressure)# subtractiontemperature=543.15*si.KELVIN-230.0*si.CELSIUSprint('temperature :',temperature)# divisionvelocity=360_000*si.METER/si.HOURprint('velocity :',velocity)# division as transformation to target unitv_cm_minute=velocity/(si.CENTI*si.METER/si.MINUTE)# this is a floating point numberprint(f'velocity (cm/min) : {v_cm_minute:8.1f}')# poweracceleration=9.81*si.METER/si.SECOND**2print('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.
[7]: # square rootspeed=np.sqrt(si.METER**2/si.SECOND**2)print('speed :',speed)# cubic rootbox_length=(27_000*si.ANGSTROM**3).cbrt()print('length :',box_length)# both numpy and methods of SINumbers workassert(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.
[8]: frommathimportsqrttry:sqrt(si.METER**2/si.SECOND**2)exceptExceptionase: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.
[9]: si.PASCAL+si.KELVIN---------------------------------------------------------------------------RuntimeError Traceback (most recent call last) CellIn[9], line 1----> 1si.PASCAL+si.KELVINRuntimeError: Inconsistent units Pa and K We can enforce correct units in our interfaces using the has_unit method:
[10]: MOL_M3=si.MOL/si.METER**3# we can create own units for future usedefideal_gas_pressure(density,temperature):ifnotdensity.has_unit(MOL_M3):raiseValueError("Please provide the molar density, e.g. in units of mol/m³.")else:returndensity*temperature*si.RGAStry: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)exceptExceptionase: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.
[11]: ps=np.array([1.0,2.0])*si.BARprint(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.
[12]: 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.
[13]: ps_np=np.linspace(1,2,5)*si.BARps_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.
[14]: 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/
>`__.
[15]: # The seven constants that inform the base unitsprint('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 [16]: # Derived constantsprint('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
.