# POD-DEIM Based Model Order Reduction for the Spherical Shallow Water Equations with Turkel-Zwas Finite Difference Discretization

Authors: Pengfei Zhao; Cai Liu; Xuan Feng
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We consider the shallow water equations (SWE) in spherical coordinates solved by Turkel-Zwas (T-Z) explicit large time-step scheme. To reduce the dimension of the SWE model, we use a well-known model order reduction method, a proper orthogonal decomposition (POD). As the computational complexity still depends on the number of variables of the full spherical SWE model, we use discrete empirical interpolation method (DEIM) proposed by Sorensen to reduce the computational complexity of the reduced-order model. DEIM is very helpful in evaluating quadratically nonlinear terms in the reduced-order model. The numerical results show that POD-DEIM is computationally very efficient for implementing model order reduction for spherical SWE.


---

## Body

## 1. Introduction

The shallow water equations (SWE) are probably the simplest model that captures the basic features of fluid flow motion in the geosciences. They describes the evolution of an incompressible and inviscid fluid in response to gravitational and rotational accelerations and their solutions represent a propagating Rossby wave along with East-West propagating gravity waves.During the last 30 years the SWE has become an important test bed for numerical methods [1–7]. Experts, especially, have to face greater challenges to the spherical SWE than the one in rectangular coordinates due to changing grid density in the spherical geometry in part due to presence of poles. Neta and Navon applied the Turkel-Zwas (T-Z) explicit large time-step scheme on a limited-area domain for the shallow water equations [8]. The most important advantage of the T-Z explicit large time-step scheme is addressing the issue of fast and slow time scales in SWE by treating the terms associated with the fast gravity-inertia waves on a coarser grid but to a higher accuracy than the terms associated with slow Rossby waves. Neta et al. also applied T-Z scheme for spherical SWE, where a staggered T-Z finite difference scheme has been used to improve the accuracy of numerical solutions [9]. Here, before performing model order reduction, the numerical solutions are generated by this staggered T-Z finite difference scheme.With the aim of obtaining more efficient description of the large scale complex system, model reduction is very helpful in reducing the computational cost and preserving numerical accuracy. Proper orthogonal decomposition (POD) is a powerful and graceful model reduction method of data analysis to obtain low-dimensional approximate descriptions of high-dimensional systems. The POD is also referred to principal component analysis (PCA), Karhunen-Loève (KL) decomposition, empirical orthogonal functions (EOF) analysis, and so forth, in many different fields [10–12]. In model reduction of SWE on rectangular coordinate, POD has been applied by Cao et al. [13], Vermeulen and Heemink [14], Daescu and Navon [15, 16], Altaf [17], and Ştefa
˘nescu and Navon [18].Discrete empirical interpolation method (DEIM) is a discrete version of the classic empirical interpolation method raised by Barrault et al. [19]. Generally, DEIM is applied to approximate a nonlinear function by combining projection with interpolation, which is very helpful in reducing the dimension of the ordinary differential equations (ODEs) generated by the POD-Galerkin process. The approach and its application was proposed by Chaturantabut and Sorensen in [20] and also widely applied in solving many problems [21–25].In this paper, we introduce the POD method to reduce the dimension of spherical SWE. The POD-Galerkin technique reduces the dimension in the sense that much fewer variables are present; however the computational complexity of the reduced SWE model is still depending on the number of variables of the full SWE. To make the model reduction more computationally efficient, we will combine POD and DEIM together to reduce the computational complexity and keep the accuracy close to the general POD results. DEIM is used due to to presence of quadratically nonlinear terms in spherical SWE model.The organization of this paper is as follows. In Section2, the spherical SWE is presented using matrix operators, and a T-Z finite difference algorithm is presented. In Section 3, the details of POD method for spherical SWE with T-Z FD scheme are introduced. In Section 4, we introduce the methodology of DEIM and the combination of POD and DEIM. Section 5 is the application of POD-DEIM approach on the spherical SWE and some numerical results.


## 2. T-Z Finite Difference Scheme on Spherical SWE

Navon and Yu proposed the T-Z explicit large time-step Fortran program for solving the spherical SWE model [26].The spherical SWE is described as follows:(1)
∂
u
∂
t
+
1
a
cos
⁡
θ
[
u
∂
u
∂
λ
+
v
cos
⁡
θ
∂
u
∂
θ
]
-
(
f
+
u
a
tan
θ
)
v
+
g
a
cos
⁡
θ
∂
h
∂
λ
=
0
,
∂
v
∂
t
+
1
a
cos
⁡
θ
[
u
∂
v
∂
λ
+
v
cos
⁡
θ
∂
v
∂
θ
]
+
(
f
+
u
a
tan
θ
)
u
+
g
a
∂
h
∂
θ
=
0
,
∂
h
∂
t
+
1
a
cos
⁡
θ
[
∂
∂
λ
(
h
u
)
+
∂
∂
θ
(
h
v
cos
⁡
θ
)
]
=
0
.Here,f is the Coriolis parameter given by
(2)
f
=
2
Ω
sin
θ
,
where Ω is the angular speed of the rotation of the earth, h is the height of the homogeneous atmosphere, u and v are the zonal and meridional wind components, respectively, θ and λ are the latitudinal and longitudinal directions, respectively, a is the radius of the earth, and g is the gravitation constant.The initial conditions are the same as those used by the previous work [9], where the height is defined as
(3)
h
(
λ
,
θ
,
0
)
=
1
g
(
Φ
-
+
2
Ω
a
u
0
sin
3
θ
cos
⁡
θ
sin
λ
)
,and the velocities are obtained geostrophically to yield(4)
u
(
λ
,
θ
,
0
)
=
-
3
u
0
sin
θ
cos
⁡
2
θ
sin
λ
+
u
0
sin
3
θ
sin
λ
,
v
(
λ
,
θ
,
0
)
=
u
0
sin
2
θ
cos
⁡
λ
,
where
(5)
g
=
9.8
m
/
sec
2
,
Φ
-
=
5.768
×
1
0
4
m
2
/
sec
2
,
u
0
=
20
m
/
sec
,
a
=
6.370
×
1
0
6
m
,
Ω
=
7.292
×
1
0
-
5
rad
/
sec
.The T-Z scheme for (1) is given as follows. Give a constant α (0
<
α
≤
1), then we have
(6)
u
k
,
j
l
+
1
=
u
k
,
j
l
-
1
-
σ
[
u
i
,
j
l
cos
⁡
θ
j
(
u
k
+
1
,
j
l
-
u
k
-
1
,
j
l
)
+
v
k
,
j
l
(
u
k
,
j
+
1
l
-
u
k
,
j
-
1
l
)
V
V
V
+
g
p
cos
⁡
θ
j
(
h
k
+
p
,
j
l
-
h
k
-
p
,
j
l
)
]
+
2
Δ
t
[
(
1
-
α
)
(
f
j
+
u
k
,
j
l
a
tan
θ
j
)
v
k
,
j
l
VVVVV
+
α
2
(
f
j
+
u
k
+
p
,
j
l
a
tan
θ
j
)
v
k
+
p
,
j
l
VVVVV
+
α
2
(
f
j
+
u
k
-
p
,
j
l
a
tan
θ
j
)
v
k
-
p
,
j
l
]
,
v
k
,
j
l
+
1

=
v
k
,
j
l
-
1

-

σ

[
u
k
,
j
l
cos
⁡
θ
j
(
v
k
+
1
,
j
l

-

v
k
-
1
,
j
l
)

+

v
k
,
j
l
(
v
k
,
j
+
1
l

-

v
k
,
j
-
1
l
)
V
V
V
V
V
V
V
+
g
q
(
h
k
+
j
+
q
l
-
h
k
,
j
-
q
l
)
u
k
,
j
l
cos
⁡
θ
j
]
-
2
Δ
t
[
(
1
-
α
)
(
f
j
+
u
k
,
j
l
a
tan
θ
j
)
u
k
,
j
l
VVVVV
+
α
2
(
f
j
+
q
+
u
k
,
j
+
q
l
a
tan
θ
j
+
q
)
u
k
,
j
+
q
l
VVVVV
+
α
2
(
f
j
-
q
+
u
k
,
j
-
q
l
a
tan
θ
j
-
q
)
u
k
,
j
-
q
l
]
,
h
k
,
j
l
+
1

=
h
k
,
j
l
-
1

-

σ

{
u
k
,
j
l
cos
⁡
θ
j
(
h
k
+
1
,
j
l

-

h
k
-
1
,
j
l
)

+

v
k
,
j
l
(
h
k
,
j
+
1
l

-

h
k
,
j
-
1
l
)
V
V
V
V
V
×
h
k
,
j
l
cos
⁡
θ
j
[
(
1
-
α
)
(
u
k
+
p
,
j
l
-
u
k
-
p
,
j
-
q
l
)
V
V
V
V
V
V
V
V
V
V
I
+
α
2
(
u
k
+
p
,
j
+
q
l

-

u
k
-
p
,
j
+
q
l

+

u
k
+
p
,
j
-
q
l
V
V
V
V
V
V
V
V
V
V
V
V
V
V
I
-
u
k
-
p
,
j
-
q
l
)
]
1
p
+
h
k
,
j
l
cos
⁡
θ
j
V
V
V
V
I
I
×
[
(
1
-
α
)
(
v
k
,
j
+
q
l
cos
⁡
θ
j
+
q

-

v
k
,
j
-
q
l
cos
⁡
θ
j
-
q
)
VVVVVVII
+
α
2
(
v
k
+
p
,
j
+
q
l
cos
⁡
θ
j
+
q
-
v
k
+
p
,
j
-
q
l
cos
⁡
θ
j
-
q
VVVVVVVVVII
×
v
k
-
p
,
j
+
q
cos
⁡
θ
j
+
q
-
v
k
-
p
,
j
-
q
v
k
+
p
,
j
+
q
l
)
V
V
V
V
V
V
I
I
×
cos
⁡
θ
j
-
q
v
k
,
j
+
q
l
]
1
q
u
k
,
j
l
c
o
s
θ
j
}
,
where
(7)
σ
=
Δ
t
a
Δ
λ
=
Δ
t
a
Δ
θ
.Now, we define the vector form of the numerical solutions as(8)
u
(
l
)
=
(
u
1,1
l
,
u
1,2
l
,
…
,
u
1
,
N
y
l
,
u
2,1
l
,
…
u
2
,
N
y
l
,
u
3,1
l
…
,
u
N
x
,
N
y
l
)
T
,
v
(
l
)
=
(
v
1,1
l
,
v
1,2
l
,
…
,
v
1
,
N
y
l
,
v
2,1
l
,
…
v
2
,
N
y
l
,
u
3,1
l
…
,
v
N
x
,
N
y
l
)
T
,
h
(
l
)
=
(
h
1,1
l
,
h
1,2
l
,
…
,
h
1
,
N
y
l
,
h
2,1
l
,
…
h
2
,
N
y
l
,
h
3,1
l
…
,
h
N
x
,
N
y
l
)
T
.Then, the T-Z scheme in operator form is(9)
u
(
l
+
1
)
=
u
(
l
-
1
)
-
N
11
(
u
(
l
)
)
-
N
12
(
u
(
l
)
,
v
(
l
)
)
+
N
13
(
u
(
l
)
,
v
(
l
)
)
+
L
11
(
v
(
l
)
)
-
L
12
(
h
(
l
)
)
,
v
(
l
+
1
)
=
v
(
l
-
1
)
-
N
21
(
v
(
l
)
)
-
N
22
(
u
(
l
)
,
v
(
l
)
)
-
N
23
(
u
(
l
)
)
-
L
21
(
u
(
l
)
)
-
L
22
(
h
(
l
)
)
,
h
(
l
+
1
)
=
h
(
l
-
1
)
-
N
31
(
u
(
l
)
,
h
(
l
)
)
-
N
32
(
v
(
l
)
,
h
(
l
)
)
-
N
33
(
u
(
l
)
,
h
(
l
)
)
-
N
34
(
v
(
l
)
,
h
(
l
)
)
.Here,N
i
,
j and L
i
,
j are the corresponding nonlinear terms and linear terms in the operator formulas, respectively. Their expressions are as follows:
(10)
N
11
(
u
)
=
σ
[
1
c
o
s
θ
,
1
c
o
s
θ
,
…
,
1
c
o
s
θ
]
T
︸
N
x
*
u
*
A
λ
1
u
N
12
(
u
,
v
)
=
σ
v
*
A
θ
1
u
,
N
13
(
u
,
v
)
=
2
Δ
t
[
1
t
a
n
θ
,
1
t
a
n
θ
,
…
,
1
t
a
n
θ
]
T
︸
N
x
*
D
p
,
α
(
u
*
v
)
,
L
11
(
v
)
=
σ
D
p
,
α
2
[
f
,
f
,
…
,
f
]
T
︸
N
x
*
v
,
L
12
(
h
)
=
g
p
[
1
c
o
s
θ
,
1
c
o
s
θ
,
…
,
1
c
o
s
θ
]
T
︸
N
x
*
A
λ
2
h
,
N
21
(
v
)
=
σ
v
*
A
θ
1
v
,
N
22
(
u
,
v
)
=
σ
[
1
c
o
s
θ
,
1
c
o
s
θ
,
…
,
1
c
o
s
θ
]
T
︸
N
x
*
u
*
A
λ
1
v
,
N
23
(
u
)
=
2
Δ
t
[
1
t
a
n
θ
,
1
t
a
n
θ
,
…
,
1
t
a
n
θ
]
T
︸
N
x
*
D
q
,
α
(
u
*
u
)
,
L
21
(
v
)
=
σ
D
q
,
α
2
[
f
,
f
,
…
,
f
]
T
︸
N
x
*
u
,
L
22
(
h
)
=
g
q
[
1
c
o
s
θ
,
1
c
o
s
θ
,
…
,
1
c
o
s
θ
]
T
︸
N
x
*
A
θ
2
h
,
N
31
(
u
,
h
)
=
σ
[
1
c
o
s
θ
,
1
c
o
s
θ
,
…
,
1
c
o
s
θ
]
T
︸
N
x
*
u
*
A
λ
1
h
,
N
32
(
v
,
h
)
=
σ
v
*
A
θ
1
u
,
N
33
(
u
,
h
)
=
σ
p
[
1
c
o
s
θ
,
1
c
o
s
θ
,
…
,
1
c
o
s
θ
]
T
︸
N
x
*
h
*
A
p
,
q
,
α
(
u
)
,
N
34
(
v
,
h
)
=
σ
q
[
1
c
o
s
θ
,
1
c
o
s
θ
,
…
,
1
c
o
s
θ
]
T
︸
N
x
*
h
*
A
p
,
q
,
α
(
c
o
s
θ
*
v
)
,
where 1
/
c
o
s
θ
=
[
1
/
cos
⁡
θ
1
,
1
/
cos
⁡
θ
2
,
…
,
1
/
cos
⁡
θ
N
y
], 1
/
t
a
n
θ
=
[
1
/
tan
θ
1
,
1
/
tan
θ
2
,
…
,
1
/
tan
θ
N
y
], A
λ
i,A
θ
i,i
=
1,2 are proper constant matrices for discrete first-order differential operators which take into account the initial and boundary conditions, D
p
,
α and D
q
,
α are the corresponding constant coefficient matrices for the weighted average of u
*
v and u
*
u, A
p
,
q
,
α is an approximate constant matrix operator with the description of first-order difference combined with a four-point difference scheme, and “*” is the componentwise multiplication symbol.


## 3. POD on Spherical SWE

The initial conditions are also obtained by multiplying the equations of the original initial conditions by theU
T, V
T, and H
T. Though the POD method is helpful to reduce the dimension of the model, the nonlinear terms are still retaining high computational complexities depending on the spatial dimension n
x
y of the original model from both evaluating the nonlinear functions and performing matrix multiplications to project on POD bases. However, by applying the DEIM method, we can remove the dependency and also obtain an expected computational efficiency. The nonlinear terms can be approximated by DEIM with the precomputing process, so the computational cost should be decreased significantly. Another reason why we use DEIM here is also that the process is very easy to be achieved; namely, only a few interpolation indices need to be selected in order to evaluate the nonlinear terms. We will introduce DEIM in the next section.To make the POD version of the spherical SWE more intelligible, we should build the POD decomposition of each variable and make the POD basis separate at first. Let us take the variableu for example. For v and h, we just need to follow the similar procedure to obtain the POD basis.Make the number of POD basisN much less than the spatial dimension n
x
y of u, that is, N
≪
n
x
y, and choose to construct the POD basis U
∈
R
n
x
y
×
n
k,k
∈
N
+ by solving the eigenvalue problem,
(11)
Y
T
Y
u
^
i
=
λ
i
u
^
i
,
i
=
1,2
,
…
,
N
T
,
and retaining the set of right singular vectors of Y corresponding to the k largest singular values; that is, U
=
{
u
i
}
i
=
1
k,u
i
=
(
1
/
λ
i
)
Y
u
^
i.Then, we can obtain the POD basis ofv and h as u, H
∈
R
n
x
y
×
k, and give the approximation of u,v, and h,
(12)
u
(
t
n
)
≈
U
u
~
(
t
n
)
,
v
(
t
n
)
≈
V
v
~
(
t
n
)
,
h
(
t
n
)
≈
H
h
~
(
t
n
)
,
u
~
(
t
n
)
,
v
~
(
t
n
)
,
h
~
(
t
n
)
∈
R
k
,
U
,
V
,
H
∈
R
n
x
y
×
k
.Now, we can use Galerkin projection method to the discrete T-Z FD model by replacingu, v, and h with their approximations U
u
~, V
v
~, and H
h
~, respectively, and then premultiplying the corresponding equations by U
T, V
T, and H
T. By applying this procedure, we transform the discrete version of original model to the POD reduced model as follows:
(13)
u
~
(
l
+
1
)
=
u
~
(
l
-
1
)
-
U
T
N
^
11
(
u
~
(
l
)
)
-
U
T
N
^
12
(
u
~
(
l
)
,
v
~
(
l
)
)
+
U
T
N
^
13
(
u
~
(
l
)
,
v
~
(
l
)
)
+
U
T
L
^
11
(
v
~
(
l
)
)
-
U
T
L
^
12
(
h
~
(
l
)
)
,
v
~
(
l
+
1
)
=
v
~
(
l
-
1
)
-
V
T
N
^
21
(
v
~
(
l
)
)
-
V
T
N
^
22
(
u
~
(
l
)
,
v
~
(
l
)
)
-
V
T
N
^
23
(
u
~
(
l
)
)
-
V
T
L
^
21
(
u
~
(
l
)
)
-
V
T
L
^
22
(
h
~
(
l
)
)
,
h
~
(
l
+
1
)
=
h
~
(
l
-
1
)
-
H
T
N
^
31
(
u
~
(
l
)
,
h
~
(
l
)
)
-
H
T
N
^
32
(
v
~
(
l
)
,
h
~
(
l
)
)
-
H
T
N
^
33
(
u
~
(
l
)
,
h
~
(
l
)
)
-
H
T
N
^
34
(
v
~
(
l
)
,
h
~
(
l
)
)
,
where N
^
11
,
N
^
12
,
N
^
13
,
N
^
21
,
N
^
22
,
N
^
23
,
N
^
31
,
N
^
32
,
N
^
33
,
N
^
34
,
L
^
11
,
L
^
12
,
L
^
21
,
L
^
22
∈
R
k. Their definitions can be easily obtained by substituting u
≈
U
u
~,v
≈
V
v
~, and h
≈
H
h
~ into the corresponding equations.


## 4. POD-DEIM on Spherical SWE

### 4.1. Discrete Empirical Interpolation Method

In this section, the application of DEIM to POD reduced spherical SWE model will be introduced. DEIM is a discrete version of the empirical interpolation method proposed by Barrault et al. [19]. The methodology and application were analyzed by Chaturantabut and Sorensen in [20]. Also, the convergence theorem and error estimate of POD-DEIM nonlinear model reduction is presented in [27, 28].The major function of DEIM is to provide a great efficient way to approximate nonlinear terms by replacing the high complex algebraic operation process in solving the reduced model. Also, it is incorporated into the reduced-basis techniques to provide a better reduced-basis treatment of nonaffine and nonlinear parameterized PDEs in terms of CPU time.The realization process of DEIM is as follows. Letf
:
D
→
R
n,D
⊂
R
n be a nonlinear function. If U
=
[
u
1
,
u
2
,
…
,
u
m
],u
i
∈
R
n,i
=
1,2
,
…
,
m, is a linearly independent set, for m
≤
n, then, for τ
∈
D, the DEIM approximation of order m for f
(
τ
) in the space spanned by {
u
l
}
l
=
1
m is given by
(14)
f
(
τ
)
≈
U
c
(
τ
)
,
U
∈
R
n
×
m
,
c
(
τ
)
∈
R
m
.Here, the matrixU is generated by doing POD on the nonlinear snapshots f
(
τ
t
i
)
,
τ
t
i
∈
D (τ may be a function defined from [
0
,
T
]
→
D, and τ
t
i is the value of τ evaluated at t
i), i
=
1,2
,
…
,
n
s,n
s
>
0. Next, interpolation is used to determine the coefficient vector c
(
τ
) by selecting m rows ρ
1
,
ρ
2
,
…
,
ρ
m
,
ρ
i
∈
N
+, of the overdetermined linear system (14) to form a m-by-m linear system
(15)
P
T
U
c
(
τ
)
=
P
T
f
(
τ
)
,
where P
=
[
e
ρ
1
,
e
ρ
2
,
…
,
e
ρ
m
]
∈
R
n
×
m, e
ρ
i
=
[
0
,
…
0
,
1
︸
ρ
i
,
0
,
…
,
0
]
∈
R
n. The DEIM approximation of f
∈
R
n becomes
(16)
f
(
τ
)
≈
U
(
P
T
U
)
-
1
P
T
f
(
τ
)
.Now the only unknowns that need to be determined are the indicesρ
1
,
ρ
2
,
…
,
ρ
m, which can be obtained in terms of the pseudo-algorithm (see Algorithm 1).Algorithm 1:DEIM pseudo-algorithm for interpolation indices [18, 27].
INPUT: {
u
l
}
l
=
1
m
⊂
R
n (linearly independent),
OUTPUT: ρ
→
=
[
ρ
1
,
ρ
2
,
…
,
ρ
m
]
∈
N
m,
Step  1. [
|
ψ
|
ρ
1
]
=
max
⁡
|
u
1
|, ψ
∈
R and ρ
1 is the component position of the
largest absolute value of u
1, with the smallest index taken in case of a tie,
Step  2. U
=
[
u
1
]
∈
R
n, P
=
[
e
ρ
1
]
∈
R
n, ρ
→
=
[
ρ
1
]
∈
N
+,
Step  3. For l
=
2
,
…
,
m do
(a) Solve (
P
T
U
)
c
=
P
T
u
l for c
∈
R
l
-
1, U, P
∈
R
n
×
(
l
-
1
),
(b) r
=
u
l
-
U
c, r
∈
R
n,
(c) [
|
ψ
|
ρ
l
]
=
max
⁡
{
|
r
|
},
(d) U
←
[
U
u
l
], P
←
[
P
e
ρ
l
], ρ
→
←
[
ρ
→
ρ
l
],
end For.The DEIM is used to select a set of indices from a linearly independent basis. In the first step of the algorithm, the first DEIM interpolation indexρ
1 is selected by searching the position of the largest value of the first POD basis |
u
1
|. Then, the remaining interpolation indices ρ
l, l
=
2,3
,
…
,
m, is determined by selecting the largest magnitude of the residual vector |
r
|. For the linear independence of the input basis {
u
l
}
l
=
1
m, r is a nonzero vector and the output indices {
ρ
i
}
i
=
1
m are not repeating in each iteration. Indeed, DEIM is an efficient method to make decisions on estimating the coordinate value of the n dimensional nonlinear functions in a m dimensional subspace. The error bound of the DEIM approximation is provided by Chaturantabut and Sorensen [28].


### 4.2. POD-DEIM Model of Spherical SWE

The DEIM approximation will be applied on POD model of spherical SWE so the computational complexity will reduce in proportion of the number of reduced variables in POD-DEIM model.LetU
N
^
11
∈
R
n
x
y
×
m, m
≤
n
x
y be the POD basis matrix of rank m for the snapshots of the nonlinear function N
^
11 (obtained from T-Z scheme). With the aid of DEIM, we select a set of m DEIM indices corresponding to U
N
^
11, and determine the matrix P
N
^
11
∈
R
n
x
y
×
m. The DEIM approximation of N
^
11 assumes the form
(17)
N
^
11
(
u
^
)
≈
U
N
^
11
(
P
N
^
11
T
U
N
^
11
)
-
1
P
N
^
11
T
N
^
11
,
so that the nonlinear term U
T
N
^
11 in POD model can be approximated as
(18)
U
T
N
^
11
(
u
^
)
≈
U
T
U
N
^
11
(
P
N
^
11
T
U
N
^
11
)
-
1
P
N
^
11
T
N
^
11
.Similarly, we obtain the DEIM approximations of the rest of the nonlinear functions in POD model,(19)
U
T
N
^
12
(
u
^
,
v
^
)
≈
U
T
U
N
^
12
(
P
N
^
12
T
U
N
^
12
)
-
1
P
N
^
12
T
N
^
12
,
U
T
N
^
13
(
u
^
,
v
^
)
≈
U
T
U
N
^
13
(
P
N
^
13
T
U
N
^
13
)
-
1
P
N
^
13
T
N
^
13
,
V
T
N
^
21
(
v
^
)
≈
V
T
V
N
^
21
(
P
N
^
21
T
V
N
^
21
)
-
1
P
N
^
21
T
N
^
21
,
V
T
N
^
22
(
u
^
,
v
^
)
≈
V
T
V
N
^
22
(
P
N
^
22
T
V
N
^
22
)
-
1
P
N
^
22
T
N
^
22
,
V
T
N
^
23
(
u
^
)
≈
V
T
V
N
^
23
(
P
N
^
23
T
V
N
^
23
)
-
1
P
N
^
23
T
N
^
23
,
H
T
N
^
31
(
u
^
,
h
^
)
≈
H
T
V
N
^
31
(
P
N
^
31
T
H
N
^
31
)
-
1
P
N
^
31
T
N
^
31
,
H
T
N
^
32
(
v
^
,
h
^
)
≈
H
T
V
N
^
32
(
P
N
^
32
T
H
N
^
32
)
-
1
P
N
^
32
T
N
^
32
,
H
T
N
^
33
(
u
^
,
h
^
)
≈
H
T
V
N
^
33
(
P
N
^
33
T
H
N
^
33
)
-
1
P
N
^
33
T
N
^
33
,
H
T
N
^
34
(
v
^
,
h
^
)
≈
H
T
V
N
^
34
(
P
N
^
34
T
H
N
^
34
)
-
1
P
N
^
34
T
N
^
34
,
where U
N
^
12, U
N
^
13, U
N
^
21, U
N
^
22, U
N
^
23, U
N
^
31, U
N
^
32, U
N
^
33, and U
N
^
34 are the corresponding POD basis matrices of rank m for the snapshots of the nonlinear functions N
^
12, N
^
13, N
^
21, N
^
22, N
^
23, N
^
31, N
^
32, N
^
33, and N
^
34. Replace the POD nonlinear terms by the DEIM estimation nonlinear terms; then we have the POD-DEIM model of spherical SWE.The initial conditions are obtained by following the similar procedure of that on POD model. The numerical examples will be shown in the next section.


## 4.1. Discrete Empirical Interpolation Method

In this section, the application of DEIM to POD reduced spherical SWE model will be introduced. DEIM is a discrete version of the empirical interpolation method proposed by Barrault et al. [19]. The methodology and application were analyzed by Chaturantabut and Sorensen in [20]. Also, the convergence theorem and error estimate of POD-DEIM nonlinear model reduction is presented in [27, 28].The major function of DEIM is to provide a great efficient way to approximate nonlinear terms by replacing the high complex algebraic operation process in solving the reduced model. Also, it is incorporated into the reduced-basis techniques to provide a better reduced-basis treatment of nonaffine and nonlinear parameterized PDEs in terms of CPU time.The realization process of DEIM is as follows. Letf
:
D
→
R
n,D
⊂
R
n be a nonlinear function. If U
=
[
u
1
,
u
2
,
…
,
u
m
],u
i
∈
R
n,i
=
1,2
,
…
,
m, is a linearly independent set, for m
≤
n, then, for τ
∈
D, the DEIM approximation of order m for f
(
τ
) in the space spanned by {
u
l
}
l
=
1
m is given by
(14)
f
(
τ
)
≈
U
c
(
τ
)
,
U
∈
R
n
×
m
,
c
(
τ
)
∈
R
m
.Here, the matrixU is generated by doing POD on the nonlinear snapshots f
(
τ
t
i
)
,
τ
t
i
∈
D (τ may be a function defined from [
0
,
T
]
→
D, and τ
t
i is the value of τ evaluated at t
i), i
=
1,2
,
…
,
n
s,n
s
>
0. Next, interpolation is used to determine the coefficient vector c
(
τ
) by selecting m rows ρ
1
,
ρ
2
,
…
,
ρ
m
,
ρ
i
∈
N
+, of the overdetermined linear system (14) to form a m-by-m linear system
(15)
P
T
U
c
(
τ
)
=
P
T
f
(
τ
)
,
where P
=
[
e
ρ
1
,
e
ρ
2
,
…
,
e
ρ
m
]
∈
R
n
×
m, e
ρ
i
=
[
0
,
…
0
,
1
︸
ρ
i
,
0
,
…
,
0
]
∈
R
n. The DEIM approximation of f
∈
R
n becomes
(16)
f
(
τ
)
≈
U
(
P
T
U
)
-
1
P
T
f
(
τ
)
.Now the only unknowns that need to be determined are the indicesρ
1
,
ρ
2
,
…
,
ρ
m, which can be obtained in terms of the pseudo-algorithm (see Algorithm 1).Algorithm 1:DEIM pseudo-algorithm for interpolation indices [18, 27].
INPUT: {
u
l
}
l
=
1
m
⊂
R
n (linearly independent),
OUTPUT: ρ
→
=
[
ρ
1
,
ρ
2
,
…
,
ρ
m
]
∈
N
m,
Step  1. [
|
ψ
|
ρ
1
]
=
max
⁡
|
u
1
|, ψ
∈
R and ρ
1 is the component position of the
largest absolute value of u
1, with the smallest index taken in case of a tie,
Step  2. U
=
[
u
1
]
∈
R
n, P
=
[
e
ρ
1
]
∈
R
n, ρ
→
=
[
ρ
1
]
∈
N
+,
Step  3. For l
=
2
,
…
,
m do
(a) Solve (
P
T
U
)
c
=
P
T
u
l for c
∈
R
l
-
1, U, P
∈
R
n
×
(
l
-
1
),
(b) r
=
u
l
-
U
c, r
∈
R
n,
(c) [
|
ψ
|
ρ
l
]
=
max
⁡
{
|
r
|
},
(d) U
←
[
U
u
l
], P
←
[
P
e
ρ
l
], ρ
→
←
[
ρ
→
ρ
l
],
end For.The DEIM is used to select a set of indices from a linearly independent basis. In the first step of the algorithm, the first DEIM interpolation indexρ
1 is selected by searching the position of the largest value of the first POD basis |
u
1
|. Then, the remaining interpolation indices ρ
l, l
=
2,3
,
…
,
m, is determined by selecting the largest magnitude of the residual vector |
r
|. For the linear independence of the input basis {
u
l
}
l
=
1
m, r is a nonzero vector and the output indices {
ρ
i
}
i
=
1
m are not repeating in each iteration. Indeed, DEIM is an efficient method to make decisions on estimating the coordinate value of the n dimensional nonlinear functions in a m dimensional subspace. The error bound of the DEIM approximation is provided by Chaturantabut and Sorensen [28].


## 4.2. POD-DEIM Model of Spherical SWE

The DEIM approximation will be applied on POD model of spherical SWE so the computational complexity will reduce in proportion of the number of reduced variables in POD-DEIM model.LetU
N
^
11
∈
R
n
x
y
×
m, m
≤
n
x
y be the POD basis matrix of rank m for the snapshots of the nonlinear function N
^
11 (obtained from T-Z scheme). With the aid of DEIM, we select a set of m DEIM indices corresponding to U
N
^
11, and determine the matrix P
N
^
11
∈
R
n
x
y
×
m. The DEIM approximation of N
^
11 assumes the form
(17)
N
^
11
(
u
^
)
≈
U
N
^
11
(
P
N
^
11
T
U
N
^
11
)
-
1
P
N
^
11
T
N
^
11
,
so that the nonlinear term U
T
N
^
11 in POD model can be approximated as
(18)
U
T
N
^
11
(
u
^
)
≈
U
T
U
N
^
11
(
P
N
^
11
T
U
N
^
11
)
-
1
P
N
^
11
T
N
^
11
.Similarly, we obtain the DEIM approximations of the rest of the nonlinear functions in POD model,(19)
U
T
N
^
12
(
u
^
,
v
^
)
≈
U
T
U
N
^
12
(
P
N
^
12
T
U
N
^
12
)
-
1
P
N
^
12
T
N
^
12
,
U
T
N
^
13
(
u
^
,
v
^
)
≈
U
T
U
N
^
13
(
P
N
^
13
T
U
N
^
13
)
-
1
P
N
^
13
T
N
^
13
,
V
T
N
^
21
(
v
^
)
≈
V
T
V
N
^
21
(
P
N
^
21
T
V
N
^
21
)
-
1
P
N
^
21
T
N
^
21
,
V
T
N
^
22
(
u
^
,
v
^
)
≈
V
T
V
N
^
22
(
P
N
^
22
T
V
N
^
22
)
-
1
P
N
^
22
T
N
^
22
,
V
T
N
^
23
(
u
^
)
≈
V
T
V
N
^
23
(
P
N
^
23
T
V
N
^
23
)
-
1
P
N
^
23
T
N
^
23
,
H
T
N
^
31
(
u
^
,
h
^
)
≈
H
T
V
N
^
31
(
P
N
^
31
T
H
N
^
31
)
-
1
P
N
^
31
T
N
^
31
,
H
T
N
^
32
(
v
^
,
h
^
)
≈
H
T
V
N
^
32
(
P
N
^
32
T
H
N
^
32
)
-
1
P
N
^
32
T
N
^
32
,
H
T
N
^
33
(
u
^
,
h
^
)
≈
H
T
V
N
^
33
(
P
N
^
33
T
H
N
^
33
)
-
1
P
N
^
33
T
N
^
33
,
H
T
N
^
34
(
v
^
,
h
^
)
≈
H
T
V
N
^
34
(
P
N
^
34
T
H
N
^
34
)
-
1
P
N
^
34
T
N
^
34
,
where U
N
^
12, U
N
^
13, U
N
^
21, U
N
^
22, U
N
^
23, U
N
^
31, U
N
^
32, U
N
^
33, and U
N
^
34 are the corresponding POD basis matrices of rank m for the snapshots of the nonlinear functions N
^
12, N
^
13, N
^
21, N
^
22, N
^
23, N
^
31, N
^
32, N
^
33, and N
^
34. Replace the POD nonlinear terms by the DEIM estimation nonlinear terms; then we have the POD-DEIM model of spherical SWE.The initial conditions are obtained by following the similar procedure of that on POD model. The numerical examples will be shown in the next section.


## 5. Numerical Examples

This section has two parts. The first part consists of the numerical examples for modeling the POD model, and the other part is the numerical examples for modeling the POD-DEIM model.The initial conditions are the same as those used by Neta et al. [9] which has been introduced in Section 1.The space domain has been discretized using a mesh of 72× 36 points, with Δ
λ
=
Δ
θ
=
5
o. Then the dimension of the full discretized model is 2592. The integration time window is 24 h and we use 109 time steps with Δ
t
=
800
s.The T-Z FD spherical SWE scheme proposed by Neta et al. in [9] is employed in order to obtain the numerical results of the spherical SWE model. The CFL condition is given in Navon and de Villiers [29].In the first step, we give the initial geopotential isolines and the geostrophic wind field in Figure1.Contour of geopotential and wind field att
=
0.
(a)
Contour of geopotential
(b)
Wind field calculated from the geopotential field by using the geostrophic approximationThe numerical solutions generated by T-Z FD scheme att
=
24 h are illustrated in Figure 2. These two figures are isolines of geopotential which are continuous closed curve. For the planes are the expanding spherical coordinate planes, the left part of the isoline with the value of 58000 is shown at the upper right of the plane.Contour of geopotential and wind field att
=
24 h.
(a)
Contour of geopotential
(b)
Wind field calculated from the geopotential field by using the geostrophic approximationThe POD basis functions are constructed using 109 snapshots obtained from the numerical solutions of the full T-Z FD model in the time interval[
0
,
24
h
]. The dimension of the POD bases for each variable is 28, capturing more than 99.9% of the system energy. Figure 3 describes the singular values of the snapshots for u, v, and geopotential. The curves of singular values for u and v have similar tendency of motion. However the curve for geopotential is above the other two curves for the values of geopotential are much larger than u and v.Figure 3
Singular values of snapshots solutions.In order to improve the efficiency of the POD approximation, you use the DEIM algorithm by selecting 60 DEIM points on each nonlinear term. Figure4 illustrates the distribution of the first 30 spatial points selected by the DEIM algorithm using the POD bases of nonlinear functions N
13 and N
23 as inputs.First 30 points selected by DEIM for the nonlinear functionsN
13 and N
23.
(a)
First 30 DEIM points for nonlinear functionN
13
(b)
First 30 DEIM points for nonlinear functionN
23Figure5 describes the scalogram of the local errors between POD, POD-DEIM, and T-Z FD spherical SWE solutions, where we used 60 DEIM points to estimate the nonlinear terms.Figure 5
Local errors between POD, POD-DEIM solutions and the TZ FD solutions of spherical SWE att
=
24 h (Δ
t
=
800
s). The number of DEIM points was taken 60.Figure6 illustrates the correlation coefficients for the SWE variables, where POD bases dimensions and DEIM points number are chosen as 28 and 60, respectively. Then, we will test the average relative errors.Figure 6
Comparison of the correlation coefficients with POD and POD-DEIM foru, v, and geopotential. The number of DEIM points was taken 60.We calculated the average relative errors in Euclidian norm for all three variables of spherical SWE modelw
=
u
,
v
,
h with the following norm:
(20)
erro
r
POD
=
1
N
T
∑
i
=
1
t
n
‍
∥
w
T
Z
(
;
,
i
)
-
w
POD
(
:
,
i
)
∥
2
∥
w
T
Z
(
:
,
i
)
∥
2
,
erro
r
DEIM
=
1
N
T
∑
i
=
1
t
n
‍
∥
w
T
Z
(
;
,
i
)
-
w
DEIM
(
:
,
i
)
∥
2
∥
w
T
Z
(
:
,
i
)
∥
2
.The comparison of the average relative errors in Euclidian norm by POD and POD-DEIM using the numerical solution of the T-Z spherical SWE model is shown in Table1.Table 1
Average relative errors for each of the model variables. The POD bases dimension is 28 capturing more than 99.9% of the system energy. 60 DEIM points are chosen.
POD solutions
DEIM solutions
Error
u
2.6050
e
-
003
4.2692
e
-
003
Error
v
3.3648
e
-
003
4.8253
e
-
003
Error
h
4.8412
e
-
005
8.2988
e
-
005The CPU time on obtaining the TZ FD solution of spherical SWE model is 49.283 s. Additionally, we applied Galerkin projection to TZ FD model, and its computational time is 17.252 s. Also, when we choose 60 DEIM points to update the POD reduced model, we only need 0.127 s to obtain the numerical results. Runge-Kutta-Fehlberg method (RKF45) was used in solving not only the POD reduced model, but also the POD-DEIM model. From the computational time and root mean square errors in Table2, we know that the DEIM algorithm is much more efficient when combined with POD in solving the model reduction problems.Table 2
CPU time gains and the root mean square errors for each of the model variables att
=
24 h. The POD bases dimensions were 28 capturing more than 99.9% of the system energy. 60 DEIM points were chosen.
TZ FD solutions
POD solutions
DEIM solutions
CPU time (s)
49.283
17.252
0.127
RMSE
u
—
2.166
e
-
004
3.210
-
004
RMSE
v
—
3.449
e
-
004
4.143
e
-
004
RMSE
h
—
1.663
e
-
003
3.079
e
-
003We also made some other numerical examples with different numbers of the spatial discretization points. Figure7 describes that the CPU time in executing the TZ FD algorithm on solving spherical SWE model is very sensitive to the number of spatial discretization points. Also, in case of executing the POD algorithm, we have similar results. Compared to the previous two algorithms, POD-DEIM is affected much more slightly by the numbers of spatial discretization points. The nonlinear term evaluation with DEIM is a precomputing process. It is extremely fast in executing the DEIM process. Then, we can easily evaluate the nonlinear terms values and solve the ODEs after applying Garlerkin projection. All the other parts are the same to POD progress. Hence, the core of calculating the approximate solutions with POD-DEIM is just solving ODEs without the expensive matrix operation for nonlinear terms. This is why we can get excellent speed-up from DEIM.Figure 7
CPU time for TZ FD, POD, and POD-DEIM.It seems that the more the number of DEIM points we selected, the lower the errors we will get by using the theoretical results for DEIM. However, it is a little difficult for us to strike an average between which number of DEIM points is good enough and what error results are low enough. Also, for the little impact on speed-up, we can just choose the number of DEIM points large enough. For example, we can choose the number of DEIM points more than half of the number of time differentiations. Here we can choose 60, which is much larger than a half of 108. There is no doubt that if we can choose the number much larger, we will get the error results a little better and nearly make no difference in consuming the CPU time.


## 6. Conclusion

In this paper, we have applied the POD and DEIM algorithm on T-Z FD spherical SWE model. The POD is used to reduce the dimensions of the original model and obtain a reduced model with lower dimensions from the full model with very high dimensions in space domain. With the aid of DEIM algorithm, the nonlinear terms in the reduced model have been estimated accurately, and also we can omit the vast majority of computational steps which helps to save the computational cost. The numerical results show that POD-DEIM combined algorithm is very efficient and accurate in doing model order reduction on spherical shallow water equations.



---

*Source: 292489-2014-07-15.xml*