Storing operators and data in modules - a tutorial

Author: PA Browne. Time-stamp: <2016-10-21 17:47:45 pbrowne>

Motivation

The operators used for a specific model and observation network may be complex. Specifically we are thinking of \(H, B, Q\) and \(R\).

It would therefore be useful to either load or compute these operators once, before the timestep process begins and the operators are called possibly a large number of times.

EMPIRE has the capacity to do just this, and in this tutorial we shall show how this can be achieved.

Location of template modules

There are a number of "template" modules that can assist the user in storing their operators and data. These are known to EMPIRE and are compiled before model_specific.f90 so that they may be accessed from within any subroutine contained in model_specific.f90.

They can be found in the directory src/user/ . They include

• src/user/Rdata.f90 - a place to store data for \(R\) and \(H\)
• src/user/Bdata.f90 - a place to store data for \(B\)
• src/user/Qdata.f90 - a place to store data for \(Q\)

Within each of these files is a module, called Rdata, Bdata and Qdata respectively. Each of these has its own "load*" subroutine. This is the place to put any setup of the operators.

Description of an example B matrix

For simplicity, let us assume that the background error covariance matrix, \(B\), is diagonal with different entries in each of its diagonal components.

\(B = \begin{bmatrix} b_{11}^2 & 0 & 0 & \cdots & 0 \\ 0 & b_{22}^2 & 0 & \cdots & 0 \\ 0 & 0 & b_{33}^2 & \cdots & 0 \\ \vdots &\vdots &\vdots &\ddots & \vdots \\ 0 &0 &0 & \cdots & b_{nn}^2 \end{bmatrix}\)

Note that it is very easy to get variances and standard deviations mixed up - with the notation above, it is clear to see that \(b_{ii}\) is a standard deviation and \(b_{ii}^2\) is a variance.

With this definition of \(B\) it is very easy to write down its (matrix) squareroot:

\(B^{1/2} = \begin{bmatrix} b_{11} & 0 & 0 & \cdots & 0 \\ 0 & b_{22} & 0 & \cdots & 0 \\ 0 & 0 & b_{33} & \cdots & 0 \\ \vdots &\vdots &\vdots &\ddots & \vdots \\ 0 &0 &0 & \cdots & b_{nn} \end{bmatrix}\)

As \(B^{1/2}\) is a very sparse matrix, we shall not store all the zeros. We shall implicitly realise that it all off-diagonal entries are zero and store only the diagonal values.

Description of the module

All the data for each operator that you wish to store and have access to elsewhere should come between the implicit none and the contains statements in the module.

So in this case, the following lines

module bdata
    implicit none
contains

would become, after declaring an array to store the diagonal entries of \(B^{1/2}\)

module bdata
    implicit none
    real(kind=kind(1.0d0)), dimension(:), allocatable :: diagbhalf
contains

At this stage, the dimensions of the state vector may not be set, hence diagBhalf is of unknown size. We do know, however, that we wish this to be a 1-dimensional array.

Description of the loading subroutine

Below is an example of a loadB subroutine which we will talk through subsequently.

subroutine loadb
    use sizes
    implicit none
    integer :: i
    allocate(diagbhalf(state_dim))
    open(2,file='Bhalfdata.txt',action='read',status='old')
    do i = 1,state_dim
        read(2,*) diagbhalf(i)
    end do
    close(2)
end subroutine loadB

So what is each line doing?

use sizes

We need to know how large the matrix is, hence we get this data from the sizes module.

implicit none

Simply good practice with Fortran - this means we must declare all the variables we use.

integer :: i

Declare an integer i to be used as a counter.

allocate(diagbhalf(state_dim))

We need to allocate the correct amount of space for diagBhalf. The size of the system is stored in the sizes module, in the variable state_dim. Hence this line sets aside the correct amount of memory to store the data for the matrix.

open(2,file='Bhalfdata.txt',action='read',status='old')

Here we are assuming that wherever EMPIRE is executed, we will have a file called Bhalfdata.txt. On line \(i\) of this file will be \(b_{ii}\).

do i = 1,state_dim

Now loop over each of the diagonal entries of the matrix

read(2,*) diagbhalf(i)

Read from the input file the value that is in the txt file and store it in the correct entry in the array diagBhalf

close(2)

Close the input file.

Calling the loadB subroutine

Now that the mechanism for reading in the data from a file has been written, we must make sure this is called. To do so, you should add a line in configure_model to do this.

subroutine configure_model
  use pf_control
  use timestep_data
  use sizes
  use bdata  !THIS IS AN ADDITIONAL LINE TO USE BDATA MODULE
  implicit none
  state_dim = 40
  obs_dim = 20
  call timestep_data_set_total(pf%time_bwn_obs*pf%time_obs)
  print*,'#################################'
  print*,'######### SANITY CHECK ##########'
  print*,'#################################'
  print*,'## STATE DIMENSION = ',state_dim
  print*,'##  OBS  DIMENSION = ',obs_dim
  print*,'## TOTAL TIMESTEPS = ',tsdata%total_timesteps
  print*,'#################################'
  call loadb !THIS IS AN ADDITIONAL LINE TO USE BDATA MODULE
end subroutine configure_model

It is important that the call to loadb happens after state_dim is set, as its value is used within loadb.

Using the data that has been read in

Firstly, let's discuss the algorithm we shall use to implement \(B^{1/2}x\), for some input vector \(x\).

Mathematically, let us write this as, given \(\mathbf{x}\), compute \(\mathbf{u} = B^{1/2}\mathbf{x}\).

\(\mathbf{u}_i = \sum_j B^{1/2}_{i,j} \mathbf{x}_j = \sum_j b_{ij} \mathbf{x}_j\).

Now as \(B^{1/2}\) is diagonal, this sum reduces so that

\(\mathbf{u}_i = b_{ii} \mathbf{x}_i\).

Hence we can do an "component wise" multiplication within the code.

We need to access Bdata within the subroutine bhalf . To do so, we put a use Bdata statement before the implicit none statement of subroutine bhalf .

The subroutine would therefore look something like the following:

subroutine bhalf(nrhs,x,bx)
    use bdata
    use sizes
    implicit none
    integer, parameter :: rk=kind(1.0d+0)
    integer, intent(in) :: nrhs
    real(kind=rk), dimension(state_dim,nrhs), intent(in) :: x
    real(kind=rk), dimension(state_dim,nrhs), intent(out) :: bx
    integer :: i
    do i = 1,nrhs
        bx(:,i) = diagbhalf*x(:,i)
    end do
end subroutine Bhalf

This way we have access to the array diagBhalf that we stored in Bdata, and we implement mulitplication by \(B^{1/2}\) by a component-wise multiplication of the arrays for each different vector that we are provided.

Similar things can be done for other operators.