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.
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 Ştefnescu 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 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.
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:
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 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 where
The T-Z scheme for (1) is given as follows. Give a constant (), 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 and , 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 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 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 , by solving the eigenvalue problem, and retaining the set of right singular vectors of corresponding to the largest singular values; that is, ,.
Then, we can obtain the POD basis of and as , , and give the approximation of ,, and ,
Now, we can use Galerkin projection method to the discrete T-Z FD model by replacing , , and with their approximations , , and , 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 . Their definitions can be easily obtained by substituting ,, and 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. Let , be a nonlinear function. If ,,, is a linearly independent set, for , then, for , the DEIM approximation of order for in the space spanned by is given by
Here, the matrix is generated by doing POD on the nonlinear snapshots ( may be a function defined from , and is the value of evaluated at ), ,. Next, interpolation is used to determine the coefficient vector by selecting rows , of the overdetermined linear system (14) to form a -by- linear system where , . The DEIM approximation of becomes
Now the only unknowns that need to be determined are the indices , 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 is selected by searching the position of the largest value of the first POD basis . Then, the remaining interpolation indices , , is determined by selecting the largest magnitude of the residual vector . For the linear independence of the input basis , is a nonzero vector and the output indices 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].
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.
Let , be the POD basis matrix of rank for the snapshots of the nonlinear function (obtained from T-Z scheme). With the aid of DEIM, we select a set of DEIM indices corresponding to , and determine the matrix . The DEIM approximation of assumes the form so that the nonlinear term in POD model can be approximated as
Similarly, we obtain the DEIM approximations of the rest of the nonlinear functions in POD model, where , , , , , , , , and are the corresponding POD basis matrices of rank for the snapshots of the nonlinear functions , , , , , , , , and . 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 . Then the dimension of the full discretized model is 2592. The integration time window is 24 h and we use 109 time steps with .
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.
(a) Contour of geopotential
(b) Wind field calculated from the geopotential field by using the geostrophic approximation
The numerical solutions generated by T-Z FD scheme at 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.
(a) Contour of geopotential
(b) Wind field calculated from the geopotential field by using the geostrophic approximation
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 . 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 , , and geopotential. The curves of singular values for and 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 .
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 and as inputs.
| (a) First 30 DEIM points for nonlinear function |
| (b) First 30 DEIM points for nonlinear function |
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 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 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 Table 2, we know that the DEIM algorithm is much more efficient when combined with POD in solving the model reduction problems.
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.
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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This paper is supported by China Scholarship Council. Pengefi Zhao expresses his special thanks to Professor I.M. Navon from Department of Scientific Computing at Florida State University for his kind guidance and also to Dr. R. Stefanescu from Department of Computer Science at Virginia Tech for providing the code of POD on SWE in rectangular coordinates. Also, many thanks are to the NSFC (Grant nos. 11171131, 61133011, 61170092, 60973088, 61202308, 11026043, and 41304087) and the Basic Research Program of Jilin University (450060481098). This work is supported in part by the National Basic Research Program of China (973) under Grant 2009CB219301, the National Public Benefit Scientific Research Foundation of China under Grant 201011078, the China Scholarship Council and the National Innovation Research Project for Exploration and Development of Oil Shale under Grant nos. OSP-02 and OSR-02. The authors gratefully acknowledge the anonymous reviewers for their hard work and good patience.