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

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 SWEmodel, we use discrete empirical interpolation method (DEIM) proposed by Sorensen to reduce the computational complexity of the reduced-ordermodel. DEIM is very helpful in evaluating quadratically nonlinear terms in the reduced-ordermodel.The numerical results show that POD-DEIM is computationally very efficient for implementing model order reduction for spherical SWE.


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][2][3][4][5][6][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][11][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 S ¸tef ȃ 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][22][23][24][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 Section 2, 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.

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: 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,  and V are the zonal and meridional wind components, respectively,  and  are the latitudinal and longitudinal directions, respectively,  is the radius of the earth, and  is the gravitation constant.
The initial conditions are the same as those used by the previous work [9], where the height is defined as and the velocities are obtained geostrophically to yield  (, , 0) = − 3 0 sin cos 2  sin  +  0 sin 3  sin , where ( The T-Z scheme for (1) is given as follows.Give a constant  (0 <  ≤ 1), then we have where Now, we define the vector form of the numerical solutions as Then, the T-Z scheme in operator form is Here,  , and  , are the corresponding nonlinear terms and linear terms in the operator formulas, respectively.Their expressions are as follows: where are proper constant matrices for discrete first-order differential operators which take into account the initial and boundary conditions,  , and  , are the corresponding constant coefficient matrices for the weighted average of u * k and u * u,  ,, 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.

POD on Spherical SWE
The initial conditions are also obtained by multiplying the equations of the original initial conditions by the   ,   , and   .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   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 variable  for example.For V and ℎ, we just need to follow the similar procedure to obtain the POD basis.
Make the number of POD basis  much less than the spatial dimension   of , that is,  ≪   , and choose to construct the POD basis  ∈ R   ×  ,  ∈ N + by solving the eigenvalue problem, and retaining the set of right singular vectors of  corresponding to the  largest singular values; that is,  = {  }  =1 ,   = (1/√  )û  .Then, we can obtain the POD basis of V and ℎ as ,  ∈ R   × , and give the approximation of u, k, and h, Now, we can use Galerkin projection method to the discrete T-Z FD model by replacing u, k, and h with their approximations ũ, Ṽ, and  h, respectively, and then premultiplying the corresponding equations by   ,   , and   .By applying this procedure, we transform the discrete version of original model to the POD reduced model as follows: where N11 , N12 , N13 , N21 , N22 , N23 , N31 , N32 , N33 , N34 , L11 , L12 , L21 , L22 ∈ R  .Their definitions can be easily obtained by substituting u ≈ ũ, k ≈ k, and h ≈  h into the corresponding equations.

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 reducedbasis 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.Let  :  → R  ,  ⊂ R  be a nonlinear function.If  = [ 1 ,  2 , . . .,   ],   ∈ R  ,  = 1, 2, . . ., , is a linearly independent set, for  ≤ , then, for  ∈ , the DEIM approximation of order  for () in the space spanned by {  }  =1 is given by Here, the matrix  is generated by doing POD on the nonlinear snapshots (   ),    ∈  ( may be a function defined from [0, ] → , and    is the value of  evaluated at   ),  = 1, 2, . . .,   ,   > 0. Next, interpolation is used to determine the coefficient vector () by selecting  rows  1 ,  2 , . . .,   ,   ∈ N + , of the overdetermined linear system (14) to form a -by- linear system where Now the only unknowns that need to be determined are the indices  1 ,  2 , . . .,   , which can be obtained in terms of the pseudo-algorithm (see Algorithm 1).
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 | 1 |.Then, the remaining interpolation indices   ,  = 2, 3, . . ., , is Step 3.For  = 2, . . .,  do end For.
determined by selecting the largest magnitude of the residual vector ||.For the linear independence of the input basis {  }  =1 ,  is a nonzero vector and the output indices {  }  =1 are not repeating in each iteration.Indeed, DEIM is an efficient method to make decisions on estimating the coordinate value of the  dimensional nonlinear functions in a  dimensional subspace.The error bound of the DEIM approximation is provided by Chaturantabut and Sorensen [28].

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.
Let  N11 ∈R   × ,  ≤   be the POD basis matrix of rank  for the snapshots of the nonlinear function N11 (obtained from T-Z scheme).With the aid of DEIM, we select a set of  DEIM indices corresponding to  N11 , and determine the matrix  N11 ∈ R   × .The DEIM approximation of N11 assumes the form so that the nonlinear term   N11 in POD model can be approximated as Similarly, we obtain the DEIM approximations of the rest of the nonlinear functions in POD model, where  N12 ,  N13 ,  N21 ,  N22 ,  N23 ,  N31 ,  N32 ,  N33 , and  N34 are the corresponding POD basis matrices of rank  for the snapshots of the nonlinear functions N12 , N13 , N21 , N22 , N23 , N31 , N32 , N33 , and N34 .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.

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  .Then the dimension of the full discretized model is 2592.The integration time window is 24 h and we use 109 time steps with Δ = 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 Figure 1.
The numerical solutions generated by T-Z FD scheme at  = 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.
The 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 , V, and geopotential.The curves of singular values for  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  and V.
In order to improve the efficiency of the POD approximation, you use the DEIM algorithm by selecting 60 DEIM points on each nonlinear term.Figure 4 illustrates the distribution of the first 30 spatial points selected by the DEIM algorithm using the POD bases of nonlinear functions  13 and  23 as inputs.
Figure 5 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 6 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.
We calculated the average relative errors in Euclidian norm for all three variables of spherical SWE model  = , V, ℎ with the following norm: 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 Table 1.
The 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 We also made some other numerical examples with different numbers of the spatial discretization points.Figure 7 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.
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.

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.Shale under Grant nos.OSP-02 and OSR-02.The authors gratefully acknowledge the anonymous reviewers for their hard work and good patience.
First 30 DEIM points for nonlinear function  23

Figure 6 :Figure 7 :
Figure 6: Comparison of the correlation coefficients with POD and POD-DEIM for , V, and geopotential.The number of DEIM points was taken 60.
and  1 is the component position of the largest absolute value of  1 , with the smallest index taken in case of a tie, Step 2.  = [ 1

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.

Table 2 :
CPU time gains and the root mean square errors for each of the model variables at  = 24 h.The POD bases dimensions were 28 capturing more than 99.9% of the system energy.60 DEIM points were chosen.