02 - Variables and HDF5 IO#

ReMKiT1D values of interest during simulations are almost always stored in Variables. This tutorial covers the following:

  1. Basic Variable construction

  2. Variable properties

  3. Variable class features - unit conversion, dimensions, xarray conversion, etc.

  4. Nodes and expression trees for derived variables

  5. Variable container

  6. Saving and loading ReMKiT1D hdf5 files

[1]:

from RMK_support import Variable,varAndDual,node,varFromNode,VariableContainer,Grid,loadVarContFromHDF5,loadVariableFromHDF5,writeRMKHDF5,loadFromHDF5

import numpy as np

In order to construct variables, we require a grid object. We reuse the grid we built in the first tutorial.

[2]:

grid = Grid(xGrid = 0.1 * np.ones(16),
            interpretXGridAsWidths = True,
            vGrid = 0.1 * np.ones(8),
            interpretVGridAsWidths = True,
            lMax = 3,
            )

The Variable constructor has many keyword argument options, only some of which will be explored here. The user is encouraged to explore other tutorials and examples, as well as the documentation to see how other options are used.

[3]:

a = Variable("a", # Name of the variable
             grid, # Grid on which the variable lives
             data = np.linspace(0,15,16), # Spatial grid vector - values of the variable
                                          # - by default variables live on the spatial grid
                                          # i.e. they are considered fluid
             units = "aUnits", # Optional name of the variable units
             unitSI = "aSI", # Optional name of the SI units corresponding to the variable
             normSI = 10.0, # normalisation constant - i.e. aUnits/aSI
             )

print(a.name+":")
print(a.dims) # Variable dimensions
print(a.units)
print(a.data)

a:
['x']
aUnits
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15.]

If we specify normSI we can convert between the units

[4]:

a.switchUnits()
print(a.units)
print(a.data)

print("And back:")
a.switchUnits()
print(a.units)
print(a.data)

aSI
[  0.  10.  20.  30.  40.  50.  60.  70.  80.  90. 100. 110. 120. 130.
 140. 150.]
And back:
aUnits
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15.]

We can interrogate the variable about its properties directly.

[5]:

print(a.isFluid) # The variable lives only on the x grid
print(a.isDerived) # The variable is derived as opposed to implicit
print(a.isScalar) # The variable is a scalar - 0D (always derived)
print(a.isDistribution) # The variable lives on x,h,v
print(a.isSingleHarmonic) # The variable lives on x,v

True
False
False
False
False

It is possible to automatically generate a variable and its dual at construction

[6]:

a,a_dual = varAndDual("a",grid,
             data = np.linspace(0,15,16),
             units = "aUnits",
             unitSI = "aSI",
             normSI = 10.0,
             ); # The warning here is just to make sure we are doing things on purpose

[7]:

print(a.name)
print(a_dual.name)

print(a_dual.isDerived) # The dual variable is derived from the regular grid variable
print(a_dual.derivation.name) # The derivation name (more on this in next tutorial)

print(a.dual.name) # Each variable is aware of its dual
print(a_dual.dual.name)

a
a_dual
True
gridToDual
a_dual
a

We can also treat the variable on the dual grid as the primary variable

[8]:

b_dual,b = varAndDual("b",grid,primaryOnDualGrid=True) # Note the ordering:
                                                       # the dual variable is returned first

print(b_dual.isDerived)
print(b.isDerived)

False
True

Derivations will be covered in the next tutorial, but a particular class of derived variable is easy enough to construct without the knowledge of derivations to warrant examining here. These are derived variables based on node derivations, i.e. converting Python expressions into Fortran-parsable form.

Any variable can be wrapped into a Node object, which can then be used in Python expressions, producing new Nodes. Finally, a Node can be converted into a derived variable.

[9]:

cNode = node(a) + 2*node(b)**2 # This is a new node representing a + 2*b**2

c = varFromNode("c",grid,cNode) # Converting a node into a variable

print(type(c.derivation)) # Node derivation
print(c.derivationArgs) # c depends on a and b

print(c.latex()) # More on LaTeX representation in a later tutorial

<class 'RMK_support.derivations.NodeDerivation'>
['a', 'b']
\text{c}= \text{a}+2\left(\text{b}\right)^{2}

We can also manually create dual variable from an existing variable by using makeDual.

This will also associate the new variable as the dual of the old.

[10]:

c_dual = c.makeDual("c_dual")

print(c.dual.name)

c_dual

Or, if we do not wish to assign the dual variable to some object, we can associate a dual with withDual

[11]:

a2 = Variable("a2",grid).withDual()

print(a2.dual.name)

a2_dual

A note about variables on the dual grid in ReMKiT1D#

These variable live on the right cell faces, and only inner cell faces usually house evolved quantities in ReMKiT1D (except for periodic grids). However, since both variables on regular and the dual grid are stored in the same objects, the lengths of the data arrays are the same. This leads to the right-most cell face entry for dual grid variables to generally be unused, and set to 0 to avoid unintended behaviour. For the treatment of boundary conditions the reader is directed to the examples.

Variable Container#

While useful on their own, variables need to be stored in various containers in order to be used fully by both ReMKiT1D and the Python interface. The main container is the VariableContainer, which we construct and explore here

[12]:

vc = VariableContainer(grid) # An empty variable container

# We can add multiple variables to the container (note that by default the VariableContainer adds the duals of any added variables)

vc.add(a,b)

# Variable containers provide accessors that require knowing the variable name

print(vc["a"].units) # This will correctly index variable a and retrieve its units

aUnits

[13]:

# Variables can also be renamed when adding them

vc.add(c.rename("c2"))

c2 = vc["c2"] # This will retrieve the variable that we just added, which is c, but renamed
print(c2.name) # Check that it is indeed renamed
print(c2.derivationArgs)

# We can achieve the same effect of renaming and addition to the container like this

vc["c3"] = c
print(vc["c3"].name) # Check that it is indeed renamed
print(vc["c3"].derivationArgs)

# Note that this does not change the original variable

print(c.name)

c2
['a', 'b']
c3
['a', 'b']
c

The coordinates of the container include both the regular and staggered grid, and can also include time (mostly used when loading multiple data files).

[14]:

print(vc.coords)

{'x': array([0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05,
       1.15, 1.25, 1.35, 1.45, 1.55]), 'x_dual': array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2, 1.3,
       1.4, 1.5, 1.6]), 'h': array([0, 1, 2, 3]), 'v': array([0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75])}

Both individual variables and containers can be converted into xarray DataArray/Dataset objects

[15]:

print(a.dataArr)

<xarray.DataArray (x: 16)>
array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,
       13., 14., 15.])
Dimensions without coordinates: x
Attributes:
    isDerived:         False
    isDistribution:    False
    units:             aUnits
    isStationary:      False
    isScalar:          False
    isOnDualGrid:      False
    priority:          0
    derivationRule:    none
    isSingleHarmonic:  False
    normSI:            10.0
    unitSI:            aSI

[16]:

print(vc.dataset)

<xarray.Dataset>
Dimensions:  (x: 16, x_dual: 16, h: 4, v: 8, dim_0: 1)
Coordinates:
  * x        (x) float64 0.05 0.15 0.25 0.35 0.45 ... 1.15 1.25 1.35 1.45 1.55
  * x_dual   (x_dual) float64 0.1 0.2 0.3 0.4 0.5 0.6 ... 1.2 1.3 1.4 1.5 1.6
  * h        (h) int64 0 1 2 3
  * v        (v) float64 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75
Dimensions without coordinates: dim_0
Data variables:
    time     (dim_0) float64 0.0
    a        (x) float64 0.0 1.0 2.0 3.0 4.0 5.0 ... 11.0 12.0 13.0 14.0 15.0
    a_dual   (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    b        (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    b_dual   (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    c2       (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    c2_dual  (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    c3       (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    c3_dual  (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0

Note that a_dual was initialised with the same values as a. For some derived variables we can evaluate them using the dataset. This depends on whether the derivation has a Python evaluation method defined. Interpolation does:

[17]:

a_dual.evaluate(vc.dataset)

#A following call to vc.dataset will now have the interpolated value of a_dual - showing that adding variables does not copy them!
print(vc.dataset)

<xarray.Dataset>
Dimensions:  (x: 16, x_dual: 16, h: 4, v: 8, dim_0: 1)
Coordinates:
  * x        (x) float64 0.05 0.15 0.25 0.35 0.45 ... 1.15 1.25 1.35 1.45 1.55
  * x_dual   (x_dual) float64 0.1 0.2 0.3 0.4 0.5 0.6 ... 1.2 1.3 1.4 1.5 1.6
  * h        (h) int64 0 1 2 3
  * v        (v) float64 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75
Dimensions without coordinates: dim_0
Data variables:
    time     (dim_0) float64 0.0
    a        (x) float64 0.0 1.0 2.0 3.0 4.0 5.0 ... 11.0 12.0 13.0 14.0 15.0
    a_dual   (x_dual) float64 0.5 1.5 2.5 3.5 4.5 ... 11.5 12.5 13.5 14.5 15.5
    b        (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    b_dual   (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    c2       (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    c2_dual  (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    c3       (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    c3_dual  (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0

The time variable is automatically added to all variable containers, and can become its own coordinate (see examples when loading simulation results)

[18]:

print(vc["time"].isScalar)
print(vc["time"].isDerived)

True
True

HDF5 IO#

We can write the variable container defined above to a ReMKiT1D-readable HDF5 file - useful when producing custom initial conditions.

[19]:

writeRMKHDF5(vc,"test.h5")

We can then load either individual variables from the h5 file, or the entire variable container if we have it

[20]:

loaded_a = loadVariableFromHDF5(a,["test.h5"]) # We need a list of files to load - allows for time series

print(loaded_a.dataArr)

<xarray.DataArray (x: 16)>
array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,
       13., 14., 15.])
Dimensions without coordinates: x
Attributes:
    isDerived:         False
    isDistribution:    False
    units:             aUnits
    isStationary:      False
    isScalar:          False
    isOnDualGrid:      False
    priority:          0
    derivationRule:    none
    isSingleHarmonic:  False
    normSI:            10.0
    unitSI:            aSI

[21]:

loaded_vc = loadVarContFromHDF5(*tuple(vc.variables), # *args should be a tuple of variables
                                                      # this way we load all variables in the container
                                filepaths=["test.h5"]
                                )
print(loaded_vc.dataset)

<xarray.Dataset>
Dimensions:  (x: 16, x_dual: 16, h: 4, v: 8, t: 1, dim_0: 1)
Coordinates:
  * x        (x) float64 0.05 0.15 0.25 0.35 0.45 ... 1.15 1.25 1.35 1.45 1.55
  * x_dual   (x_dual) float64 0.1 0.2 0.3 0.4 0.5 0.6 ... 1.2 1.3 1.4 1.5 1.6
  * h        (h) int64 0 1 2 3
  * v        (v) float64 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75
  * t        (t) float64 0.0
Dimensions without coordinates: dim_0
Data variables:
    time     (dim_0) float64 0.0
    a        (x) float64 0.0 1.0 2.0 3.0 4.0 5.0 ... 11.0 12.0 13.0 14.0 15.0
    a_dual   (x_dual) float64 0.5 1.5 2.5 3.5 4.5 ... 11.5 12.5 13.5 14.5 15.5
    b        (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    b_dual   (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    c2       (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    c2_dual  (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    c3       (x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    c3_dual  (x_dual) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0

/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:807: UserWarning: Variable time already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:807: UserWarning: Variable a_dual already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:818: UserWarning: Variable a already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:807: UserWarning: Variable b_dual already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:818: UserWarning: Variable b already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:807: UserWarning: Variable c2_dual already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:818: UserWarning: Variable c2 already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:807: UserWarning: Variable c3_dual already in VariableContainer. Overwriting.
  warnings.warn(
/home/stefan/.local/lib/python3.8/site-packages/RMK_support/variable_container.py:818: UserWarning: Variable c3 already in VariableContainer. Overwriting.
  warnings.warn(

We can also load only using variable names, and the code will try its best to infer the properties of the variables (will lose a lot of information, but should keep dimensionality info).

[22]:

loaded_vc_2 = loadFromHDF5(grid,["a","a_dual"],filepaths=["test.h5"])

print(loaded_vc_2.dataset)

<xarray.Dataset>
Dimensions:  (x: 16, x_dual: 16, h: 4, v: 8, t: 1, dim_0: 1)
Coordinates:
  * x        (x) float64 0.05 0.15 0.25 0.35 0.45 ... 1.15 1.25 1.35 1.45 1.55
  * x_dual   (x_dual) float64 0.1 0.2 0.3 0.4 0.5 0.6 ... 1.2 1.3 1.4 1.5 1.6
  * h        (h) int64 0 1 2 3
  * v        (v) float64 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75
  * t        (t) float64 0.0
Dimensions without coordinates: dim_0
Data variables:
    time     (dim_0) float64 0.0
    a        (x) float64 0.0 1.0 2.0 3.0 4.0 5.0 ... 11.0 12.0 13.0 14.0 15.0
    a_dual   (x) float64 0.5 1.5 2.5 3.5 4.5 5.5 ... 11.5 12.5 13.5 14.5 15.5

[23]:

print(loaded_vc_2["a_dual"].isDerived) # Does not know that a_dual is derived

print(loaded_vc["a_dual"].isDerived) # Correct because it was loaded using the Variable objects

False
True