 # Validation of a PETSc based software implementing a 4DVAR Data Assimilation algorithm: a case study related with an Oceanic Model based on Shallow Water equation

In this work are presented and discussed some results related to the validation process of a software module based on PETSc which implements a Data Assimilation algorithm.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Definition of the DA problem and description of the PETSc based software implementation

The considered software intends to solve the problem defined on a suitable domain decomposition

 DD(Δ×Ω)={Δj×Ωh}j;h

of time-space domain as described in Definition 1 (all the needed notations can be founded in ).

###### Definition 1 (The 4D-VAR DA problem defined on the domain decomposition DD(Δ×Ω) - the 4D-VAR DD-DA problem).

The 4D Variational DD-DA problem consists in computing the vector

such that

 ~xDA = ∑j;hEOjh(xDAjh) (1)

where

 xDAjh = argminxjhJjh(ujh), (2)

where the operator (the local 4D-VAR regularization functional) is defined as follows:

 Jjh(xjh) = ROjh[J]+μ O(ih)(jk)[x(ih)(jk)] (3)

and where is a suitably defined operator on the overlapped domain . Parameter is a regularization parameter. The 4D-VAR regularization functional is defined as:

 J(x)=∥x−xb∥2B−1+λntobs−1∑k=0∥Htk(Mt0→tk[x])−vk∥2R−1k (4)

where is a regularization parameter, and () are the covariance matrices of the errors on the background and the observations respectively, while and denote the weighted euclidean norm.

The 4D-VAR DD-DA problem solution is computed performing the following steps on each subdomains (the so called 4D-VAR DD-DA algorithm):

• Locally compute all the parameters that define the local 4D-VAR regularization functional

• Locally compute the minimum (needed values for overlapping regions are obtained when necessary - i.e., for the model local evolution)

• Globally contribute to computation of

In order to compute the minimum of all the functionals , the DD-DA algorithm has to face with some issues. In more details, we have to address:

• the linearization of the operator , let us say , used for the evaluations of required by the minimisation algorithm;

• the evaluation of the adjoint operator of , let us say , used for the evaluation of required by the minimisation algorithm;

Both the points above should require the computation of the discretization of the Jacobian of

 ∇Mt−Δt→t

Following some details about software implementation in PETSc (Portable, Extensible Toolkit for Scientific Computation) environment. To implement the entire algorithm we plan to use:

1. the PETSc time steppers TS module for solving time-dependent (nonlinear) PDEs, including the computation of adjoint;

2. the PETSc DM module wich is a powerfull tool for the managment of all mesh data related with domain decomposition;

3. The TAO software library  for the computation of (2). The Toolkit for Advanced Optimization (TAO) is aimed at the solution of large-scale optimization problems on high-performance architectures. TAO is suitable for both single-processor and massively-parallel architectures. The current version of TAO has algorithms for unconstrained and bound-constrained optimization.

4. The SLEPc software library  for the computation of spectral decomposition usefull to compute a preconditioner of the error covariance matrices (i.e., see approach used in 

). The Scalable Library for Eigenvalue Problem Computations (SLEPc) is a software library for the solution of large scale sparse eigenvalue problems on parallel computers. It can also be used for computing a partial SVD of a large, sparse, rectangular matrix, and to solve nonlinear eigenvalue problems.

All the abobe mentioned software are integrated or based on PETSc (see figure 1 for a representation of the Software stack and algorithm implementation).

To represent the Jacobian of in PETSc we decided to follow a “matrix-free approach”: we used a MATHSHELL type for PETSc Mat object to represent just defining its way of operating

## 2 Case study description

The case study is based on the Shallow Water Equations (SWEs) on the sphere. The SWE have been used extensively as a simple model of the atmosphere or ocean circulation since they contain the essential wave propagation mechanisms found in general circulation models . The SWEs in spherical coordinates are:

 ∂u∂t = −1acosθ(u∂u∂λ+vcosθ∂u∂θ)+(f+utanθa)v−gacosθ∂h∂λ (5) ∂v∂t = −1acosθ(u∂v∂λ+vcosθ∂v∂θ)+(f+utanθa)u−ga∂h∂θ (6) ∂h∂t = −1acosθ(∂(hu)∂λ+∂(hucosθ)∂θ) (7)

Here is the Coriolis parameter given by , where is the angular speed of the rotation of the Earth, is the height of the homogeneous atmosphere (or of the free ocean surface), and are the zonal and meridional wind (or the ocean velocity) components, respectively, and are the latitudinal and longitudinal directions, respectively, is the radius of the earth and is the gravitational constant.

We express the system of equations (5)-(7) using a compact form, i.e.:

 ∂Z∂t=Mt−Δt→t(Z) (8)

where

 Z=⎛⎜⎝uvh⎞⎟⎠ (9)

and

 Mt−Δt→t(Z) = ⎛⎜ ⎜ ⎜ ⎜ ⎜⎝−1acosθ(u∂u∂λ+vcosθ∂u∂θ)+(f+utanθa)v−gacosθ∂h∂λ−1acosθ(u∂v∂λ+vcosθ∂v∂θ)+(f+utanθa)u−ga∂h∂θ−1acosθ(∂(hu)∂λ+∂(hucosθ)∂θ)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (10) = ⎛⎜⎝F1F2F3⎞⎟⎠

We discretize (8) just in space using an un-staggered Turkel-Zwas scheme [5, 6], and we obtain:

 ∂Zdisc∂t=Mt−Δt→tdisc(Zdisc) (11)

where

 Zdisc=⎛⎜ ⎜ ⎜⎝(ui,j)i=0,…,nlon−1;j=0,…,nlat−1(vi,j)i=0,…,nlon−1;j=0,…,nlat−1(hi,j)i=0,…,nlon−1;j=0,…,nlat−1⎞⎟ ⎟ ⎟⎠ (12)

and

 Mt−Δt→tdisc(Zdisc)=⎛⎜ ⎜ ⎜⎝(Ui,j)i=0,…,nlon−1;j=0,…,nlat−1(Vi,j)i=0,…,nlon−1;j=0,…,nlat−1(Hi,j)i=0,…,nlon−1;j=0,…,nlat−1⎞⎟ ⎟ ⎟⎠ (13)
 Ui,j = −σlonui,jcosθj(ui+1,j−ui−1,j) −σlat vi,j(ui,j+1−ui,j−1) −σlongpcosθj(hi+p,j−hi−p,j) +2[(1−α)(2Ωsinθj+ui,jatanθj)vi,j +α2(2Ωsinθj+ui+p,jatanθj)vi+p,j +α2(2Ωsinθj+ui−p,jatanθj)vi−p,j] Vi,j = −σlonui,jcosθj(vi+1,j−vi−1,j) −σlat vi,j(ui,j+1−ui,j−1) −σlatgq(hi,j+q−hi,j−q) −2[(1−α)(2Ωsinθj+ui,jatanθj)ui,j +α2(2Ωsinθj+q+ui,j+qatanθj+q)ui,j+q +α2(2Ωsinθj−q+ui,j−qatanθj−q)ui,j−q] Hi,j = −α{ui,jcosθj(hi+1,j−hi−1,j) + vi,j(hi,j+1−hi,j−1) +α2(ui+p,j+q−ui−p,j+q+ui+p,j−q−ui−p,j−q)]1p +[(1−α)(vi,j+qcosθj+q−vi,j−qcosθj−q) +α2(vi+p,j+qcosθj+q−vi+p,j−qcosθj−q) +α2(vi−p,j+qcosθj+q−vi−p,j−qcosθj−q)]1q}

If we define as follows:

 M0,MstepsΔt(⋅)=Mt−Δt→tdisc(Mt−Δt→tdisc(⋯Mt−Δt→tdisc(⋅)))Msteps, (14)

we note that the symbol represents the model “applied” times.

We also note that the numerical model is defined by the following parameters:

• discretization step in time domain,

• , discretization step in space domain,

• parameter of the Turkel-Zwas schema,

• , parameters of the Turkel-Zwas schema.

To verify the correct operation of the software module which implements the model, we tested the computed values of , when (i.e., when domain is not decomposed), where:

1. ,

2. ,

3. and ,

4. is a syntetic vector containing all the considered fields: the see-level field is generated by a Gaussian stochastic process; both velocity fields and are set to zero,

5. e defined on the basis of discretization grid used by data available at repository Ocean Synthesis/Reanalysis Directory of Hamburg University (see ).

We note that the values for parameters at above mentioned points 1, 2 and 3 were chosen on the basis of the considerations and results described in .

In figure 2-(a) is represented. In figure 2-(b)-(d) are represented respectively where . We used different values for with the aim to empirically determine the “best value” for . The considered values for were chosen taking into account the considerations about CFL condition for Turkel-Zwas methods described in .

To give a measure of how differs from depending from , in table 1 we show

 ∥M0,Msteps=30Δt(x0)−x0∥2/∥x0∥2

for the considered values of .

## 3 Test results related with the DA software module operation

To verify the correct operation of the software module which implements the DA process we tested the computed values of , when (i.e., when domain is not decomposed), starting:

• from the () diagonal matrix representing the covariance matrix of the errors on all the observations vectors

• from the diagonal matrix representing the observational operator and

• from the backgroud and

• from the set of observations vectors ,

• and by using, as preconditioner of covariance matrix , its Truncated SVD (the first singular values of are considered),

where

1. ,

2. Problem 1 For each ,

(each elements of were obtained by rounding, on the third significant digit, the respective elements of ).

is the identity matrix.

Problem 2 For each , if is a multiple of
else
( are sparse vectors whose lenght is and whose non-zero elements (the of the total) were obtained by adding a scaled number

from normal distribution to the respective elements of

).
Problem 3 For each ,
(each elements of were obtained by rounding, on the second significant digit, the respective elements of ). is the identity matrix.
Problem 4 For each , if is a multiple of
else
( are sparse vectors whose lenght is and whose non-zero elements (the of the total) were obtained by adding, to the respective elements of , a number from normal distribution scaled by a factor related with order of magnitude of each elements ).

3. ,

4. .

5. .

6. , where

 (xΔterr)i=(xΔtb)i−νΔt

and where

 νΔt=∑j=0,ldots,3∗nlat∗nlon−1(xΔtb)j3∗nlat∗nlon

Figure 3 shows how the background (the red circle in the image) and observations (the blue circle in the image) were chosen/built from data: in particular, for each values of , the image intends to show which subset of is considered to generate background and observations. In particular, the procedure used to build the input data for DA problem performs the following steps:

1. “application” of the model to the starting point to obtain the vector where ;

2. from the vector computed at above point 1 we obtain both the background vector and, by using one of the definitions for Problem 1 or Problem 2, the first observation vector ;

3. further “application” of the model to compute the set of vectors from which obtain, by using one of the definitions for Problem 1 or Problem 2, the remaining observation vectors .

Then, the “assimilation window” in time domain, when the value of is fixed, is the intervall defined as:

 AWntobsΔt=[(30−ntobs+1)Δt,30Δt]. Figure 3: How xΔtb and Δtxnobs were chosen/built from M0,Msteps=mΔt(x0) data

Three sets of tests are performed: the first set intends to evaluate how different values for and influence the behavior of DA software module; the second and third sets of tests intend to give elements to evaluate how Data Assimilation used during model application (i.e., ) improve the quality of the model computed values without DA (i.e., ) with respect to observations 111We note that the symbols and represent the model “applied” times to the background and to the solution of DA problem respectively.

Tests Set 1 In order to evaluate how different values for and influence the behavior of DA software module, in tables 3 (for Problem 1), 4 (for Problem 2), 5 (for Problem 3) and 6 (for Problem 4) the values of and are showed for the above listed values of and and for the four considered problems where

 errΔtb = ∥∥HxΔtb−Δtxn=1obs∥∥2/∥∥Δtxn=1obs∥∥2 (15) errΔtDA = ∥∥HxΔtDA−Δtxn=1obs∥∥2/∥∥Δtxn=1obs∥∥2 (16)

and where . We note that values of are not reported when the algorithm failed (i.e., when the Truncated SVD computation failed). The table 2 show, as examples, the values of when and .

We also note that:

1. the smaller values for is (i.e., where ), more the DA software module is able to effectively compute the solution of the DA problem;

2. the larger values for is, more accurate is the solution of the DA problem computed by DA software module (if it is successful).

This behavior could be explained considering that:

1. the smaller the value of is, better conditioned is the matrix , but also

2. the larger values for is, “closer” the matrices and are.