Structure Preserving Model Order Reduction of Shallow Water Equations

TL;DR

This paper introduces two structure-preserving reduced-order modeling techniques for the 2D shallow water equations, ensuring long-term stability and conservation of invariants like energy and enstrophy.

Contribution

It presents novel ROM construction methods that preserve Hamiltonian and quadratic structures of the SWE, enhancing stability and accuracy over long simulations.

Findings

Both methods preserve invariants such as energy and enstrophy.

The ROMs demonstrate high accuracy and computational efficiency.

Long-term stability of the solutions is achieved.

Abstract

In this paper, we present two different approaches for constructing reduced-order models (ROMs) for the two-dimensional shallow water equation (SWE). The first one is based on the noncanonical Hamiltonian/Poisson form of the SWE. After integration in time by the fully implicit average vector field method, ROMs are constructed with proper orthogonal decomposition/discrete empirical interpolation method (POD/DEIM) that preserves the Hamiltonian structure. In the second approach, the SWE as a partial differential equation with quadratic nonlinearity is integrated in time by the linearly implicit Kahan's method and ROMs are constructed with the tensorial POD that preserves the linear-quadratic structure of the SWE. We show that in both approaches, the invariants of the SWE such as the energy, enstrophy, mass, and circulation are preserved over a long period of time, leading to stable solutions. We conclude by demonstrating the accuracy and the computational efficiency of the reduced solutions by a numerical test problem.

Tables

Table 1. Table 1: Time averaged relative L 2 superscript 𝐿 2 L^{2} -errors

		${\Vert 𝐮-\widehat{𝐮}\Vert }_{Rel}$	${\Vert 𝐯-\widehat{𝐯}\Vert }_{Rel}$	${\Vert 𝐡-\widehat{𝐡}\Vert }_{Rel}$
POD-AVF	30 POD modes	1.192e-01	1.192e-01	1.473e-02
POD-DEIM-AVF	30 POD & 90 DEIM modes	2.874e-01	2.874e-01	3.478e-02
TPOD-Kahan	30 POD modes	1.265e-01	1.265e-01	1.567e-02

Table 2. Table 2: Time averaged absolute errors for the conserved quantities

		Energy	Enstrophy	Vorticity
POD-AVF	30 POD modes	3.768e-05	7.871e-03	2.123e-05
POD-DEIM-AVF	30 POD& 90 DEIM modes	1.968e-03	5.137e-03	1.042e-03
TPOD-Kahan	30 POD modes	2.901e-05	3.108e-03	3.454e-05

Table 3. Table 3: CPU time (in seconds) and speed-up factors

POD-DEIM-AVF				TPOD-Kahan
		CPU time	speed-up			CPU time	speed-up
FOM		529.80		FOM		788.89
POD	basis computation	33.46		POD	basis computation	31.02
	online computation	291.69			online computation	289.01
	total	325.15	1.63		total	321.03	2.45
DEIM	basis computation	33.89		TPOD	tensor computation [5](MP)	20.78 (9.43)
	online computation	24.10			online computation	6.42
	total	57.99	9.13		total [5](MP)	27.20 (15.89)	29.00 (49.64)

Equations

$$\frac{\partial u}{\partial t}$$

$$\frac{\partial v}{\partial t}$$

$$\frac{\partial h}{\partial t}$$

$$\xi (a,y,t)=\xi (b,y,t),\xi (x,c,t)=\xi (x,d,t),\xi \in \{u,v,h\},$$

$$u(\mathbf{x},0)=u_0(\mathbf{x}),v(\mathbf{x},0)=v_0(\mathbf{x}),h(\mathbf{x},0)=h_0(\mathbf{x}),$$

$$\frac{\partial z}{\partial t}=\mathcal{J}(z)\frac{\delta \mathcal{H}}{\delta z}=\left(\begin{array}{ccc}0& q& -{\partial }_x\\ -q& 0& -{\partial }_y\\ -{\partial }_x& -{\partial }_y& 0\end{array}\right)\left(\begin{array}{c}uh\\ vh\\ \frac{1}{2}(\upsilon \cdot \upsilon )+gh\end{array}\right),$$

$$\mathcal{H}(z)=\frac{1}{2}\iint h(\upsilon \cdot \upsilon +gh)d\mathbf{x}.$$

$$\{\mathcal{A},\mathcal{B}\}=\iint \left(q\frac{\delta ((\mathcal{A},\mathcal{B})}{\delta (u,v)}-\frac{\delta \mathcal{A}}{\delta \upsilon }\cdot \nabla \frac{\delta \mathcal{B}}{\delta h}+\frac{\delta \mathcal{B}}{\delta \upsilon }\cdot \nabla \frac{\delta \mathcal{A}}{\delta h}\right)d\mathbf{x},$$

$$\frac{\delta (\mathcal{A},\mathcal{B})}{\delta (u,v)}=\frac{\delta \mathcal{A}}{\delta u}\frac{\delta \mathcal{B}}{\delta v}-\frac{\delta \mathcal{B}}{\delta u}\frac{\delta \mathcal{A}}{\delta v}.$$

$$\{\mathcal{A},\{\mathcal{B},\mathcal{D}\}\}+\{\mathcal{B},\{\mathcal{D},\mathcal{A}\}\}+\{\mathcal{A},\{\mathcal{B},\mathcal{D}\}\}=0,$$

$$\frac{d\mathcal{H}}{dt}=\{\mathcal{H},\mathcal{H}\}=0.$$

$$\mathcal{C}=\iint hG(q)d\mathbf{x},$$

$$\{\mathcal{A},\mathcal{C}\}=0,\forall \mathcal{A}(\mathbf{z})\text{or}{\mathcal{J}}^{ij}\frac{\partial \mathcal{C}}{\partial {\mathbf{z}}^j}=0.$$

$$\mathcal{Z}=\frac{1}{2}\int \int h{q}^{2}d\mathbf{x}=\frac{1}{2}\int \int \frac{1}{h}{\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}+f\right)}^{2}d\mathbf{x},\mathcal{V}=\iint hqd\mathbf{x},\mathcal{M}=\iint hd\mathbf{x}.$$

$$\frac{\partial u}{\partial t}$$

$$\frac{\partial v}{\partial t}$$

$$\frac{\partial h}{\partial t}$$

$$\mathbf{u}(t)$$

$$\mathbf{v}(t)$$

$$\mathbf{h}(t)$$

$${\tilde{D}}_s=\left(\begin{array}{ccccc}0& 1& & & -1\\ -1& 0& 1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 0& 1\\ 1& & & 1& 0\end{array}\right)\in {\mathbb{R}}^{s\times s}.$$

$$D_x=\frac{1}{2\Delta x}{\tilde{D}}_{N_x}\otimes I_{N_y}\in {\mathbb{R}}^{N\times N},D_y=\frac{1}{2\Delta y}{I}_{N_x}\otimes {\tilde{D}}_{N_y}\in {\mathbb{R}}^{N\times N},$$

$$H^k=H({\mathbf{z}}^k)=\frac{1}{2}\sum _{i=1}^N(({\mathbf{u}}_i^k{)}^2+({\mathbf{v}}_i^k{)}^2+g{\mathbf{h}}_i^k){\mathbf{h}}_i^k\Delta x\Delta y,$$

$$Z^k=Z({\mathbf{z}}^k)=\frac{1}{2}\sum _{i=1}^N\frac{((D_x{\mathbf{v}}^k{)}_i-(D_y{\mathbf{u}}^k{)}_i+{\mathbf{f}}_i^k{)}^2}{{\mathbf{h}}_i^k}\Delta x\Delta y,$$

$$V^k=V({\mathbf{z}}^k)=\sum _{i=1}^N((D_x{\mathbf{v}}^k{)}_i-(D_y{\mathbf{u}}^k{)}_i+{\mathbf{f}}_i^k)\Delta x\Delta y,M^k=M({\mathbf{z}}^k)=\sum _{i=1}^N{\mathbf{h}}_i^k\Delta x\Delta y.$$

$$\frac{d\mathbf{z}}{dt}$$

$${\mathbf{z}}^{k+1}$$

$${\int }_0^1\nabla H(\xi ({\mathbf{z}}^{k+1}-{\mathbf{z}}^k)+{\mathbf{z}}^k)d\xi$$

$$=\left(\begin{array}{c}\frac{1}{3}\left({\mathbf{u}}^{k+1}\circ {\mathbf{h}}^{k+1}+{\mathbf{u}}^k\circ {\mathbf{h}}^k\right)+\frac{1}{6}\left({\mathbf{u}}^{k+1}\circ {\mathbf{h}}^k+{\mathbf{u}}^k\circ {\mathbf{h}}^{k+1}\right)\\ \frac{1}{3}\left({\mathbf{v}}^{k+1}\circ {\mathbf{h}}^{k+1}+{\mathbf{v}}^k\circ {\mathbf{h}}^k\right)+\frac{1}{6}\left({\mathbf{v}}^{k+1}\circ {\mathbf{h}}^k+{\mathbf{v}}^k\circ {\mathbf{h}}^{k+1}\right)\\ \frac{1}{6}\left({\mathbf{u}}^{k+1}\circ {\mathbf{u}}^{k+1}+{\mathbf{u}}^{k+1}\circ {\mathbf{u}}^k+{\mathbf{u}}^k\circ {\mathbf{u}}^k\right)+\frac{1}{6}\left({\mathbf{v}}^{k+1}\circ {\mathbf{v}}^{k+1}+{\mathbf{v}}^{k+1}\circ {\mathbf{v}}^k+{\mathbf{v}}^k\circ {\mathbf{v}}^k\right)+\frac{g}{4}\left({\mathbf{h}}^{k+1}+{\mathbf{h}}^k\right)\end{array}\right).$$

$$\frac{d\mathbf{u}}{dt}$$

$$\frac{d\mathbf{v}}{dt}$$

$$\frac{d\mathbf{h}}{dt}$$

$$\frac{d\mathbf{z}}{dt}=F(\mathbf{z})=R_1(\mathbf{z})+R_2(\mathbf{z})+L(\mathbf{z}),$$

$$R_1(\mathbf{z})=\left(\begin{array}{c}-\mathbf{u}\circ (D_x\mathbf{u})\\ -\mathbf{u}\circ (D_x\mathbf{v})\\ -D_x(\mathbf{u}\circ \mathbf{h})\end{array}\right),R_2(\mathbf{z})=\left(\begin{array}{c}-\mathbf{v}\circ (D_y\mathbf{u})\\ -\mathbf{v}\circ (D_y\mathbf{v})\\ -D_y(\mathbf{v}\circ \mathbf{h})\end{array}\right),L(\mathbf{z})=\left(\begin{array}{c}-g{D}_x\mathbf{h}+\mathbf{f}\circ \mathbf{v}\\ -g{D}_y\mathbf{h}-\mathbf{f}\circ \mathbf{u}\\ 0\end{array}\right).$$

$$\frac{{\mathbf{z}}^{k+1}-{\mathbf{z}}^k}{\Delta t}=R_1^f({\mathbf{z}}^k,{\mathbf{z}}^{k+1})+R_2^f({\mathbf{z}}^k,{\mathbf{z}}^{k+1})+\frac{1}{2}L({\mathbf{z}}^k+{\mathbf{z}}^{k+1}),$$

Full text

Keywords: Finite difference methods, linearly implicit methods, preservation of invariants, proper orthogonal decomposition, discrete empirical interpolation, tensorial proper orthogonal decomposition.

MSC classification: 65M06; 65P10; 37J05; 76B15; 15A69

1 Introduction

The shallow water equations (SWEs) consist of a set of two-dimensional partial differential equations (PDEs) describing a thin inviscid fluid layer flowing over the topography in a frame rotating about an arbitrary axis. SWEs are widely used in modeling large-scale atmosphere/ocean dynamics and numerical weather prediction. Energy and enstrophy are the most important conserved quantities of the SWEs, whereas the energy cascades to large scales whilst enstrophy cascades to small scales [4, 44]. Therefore numerical schemes that preserve the energy and the enstrophy of the SWEs lead to stable solutions in the long term integration [1]. Many geophysical flows can be written in Hamiltonian form [38]. The non-canonical Hamiltonian/Poisson form of the SWE in the rotational frame with constant Coriolis force was introduced first in [41]. Later on the Nambu formulation of the SWE [43], SWE with complete Coriolis force [20, 44], and multi-layer SWE [44] are developed. The discrete energy conservation follows from antisymmetry of the discrete Poisson bracket. Other conserved quantities are potential enstrophy, mass, and vorticity.

In this paper we consider two different formulations of the SWEs; as non-canonical Hamiltonian/Poisson PDE and as a PDE with quadratic nonlinear terms in the f-plane. Both versions of the SWEs are discretized in space using finite-differences by preserving the skew-symmetry in the Poisson matrix. The system of the semi-discretized ordinary differential equation (ODE) is also in Poisson form and has quadratic nonlinear terms. For time discretization of the SWE in Poisson form, we consider fully implicit average vector field (AVF) method [13, 15] that is energy preserving. On the other hand, the ODE system with quadratic nonlinear terms obtained from the SWE in the f-plane is discretized in time by the linearly implicit Kahan’s method [27, 12]. Both time integrators are second-order convergent in time as the finite-difference approximation in space. The fully implicit AVF method requires iterative solvers like Newton’s method at each time step for solving the nonlinear systems arising from time discretization. However, the linearly implicit Kahan’s method requires to solve only one linear system of equations in each time step. Both methods preserve well the conserved quantities of the SWEs, like the energy, enstrophy, mass, and vorticity in long-time integration.

Numerical methods for PDEs like the SWEs are computationally expensive and require a large amount of computer memory and computing time in real-time simulations. During the last decades, reduced-order models (ROMs) have emerged as a powerful approach to reduce the cost of evaluating large systems of PDEs by constructing a low-dimensional linear subspace (reduced space), that approximately represents the solution to the system of PDEs with a significantly reduced computational cost. The solutions of the high fidelity full-order model (FOM), generated by space-time discretization of PDEs, are projected onto reduced space using the proper orthogonal decomposition (POD). The POD has been widely used and is a computationally efficient reduced-order modeling technique in large-scale numerical simulations of nonlinear PDEs. Applying POD Galerkin projection, dominant modes of the PDEs are extracted from the snapshots of the FOM solutions. The computation of the FOM solutions is performed in the ”offline” stage, whereas the reduced system from the low-dimensional subspace is solved in the ”online” stage. Various ROMs are constructed for the SWEs, in conservative form [34, 33], in the $\beta$-plane [16, 17], in the f-plane [22] with POD and in the $\beta$-plane [8, 9] with dynamic mode decomposition (DMD). However, in these articles, preservation of the energy and other conservative quantities of the SWEs in the reduced space are not discussed.

The primary challenge in producing the low-dimensional models of the high-dimensional discretized PDEs is the efficient evaluation of the nonlinearities on the POD basis. The computational cost is reduced by sampling the nonlinear terms and interpolating, known as hyper-reduction techniques [3, 14, 2, 47, 39, 10]. Our method of choice is the discrete empirical interpolation method (DEIM) [14], which is one of the most frequently used hyper-reduction methods. The number of sampling points used in hyper-reduction methods often scales with the dimension of the reduced-order model, and rigorous error and stability analysis are in general not available. For PDEs like SWEs in the f-plane with polynomial nonlinearities, ROMs do not require approximating the nonlinear terms through sampling. Reduced-order operators can be precomputed in the ”offline” stage. Projection of FOM onto the reduced space yields low-dimensional matrix operators that preserve the polynomial structure of the FOM. This is an advantage because the offline-online computation is separated in contrast to the hyper-reduction methods. Recently, for PDEs with polynomial nonlinearities, the computationally efficient reduced-order models are constructed by the use of some tools from tensor theory and by matricizations of tensors [5, 6, 7, 25, 31]. In this paper, we have constructed two different ROMs for the SWEs; in Hamiltonian form and in f-plane. Our main contribution is twofold:

• •

ROMs are constructed that preserve the reduced skew-symmetry in the Poisson matrix of the SWEs with the energy preserving time integrator AVF method as in [23, 29]. We show that the reduced discrete energy (Hamiltonian) and other conserved quantities like enstrophy, mass and vorticity are well-preserved in the long-term using POD and DEIM.

• •

Semi-discrete form of the SWE in the f-plane results in an ODE with quadratic nonlinearities, which is solved in time by linearly implicit Kahan’s method. The computational complexity of the tensor products is ${\mathcal{O}}(n^{2}N^{2})$, where $n$ and $N$ are dimensions of reduced and full order models. Applying POD in a tensorial framework (TPOD) by exploiting matricizations of tensors, the computational cost is reduced. Consequently, the POD ROM for the quadratic system recovers an efficient offline-online decomposition and does not require an extra step of hyper-reduction, like DEIM. Here we make use of the sparse matrix technique MULTIPROD [32] to accelerate the tensor calculations.

We compare both approaches on a numerical test example. A large number of POD modes (about 50) are required for accurate ROM solutions, which is characteristic for problems with wave-type solutions like the SWEs. Numerical results demonstrate that the conserved quantities of the SWEs, energy, enstrophy, mass, and vorticity are well preserved by the FOMs and ROMs in long-time integration. Compared with the POD/DEIM ROMs with AVF discretization, the TPOD ROMs with Kahan’s discretization is much faster at the same level of accuracy.

The paper organized as follows. In the next section, we briefly review the SWE in Hamiltonian form and in the f-plane. The finite-difference discretization in space and time integration by the AVF method and Kahan’s method are presented in Section 3. POD/DEIM reduced models are constructed for the SWE with AVF semi-discretized Hamiltonian form, and TPOD reduced models for the semi-discrete ODE system with quadratic nonlinearities are given in Section 4. In Section 5, we compare POD, POD/DEIM, and TPOD for two-dimensional rotating SWE concerning the accuracy of the reduced solutions, preservation of the energy, enstrophy, vorticity, and computational efficiency. The paper ends with some conclusions.

2 Rotational shallow water equation

The two-dimensional rotational SWE on a rectangular space domain $\Omega=[a,b]\times[c,d]\subset\mathbb{R}^{2}$ ($a,b,c,d\in\mathbb{R}$) and with the spatial element ${\mathbf{x}}=(x,y)^{T}\in\mathbb{R}^{2}$, is given as [41]

[TABLE]

where $u({\mathbf{x}},t)$ and $v({\mathbf{x}},t)$ are the components of the (particle) velocity in $x$ and $y$-directions, respectively, $h({\mathbf{x}},t)$ is the fluid depth, and $q({\mathbf{x}},t)=(v_{x}({\mathbf{x}},t)-u_{y}({\mathbf{x}},t)+f({\mathbf{x}},t))/h({\mathbf{x}},t)$ is the potential vorticity with the Coriolis force $f({\mathbf{x}},t)$. The subscripts $x$ and $y$ denote the partial derivatives with respect to $x$ and $y$ components, respectively. Moreover, for the gravity constant $g$, it is given that $\Phi({\mathbf{x}},t)=(1/2)u({\mathbf{x}},t)^{2}+(1/2)v({\mathbf{x}},t)^{2}+gh({\mathbf{x}},t)$. The rotational SWE (1) is considered on a time interval $[0,T]$ for a final time $T>0$, and with periodic boundary conditions

[TABLE]

and with initial conditions

[TABLE]

where $u_{0}({\mathbf{x}})$, $v_{0}({\mathbf{x}})$ and $h_{0}({\mathbf{x}})$ are given functions.

2.1 SWE in Hamiltonian form

The non-canonical Hamiltonian/Poisson form of the SWE (1) is given by

[TABLE]

where $z({\mathbf{x}},t)=(u({\mathbf{x}},t),v({\mathbf{x}},t),h({\mathbf{x}},t))^{T}$ is the solution vector and $\mathbf{\upsilon}({\mathbf{x}},t)=(u({\mathbf{x}},t),v({\mathbf{x}},t))^{T}$ is the velocity field. The (energy) functional $\mathcal{H}(z)$ denotes the Hamiltonian given by

[TABLE]

The skew-symmetric Poisson bracket is defined for any two functionals $\mathcal{A}$ and $\mathcal{B}$ [35, 41] as

[TABLE]

where $\nabla=(\partial x,\partial y)^{T}$, and $\delta\mathcal{A}/\delta{\mathbf{\upsilon}}$ is the functional derivative of $\mathcal{A}$ with respect to ${\mathbf{\upsilon}}$. The functional Jacobian is given by

[TABLE]

Although the matrix ${\mathcal{J}}$ in (2) is not skew-symmetric, the skew-symmetry of the Poisson bracket appears after integrations by parts [35], and the Poisson bracket satisfies the Jacobi identity

[TABLE]

for any three functionals $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{D}$. Conservation of the Hamiltonian (3) follows from the antisymmetry of the Poisson bracket (4)

[TABLE]

Other conserved quantities of the SWE [41] are the Casimirs of the form

[TABLE]

where $G$ is an arbitrary function of the potential vorticity $q$. The Casimirs are additional constants of motion which commute with any functional $\mathcal{A}$, i.e., the Poisson bracket vanishes

[TABLE]

Important special cases of the Casimirs are potential enstrophy ${\mathcal{Z}}$, vorticity ${\mathcal{V}}$ and mass ${\mathcal{M}}$, given by

[TABLE]

2.2 SWE in f-plane

In almost all models of the ocean and atmosphere, for the Coriolis force, only the component of the planetary rotation vector that is locally normal to geopotential surfaces is retained. This approach is also known as ”traditional approximation” [44], when the earth’s surface is considered locally flat with latitudinal variation in the Coriolis force. This leads to the so-called $\beta$-plane approximation with $f\approx f_{0}\beta y$ which is used frequently in the literature (see for example [8, 16]). When the region is assumed to be small enough, the latitudinal variation in the Coriolis parameter can be ignored. Then the Coriolis force $f$ is replaced by a constant representative value, which corresponds to the f-plane formulation of the SWE.

Inserting the potential vorticity $q$ in (1), the SWE can be written in the f-plane as

[TABLE]

3 Full-order models

For the space discretization of the SWEs (2) and (5), we form the uniform grid of the spatial domain $\Omega=(a,b)\times(c,d)$ with the grid nodes ${\mathbf{x}}_{ij}=(x_{i},y_{j})^{T}$, where $x_{i}=a+(i-1)\Delta x$ and $y_{j}=c+(j-1)\Delta y$, $i=1,\ldots,N_{x}+1$, $j=1,\ldots,N_{y}+1$, i.e., we divide the domain into $N_{x}$ and $N_{y}$ equidistant parts in $x$ and $y$-directions, respectively, with the mesh sizes $\Delta x=(b-a)/N_{x}$ and $\Delta y=(d-c)/N_{y}$. Then, we define the time-dependent semi-discrete solution vectors in the following order

[TABLE]

where each $u_{ij}(t)$, $v_{ij}(t)$ and $h_{ij}(t)$ denotes the approximation of the solutions $u({\mathbf{x}},t)$, $v({\mathbf{x}},t)$ and $h({\mathbf{x}},t)$, respectively, at the grid nodes ${\mathbf{x}}_{ij}$ at time $t$, $i=1,\ldots,N_{x}$, $j=1,\ldots,N_{y}$. Note that each semi-discrete solution vector in (6) has $N:=N_{x}N_{y}$ components since we omit the solutions on the most right and most top boundary because of the periodic boundary conditions. Throughout the paper, we do not explicitly represent the time dependency of the semi-discrete solutions for simplicity, and we write ${\mathbf{u}}$, ${\mathbf{v}}$ and ${\mathbf{h}}$, and the semi-discrete vector for the solution vector $z$ is defined by ${\mathbf{z}}=({\mathbf{u}};{\mathbf{v}};{\mathbf{h}})\in\mathbb{R}^{3N}$. Similarly, the semi-discrete vector for the potential vorticity $q$ and the Coriolis force $f$ are defined by ${\mathbf{q}}$ and ${\mathbf{f}}$, respectively.

For the approximation of the first-order partial derivative terms, we use one-dimensional central finite-differences to the first-order derivative terms in either $x$ or $y$-direction, and we extend them to two dimensions utilizing the Kronecker product. For a positive integer $s$, let $\widetilde{D}_{s}$ denotes the matrix of central finite-differences to the first-order ordinary differential operator under periodic boundary conditions

[TABLE]

Then, on the two-dimensional mesh, the central finite-difference matrices corresponding to the first-order partial derivative operators $\partial_{x}$ and $\partial_{y}$ are given respectively by

[TABLE]

where $\otimes$ denotes the Kronecker product, and $I_{N_{x}}$ and $I_{N_{y}}$ are the identity matrices of size $N_{x}$ and $N_{y}$, respectively.

We further partition the time interval $[0,T]$ into $N_{t}$ uniform intervals with the step size $\Delta t=T/N_{t}$ as $0=t_{0}<t_{1}<\ldots<t_{N_{t}}=T$, and $t_{k}=k\Delta t$, $k=0,1,\ldots,N_{t}$. Then, we denote by ${\mathbf{u}}^{k}={\mathbf{u}}(t_{k})$, ${\mathbf{v}}^{k}={\mathbf{v}}(t_{k})$ and ${\mathbf{h}}^{k}={\mathbf{h}}(t_{k})$ the full discrete solution vectors at time $t_{k}$.

By the above formulations, the full discrete form of the conserved quantities are given as:

• •

energy

[TABLE]

• •

potential enstrophy

[TABLE]

• •

vorticity and mass

[TABLE]

We can now give the FOM formulations for the SWEs (2) and (5).

3.1 FOM for the SWE in Hamiltonian form

The semi-discrete formulation of the SWE (2) is given as the following $3N$-dimensional system of Hamiltonian ODEs

[TABLE]

where $\circ$ denotes element-wise or Hadamard product. In (10), the matrix function ${\mathbf{q}}^{d}$ of size $N\times N$ stands for the diagonal matrix with the diagonal elements ${\mathbf{q}}^{d}_{ii}={\mathbf{q}}_{i}$ where ${\mathbf{q}}$ is the semi-discrete vector of the potential vorticity $q({\mathbf{x}},t)$, $i=1,\ldots,N$. Time integration of the ODE system (10) by the AVF method [15] yields

[TABLE]

where the integral term is computed explicitly as

[TABLE]

The AVF method [15] preserves higher-order polynomial Hamiltonians, including the cubic Hamiltonian $H$ of the SWE (2). Quadratic Casimir functions like mass and circulation are preserved exactly by AVF method. But higher-order polynomial Casimirs like the enstrophy (cubic) can not be preserved. Practical implementation of the AVF method requires the evaluation of the integral on the right-hand side of (11). Since the discrete Hamiltonian $H$ and the discrete form of the Casimirs, potential enstrophy, mass and circulation are polynomial, they can be exactly integrated with a Gaussian quadrature rule of the appropriate order. The AVF method is used with finite element discretization of the rotational SWE [4, 46] and for thermal SWE [21] in Poisson form. Although the fully implicit integrator gives the desired properties such as conservation of the energy, it is computationally expensive. Semi-implicit implementation of the AVF method with a simplified Jacobian with a quasi-Newton solver is used for the thermal SWE in [21]. In [44], the SWE as a Poisson system is discretized in space following [41] by the Arakawa-Lamb discretization [1] in space. It was shown that by using the time discretization with the non-structure preserving time integrators like the Adams-Bashforth and Runge-Kutta methods, drifts occur in the energy and enstrophy [44]. A general approach for constructing schemes that conserve energy and potential enstrophy formulates the SWEs Nambu brackets [42], which is computationally expensive. Also, it is not possible to preserve multiple integrals like the enstrophy at the same time [19] by geometric integrators like the AVF method. The SWE in a different Poisson form is discretized in space by compatible finite elements and is solved in time using the AVF method in [4, 46].

3.2 FOM for the SWE in f-plane

Semi-discretization of the SWE (5) leads to the following ODE system

[TABLE]

The ODE system (12) consists of quadratic and linear parts, and can be rewritten as

[TABLE]

where the quadratic vector fields $R_{1}(\mathbf{z})$ and $R_{2}(\mathbf{z})$ consist of the quadratic terms, while $L(\mathbf{z})$ contains the linear terms, given by

[TABLE]

For the linear-quadratic autonomous ODE systems (13), Kahan introduced an ”unconventional” discretization [27]

[TABLE]

where the symmetric bilinear forms $R_{1}^{f}(\cdot,\cdot)$ and $R_{2}^{f}(\cdot,\cdot)$ are computed by the polarization [11] of the quadratic vector fields $R_{1}(\cdot)$ and $R_{2}(\cdot)$, respectively, defined by

[TABLE]

Kahan’s method is second-order and time-reversal and linearly implicit, i.e., it requires only one Newton iteration per time step [28]

[TABLE]

where $F^{\prime}({\mathbf{z}}^{k})$ stands for the Jacobian matrix of $F({\mathbf{z}})$ evaluated at ${\mathbf{z}}^{k}$, and $I_{3N}$ is the identity matrix of size $3N$. The Jacobian matrix $F^{\prime}({\mathbf{z}})$ can be computed analytically as

[TABLE]

where the superscript $d$ stands for the matricization of a vector such that for any vector ${\mathbf{a}}\in\mathbb{R}^{N}$, the matrix ${\mathbf{a}}^{d}\in\mathbb{R}^{N\times N}$ is the diagonal matrix with the diagonal elements ${\mathbf{a}}^{d}_{ii}={\mathbf{a}}_{i}$.

A different polarization for polynomial Hamiltonians was introduced in [18] leading to linearly implicit discrete gradient methods. But this method is restricted to non-canonical Hamiltonian systems with constant Poisson matrix. Therefore it can not be applied to the rotational SWE (2) with a state-dependent Poisson matrix. Recently in [37] a two-step version of Kahan’s method is developed for non-canonical Hamiltonians with constant Poisson matrix that preserves Hamiltonian. But this method is also not applicable to the rotational SWE due to state-dependent Poisson matrix.

4 Reduced-order models

In this section, we give the construction of two ROMs for the rotational SWE. The first one preserves the Hamiltonian for the SWE (2) in Poisson form using POD/DEIM. The second one preserves the linear-quadratic structure of the SWE (5) in the f-plane representation by applying TPOD.

4.1 Proper Orthogonal Decomposition

Construction of ROMs is based on approximately representing the solutions of high fidelity FOMs ((10) or (12)) from a low-dimensional linear subspace. In this study, we form the low-dimensional linear subspace from the so-called snapshot matrix of the FOM solutions applying POD, i.e., from the column space of the matrix whose columns are the discrete solution vectors at discrete time instances. In reduced-order modeling of fluid dynamics problems, to avoid the fluctuations, the mean centered snapshots are collected at time instances $t_{k}$, $k=1,\cdots,N_{t}$, in the snapshot matrices $S_{u}$, $S_{v}$ and $S_{h}$

[TABLE]

where $\overline{\mathbf{u}}$, $\overline{\mathbf{v}}$, $\overline{\mathbf{h}}\in\mathbb{R}^{N}$ denote the mean of the snapshots defined by

[TABLE]

For a positive integer $n\ll N_{t}$, the POD aims to form the $n$-POD basis matrices $V_{u,n},V_{v,n},V_{h,n}\in\mathbb{R}^{N\times n}$ which approximately span the column space of the snapshot matrices $S_{u}$, $S_{v}$ and $S_{h}$, respectively. The $n$-POD basis matrices are computed through the application of the singular value decomposition (SVD) to the snapshot matrices

[TABLE]

where the columns of the orthonormal matrices $V_{i}\in\mathbb{R}^{N\times N_{t}}$ and $W_{i}\in\mathbb{R}^{N_{t}\times N_{t}}$ are the left and right singular vectors of the snapshot matrix $S_{i}$, respectively, and the diagonal matrix $\Sigma_{i}\in\mathbb{R}^{N_{t}\times N_{t}}$ with the diagonal elements $(\Sigma_{i})_{jj}=\sigma_{i,j}$, $j=1,\ldots,N_{t}$, contains the singular values of $S_{i}$, $i\in\{u,v,h\}$. Then, the $n$-POD basis matrix $V_{i,n}$ consists of the first $n$ left singular vectors (POD modes) from $V_{i}$, $i\in\{u,v,h\}$. The $n$-POD basis matrix minimizes the following least squares error

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm. Thus the error in the snapshot representation is given by the sum of the squares of the singular values corresponding to those left singular vectors which are not included in the $n$-POD basis matrix. The singular values provide quantitative guidance for choosing the size of the POD basis, that accurately represents the given snapshot data. Usually, the following relative ”cumulative energy” criterion is used

[TABLE]

where $\kappa$ is a user-specified tolerance, typically taken to be $\kappa\leq 10^{-3}$ to catch at least $99.9\%$ of data information.

After computing the $n$-POD basis matrices $V_{u,n}$, $V_{v,n}$ and $V_{h,n}$, we can construct reduced approximations $\widehat{\mathbf{u}}$, $\widehat{\mathbf{v}}$ and $\widehat{\mathbf{h}}$ to the FOM solutions ${\mathbf{u}}$, ${\mathbf{v}}$ and ${\mathbf{h}}$, respectively, which are given by

[TABLE]

where ${\mathbf{u}}_{r}$, ${\mathbf{v}}_{r}$, ${\mathbf{h}}_{r}\in\mathbb{R}^{n}$ are the solution vectors of the reduced system of dimension $n\ll N$. For convenience, we also define the following vectors and matrix

[TABLE]

and then we get the approximation ${\mathbf{z}}\approx\widehat{\mathbf{z}}=\overline{\mathbf{z}}+V_{z,n}{\mathbf{z}}_{r}$. We also note that since the columns of the $n$-POD matrices are orthonormal, the reverse of the approximation identity holds as ${\mathbf{z}}_{r}=V_{z,n}^{T}(\widehat{\mathbf{z}}-\overline{\mathbf{z}})$.

4.2 ROM for the SWE in Hamiltonian form

The semi-discretized SWE (10) in Hamiltonian form can be written as a nonlinear ODE of the form

[TABLE]

By substituting the identity ${\mathbf{z}}\approx\widehat{\mathbf{z}}=\overline{\mathbf{z}}+V_{z,n}{\mathbf{z}}_{r}$ into (15), applying the Galerkin projection onto $V_{z,n}$, and using the orthonormality fact $V_{z,n}^{T}V_{z,n}=I_{3n}$, the following $3n$-dimensional POD reduced model is obtained

[TABLE]

Although the dimension of the ROM (16) is $3n\ll 3N$, the cost for evaluating the nonlinear vector $F(\cdot)$ scales not only with the reduced dimension $n$ but also with the dimension $N$ of the FOM. The computational cost is reduced by sampling the nonlinearity $F(\cdot)$ and then interpolating, known as hyper-reduction technique. Several hyper-reduction methods are developed to reduce the computational cost of evaluating the reduced nonlinear terms at selected points. We use DEIM [14] which is the most frequently used hyper-reduction methods. By the DEIM, full discrete nonlinear vectors are collected in the snapshot matrices defined by

[TABLE]

where $F_{i}^{k}\in\mathbb{R}^{N}$ denotes the $i$-th component of the nonlinear vector $F({\mathbf{z}})$ in (15) at time $t_{k}$, $k=1,\ldots,N_{t}$. Then, we can approximate each $F_{i}({\mathbf{z}})$ in the column space of the snapshot matrix $G_{i}$. We first apply POD to the snapshot matrix $G_{i}$ and find the basis matrix $V_{F_{i},m}\in\mathbb{R}^{N\times m}$ whose columns are the basis vectors (DEIM modes) spanning the column space of the snapshot matrix $G_{i}$. Then, we apply the DEIM algorithm [14] to find a projection matrix $P_{i}\in\mathbb{R}^{N\times m}$ so that we have the approximation

[TABLE]

and then we get the DEIM approximation to the reduced nonlinearities in (16) as

[TABLE]

where $\mathcal{V}_{u,1}=V_{u,n}^{T}V_{F_{1},m}(P_{1}^{T}V_{F_{1},m})^{-1}$, $\mathcal{V}_{v,2}=V_{v,n}^{T}V_{F_{2},m}(P_{2}^{T}V_{F_{2},m})^{-1}$ and $\mathcal{V}_{h,3}=V_{h,n}^{T}V_{F_{3},m}(P_{3}^{T}V_{F_{3},m})^{-1}$ are all the matrices of size $n\times m$, and they are precomputed in the offline stage. Note that the terms $P_{i}^{T}F_{i}(\widehat{\mathbf{z}})\in\mathbb{R}^{m}$ means that we need to compute only $m\ll N$ entry of the nonlinear vector $F_{i}(\widehat{\mathbf{z}})$, $i=1,2,3$. In addition, the computational complexity for the Jacobian matrix reduces from ${\mathcal{O}}(N^{2})$ to ${\mathcal{O}}(nN)$. By using DEIM approximation, the semi-discrete ROM (16) with POD/DEIM becomes

[TABLE]

The ROM (17) is also integrated in time by the AVF method similar to the scheme used for the FOM (10).

In the case of the selection of the number $m$ of DEIM modes, the ”cumulative energy” criterion (14) is used. But, because of the nature of the nonlinearity, snapshot matrices are more sensitive and a larger number of modes are needed for accurate approximation. We take $\kappa\leq 10^{-4}$ to catch at least $99.99\%$ of data.

4.3 ROM for the SWE in f-plane

For PDEs and ODEs with polynomial nonlinearity, ROMs do not require approximating the nonlinear function $F(\cdot)$ through sampling, the reduced-order operators can be precomputed in the offline stage. This is beneficial because the offline-online computation is separated in contrast to the hyper-reduction methods. In the past, for the Navier-Stokes [24, 26] equation, the quadratic polynomial forms of the FOMs are exploited by constructing reduced models. For polynomial nonlinearity, this avoids the approximation of the nonlinear terms by hyper-reduction and allows separation of offline and online computation of FOM and ROM, and the reduced model preserves the linear-quadratic structure of the original system.

In order to avoid the approximation of the nonlinear terms in the SWE (13) by hyper-reduction, such as DEIM, we represent the system in terms of tensors and Kronecker product. On the other hand, mathematical operations with tensors can be easily performed using their corresponding matrix representations. Unfolding a tensor into a matrix is called matricization of a tensor. A common matricization of a tensor ${\mathcal{Q}}$ is the $\mu$-mode matricization $Q^{(\mu)}$ [5]. For a third-order tensor in quadratic systems like (13), there are three different ways of unfolding, depending on the $\mu$-mode that are used for the unfolding. For example $Q^{(1)}$ is called $1$-mode matricization of ${\mathcal{Q}}$. $2$-mode and $3$-mode matricizations can be constructed similarly [5]. We refer to [30] for more details on these basic concepts of tensors. Another advantage of matricizations is that tensor-matrix multiplications can be performed by matrix-matrix products. The special structure of the matrix $Q$ which represents the Hessian of the right-hand side of (13) has been exploited numerically in an efficient way to construct the reduced-order Hessian [5]. The Hessian matrix $Q$ is an unfolding of a $3$-tensor ${\mathcal{Q}}\in\mathbb{R}^{3N\times 3N\times 3N}$. It is sparse due to the local structure of common discretization like finite-differences, finite elements for polynomial nonlinearities. For the SWE with quadratic nonlinear terms, the cross terms $\mathbf{z}_{i}\cdot\mathbf{z}_{j}$ vanish for $|i-j|>3$ in the semi-discretized ODE (13). Therefore the number of nonzero columns of $Q$ is only $2N$.

The semi-discrete system (13) with quadratic and linear parts can be rewritten using the Kronecker product as the following

[TABLE]

with the following redefined quadratic parts

[TABLE]

where $\otimes$ denotes the Kronecker product, and the matrices $A^{x},A^{y},B^{x},B^{y}\in\mathbb{R}^{3N\times 3N}$ are given by

[TABLE]

where $I_{N}$ denotes the identity matrix of size $N\times N$. In (19), the matrix $Q\in\mathbb{R}^{3N\times(3N)^{2}}$ represents the matricized $3$-tensor such that $Q(\mathbf{z}\otimes\mathbf{z})=\mathbf{z}\circ\mathbf{z}$ is satisfied. By substituting the identity ${\mathbf{z}}\approx\widehat{\mathbf{z}}=\overline{\mathbf{z}}+V_{z,n}{\mathbf{z}}_{r}$ into (18), and applying Galerkin projection onto $V_{z,n}$, we obtain the following reduced linear-quadratic equation

[TABLE]

where we set

[TABLE]

In the reduced linear-quadratic equation (20), the tricky part is the computation of the reduced quadratic parts $F_{r}^{u}(\widehat{\mathbf{z}})$ and $F_{r}^{v}(\widehat{\mathbf{z}})$. We next discuss the computation of the reduced quadratic part $F_{r}^{u}(\widehat{\mathbf{z}})$, and then the other reduced quadratic part $F_{r}^{v}(\widehat{\mathbf{z}})$ can be computed in a similar way. Substituting the identity $\widehat{\mathbf{z}}=\overline{\mathbf{z}}+V_{z,n}{\mathbf{z}}_{r}$ into the term $F_{r}^{u}(\widehat{\mathbf{z}})$, and using the properties of the Kronecker product operation, we obtain

[TABLE]

In (21), only the last term in the last line is quadratic and the other terms are at most linear. The quadratic term is computed by the Kronecker product as follows

[TABLE]

where the matrix $Q_{u,r}=-V_{z,n}^{T}A^{x}Q\left(V_{u,n}^{*}\otimes(B^{x}V_{z,n})\right)\in\mathbb{R}^{3n\times(3n)^{2}}$ can be computed in the offline stage, and

[TABLE]

The main computational burden is computation of the Kronecker product $V_{u,n}^{*}\otimes(B^{x}V_{z,n})\in\mathbb{R}^{(3N)^{2}\times(3n)^{2}}$ in $Q_{u,r}$, which has complexity of order ${\mathcal{O}}(n^{2}N^{2})$ for quadratic nonlinearity, and computation of $Q_{u,r}$ in (22) is costly due to the dense structure of the POD basis matrices. In [5], an algorithm is developed to construct the reduced matricized tensor $Q_{u,r}$ for quadratic nonlinearity which avoids the computation of the Kronecker product $V_{u,n}^{*}\otimes(B^{x}V_{z,n})$, having a complexity of order ${\mathcal{O}}(nN^{2})$. Therein, using the $\mu$-mode (matrix) product the reduced matricized tensor can be efficiently computed as following

• •

Compute $\mathcal{Y}^{3n\times 3N\times 3N}$ by $Y^{(1)}=-V_{z,n}^{T}A^{x}Q$,

• •

Compute $\mathcal{Z}^{3n\times 3n\times 3N}$ by $Z^{(2)}=V_{z,n}^{T}(B^{x})^{T}Y^{(2)}$,

• •

Compute $\mathcal{Q}_{u,r}^{3n\times 3n\times 3n}$ by $Q_{u,r}^{(3)}=(V_{u,n}^{*})^{T}Z^{(3)}$.

Although the $\mu$-mode (matrix) product decreases the complexity of evaluating the reduced matrix ${Q}_{u,r}$, still the matrix ${Q}$ has to be built for each different polynomial nonlinearity. Recently two new algorithms [6, 7] are developed for more efficient computation of the reduced matrix ${Q}_{u,r}$ using the particular structure of Kronecker product. A more compact form of the evaluation of the reduced matrix ${Q}_{u,r}$ is given in MatLab notation as follows [6, 7]

[TABLE]

where $G=B^{x}V_{z,n}\in\mathbb{R}^{3N\times 3n}$ and the complexity of this operation is $\mathcal{O}(Nn^{3})$. Thus, the reduced matrix ${Q}_{u,r}$ can be constructed without explicitly defining the matrix $Q$. The transpose of the Kronecker products of any given two vectors $\mathbf{a}$ and $\mathbf{b}$ can be represented as follows

[TABLE]

where vec $(\cdot)$ denotes vectorization of a matrix. Using (24), the matrix $N:=Q(V_{u,n}^{*}\otimes G)\in\mathbb{R}^{3N\times(3n)^{2}}$ in (23) can be constructed as follows

[TABLE]

In [7], a pseudo-skeletal matrix decomposition [36], CUR, is used to increase further computational efficiency of the algorithm above. Instead of the CUR matrix decomposition, here we use ”MULTIPROD” [32] in order to reduce the complexity of the reduced nonlinear terms. MULTIPROD uses virtual array expansion to perform multiple matrix products. When the matrix $V_{u,n}^{*}\in\mathbb{R}^{3N\times 3n}$ is reshaped as $V_{u,n}^{*}\in\mathbb{R}^{3N\times 1\times 3n}$, then MULTIPROD is applied to $G$ and $V_{u,n}^{*}$ in $2$ and $3$ dimensions. MULTIPROD assigns virtually a singleton to the third dimension of $G$, and we get the $3$-dimensional array (tensor) $\mathcal{N}:=\text{MULTIPROD}(G,V_{u,n}^{*})\in\mathbb{R}^{3N\times 3n\times 3n}$. Thus, we can represent (23) as $Q_{u,r}=-V_{z,n}^{T}A^{x}N^{(1)}$, where $N^{(1)}\in\mathbb{R}^{3N\times(3n)^{2}}$ is the matricization of $\mathcal{N}$. In Section 5, we compare the computational efficiency of computing the reduced matricized tensor $Q_{u,r}$ by the algorithm in [5] with the algorithms in [6, 7] improved by the use of MULTIPROD.

5 Numerical results

We consider a test example for the SWE (1) on the spatial domain $\Omega=[0,1]^{2}$, with $g=1$, $f=0$, and with the initial conditions [45]

[TABLE]

The initial wave satisfy the periodic boundary conditions. The final time is set to $T=50$, and spatial and temporal mesh sizes are taken as $\Delta x=0.01$ and $\Delta t=4\Delta x$, respectively. This leads to a spatial grid of size $N=10000$, and $N_{t}=1250$ time intervals, so that each snapshot matrix $S_{u}$, $S_{v}$ and $S_{h}$ has size $10000\times 1250$. For the FOM simulations, we consider the solutions of the SWE (10) in the Hamiltonian form with AVF time integrator (SWE-AVF), and the SWE (12) in the f-plane with Kahan’s time integrator (SWE-Kahan). In case of ROMs, related to the FOM (10) and by AVF time integrator, we consider the ROM (16) without DEIM approximation (POD-AVF) and the ROM (17) with DEIM approximation (POD-AVF-DEIM). On the other hand, related to the FOM (12) and by Kahan’s time integrator, we consider the ROM (20) without tensorial framework (POD-Kahan) and with tensorial framework (TPOD-Kahan). All the simulations are performed on a machine with Intel*®* CoreTM i7 2.5 GHz 64 bit CPU, 8 GB RAM, Windows 10, using 64 bit MatLab R2014.

The full order solutions of SWE-AVF and SWE-Kahan are given in Figure 1, depicting that the solutions by both approaches are in good agreement. In Figure 2, we give the normalized singular values of the snapshot matrices $S_{u}$, $S_{v}$ and $S_{h}$ related to the velocity components $\mathbf{u}$, $\mathbf{v}$ and the height $\mathbf{h}$, respectively. The singular values decay slowly for each snapshot matrix obtained by both SWE-AVF and SWE-Kahan, which is the characteristic of the problems with wave phenomena in fluid dynamics [40]. They pose a challenge for reduced-order methods, since their dynamical behavior cannot be captured accurately by the linear combination of a few POD modes.

In order to show conservation of the discrete energy (7), discrete enstrophy (8), and discrete vorticity (9) for a FOM solution vector $\mathbf{z}$ (or ROM solution vector $\widehat{\mathbf{z}}$), we use the time-averaged absolute errors $\|\cdot\|_{H}$, $\|\cdot\|_{Z}$ and $\|\cdot\|_{V}$ defined by

[TABLE]

On the other hand, to measure the global error between a discrete FOM solution vector and a discrete ROM approximation (FOM-ROM error), we define the following time averaged relative $L_{2}$ errors for the state variables ${\mathbf{w}}\equiv{\mathbf{u}},{\mathbf{v}},{\mathbf{h}}$

[TABLE]

and to plot the FOM-ROM error at a specific time instance, we simply take the node-wise difference between the full discrete FOM solution and ROM approximation.

Due to the slow decay of the singular values, FOM-ROM errors for all state variable with varying number of POD modes in Figure 3 decrease with small oscillations. Therefore, for all the simulations, the number of POD modes is set to a relatively large value, $n=50$, according to the relative energy criteria (14) with $\kappa=10^{-4}$. Using the same energy criteria, the number of DEIM modes is set to $m=90$.

In Figure 4, the FOM-ROM errors at the final time for the height are at the same level of accuracy for both approaches, whereas the FOM-ROM errors of the velocity components by the TPOD-Kahan are smaller than the ones by the POD-DEIM-AVF, which might be due to differences in the solution algorithms.

In Figures 5-7, the discrete energy error $|H^{k}-H^{0}|$, the discrete enstrophy error $|Z^{k}-Z^{0}|$, and the discrete vorticity error $|V^{k}-V^{0}|$ are plotted, $k=1,\ldots,N_{t}$. The discrete energy, enstrophy, and vorticity are well preserved over the time interval. Because the vorticity is a quadratic conserved quantity, it is well preserved by both methods in Figure 7. The energy and the enstrophy errors in Figures 5-6 show small drifts for POD-DEIM-AVF, but all they have bounded oscillations over the time, i.e., they are preserved approximately at the same level of accuracy. The discrete energy and the discrete enstrophy is better preserved by TPOD-Kahan, since the reduced nonlinearity is not approximated by hyper-reduction. The mass is preserved up to machine precision since it is a linear conserved quantity, and it is not shown here.

The time-averaged relative $L_{2}$-errors between FOM and ROM solutions in Table 1 are at the same level of accuracy for both versions of ROMs. The ROMs for the height are more accurately resolved than for the velocity components. The conserved quantities are more accurately preserved by the POD-AVF and by the TPOD-Kahan than by the POD-DEIM-AVF in Table 2. This indicates that using DEIM, the stability of the ROM solutions can not be guaranteed in long term integration, which is reflected in the preservation of the conserved quantities.

The CPU times and speedup factors in Table 3 for $n=50$ POD and $m=90$ DEIM modes, show the computational efficiency of the TPOD-Kahan over POD-DEIM-AVF due to the separation of offline-online computation. In Table 3, basis computation includes SVD computation, and online computation consists of time needed for projection and solution of reduced system. As pointed out in [31] that for some problems, a large number of DEIM interpolation points are required to achieve accurate solutions, which increases the computational cost of the ROMs in online computation, as in Table 3. The computational efficiency is further increased by exploiting the sparse matrix structure of the discretized SWE using MULTIPROD (MP) in the algorithm of [6, 7] over [5] as shown in Table  3.

For the reduced matricized tensor $Q_{u,r}=-V_{z,n}^{T}A^{x}Q\left(V_{u,n}^{*}\otimes(B^{x}V_{z,n})\right)$, the computational time by the two-sided projection method in [5] (TS) is compared with the algorithms in [6, 7] using MULTIPROD. In Figure 8, left, we give the required computational time versus the number of POD modes by fixing the number of grid points $N=10000$, whereas the required computational time for varying number of grid points using a fixed number of POD mode $n=50$ are given in Figure 8, right. Both figures show the computational efficiency using MULTIPROD by increasing the size of the FOM and ROMs.

6 Conclusions

We have shown that the two kind of ROMs preserve the conserved quantities of the SWE and yield stable reduced solutions in long time integration. For both approaches, the reduced solutions and the reduced conserved quantities have almost the same level of accuracy as the full ones. Using TPOD-Kahan and exploiting the quadratic structure of the SWE, the online computation time of ROMs can be much reduced compared to the POD-DEIM-AVF. The online computation time of the ROMs is further reduced by applying the algorithms in [6, 7] in combination with the MULTIPROD [32]. In a future work, we will investigate both approaches for the SWE with full Coriolis force [44] and thermal SWE [21].

ACKNOWLEDGEMENT

The authors thank for the constructive comments of the referees which helped much to improve the paper. This work has been supported by 100/2000 Ph.D. Scholarship Program of the Turkish Higher Education Council.