The shallow water equations (SWEs) consist of a set of two-dimensional equations describing a thin inviscid fluid layer flowing over the topography in a frame rotating about an arbitrary axis. SWEs are widely used in modeling of large-scale atmosphere/ocean dynamics, in 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 Cotter18 ; Stewart16 . Therefore numerical schemes that preserve energy and enstrophy lead to stable solutions in long term integration. Salmon Salmon04 introduced the non-canonical Hamiltonian/Poisson form of the SWE in the rotation frame with constant Coriolis force and constructed an energy preserving scheme in space using Arakawa-Lambs discretization Arakawa81 . The SWEs with the complete Coriolis force are introduced in Salmon05a
and share the same continuous skew-symmetric structure and Poisson bracket. The discrete energy conservation follows from antisymmetry of the discrete Poisson bracket, other conserved quantities are potential enstrophy, mass and vorticity, i.e. the Casimirs.
Time discretization schemes that preserve geometric structure like the symplecticness, energy or entropy lead to stable approximate solutions of dynamical systems in long time integration Hairer10 ; Reich05 . In this paper, we discretize SWE using finite differences in space by preserving the skew-symmetry Poisson matrix. We apply two different time discretization methods for the semi-discretized SWE. The first one is the energy preserving fully implicit, second order convergent average vector field (AVF) method Celledoni12ped ; Cohen11 . The AVF method belongs to the class of the discrete gradient method for ODEs Dahlby11 ; McLachlan99 ; Quispel96
. Alternatively the semi-discretization of the SWE in the f-plane form leads to system of ordinary differential equations (ODEs) with quadratic nonlinearities. We use second order linearly implicit Kahan’s methodKahan93 for time integration. SWE in the f-plane Celledoni13 is a linear-quadratic ordinary differential equation (ODE) after semi-discretization by finite differences in space. Kahan’s ”unconventional” discretization method applied to ODEs with quadratic vector fields, preserves the integrals or conserved quantities of many Hamiltonian and integrable systems. Fully implicit schemes like the AVF method require iterative solvers like Newton’s method at each time step for solving the nonlinear systems. Linearly implicit schemes like Kahan’s method require the solution of precisely one linear system of equations in each time step.
Numerical methods for solving partial differential equations (PDEs) are computationally costly and require a large amount of computer memory for applications in real-time or for evaluations for many different parameters. During the last decades, reduced order methods (ROMs) have emerged as a powerful approach to reduce the cost of evaluating large systems of PDEs by constructing a low-dimensional linear reduced subspace, that approximately represents the solution to the system of PDEs with a significantly reduced computational cost. The proper orthogonal decomposition (POD) with the Galerkin projection into lower dimensional space, has been widely used as a computationally efficient surrogate model in large-scale numerical simulations of nonlinear PDEs. The high fidelity solutions are generated by numerical simulations of the discretized high dimensional full order model (FOM). The POD is then applied with the Galerkin projection to compute an optimal subspace to fit the high fidelity data. In the online stage, the reduced system is solved in the low-dimensional subspace. The primary challenge in producing the low dimensional models of the high dimensional discretized PDEs is the efficient evaluation of the nonlinearities (inner products) on the POD basis. The computational cost is reduced by sampling the nonlinear terms and interpolating, known as hyper-reduction techniques. Several hyper-reduction methods are developed to reduce the computational cost of evaluating the reduced nonlinear terms, empirical interpolation method (EIM) Barrault04 , discrete empirical interpolation method (DEIM) chaturantabut10nmr
, missing point estimationAstrid08 ; Zimmermann16 , best points interpolation Nguyen08 , gappy POD Carlberg13 . We use DEIM, which is one of the most frequently used hyper-reduction methods. Rigorous error and stability analysis are in general not amenable for subsampling procedures used in hyper-reduction methods. For PDEs and ODEs with polynomial nonlinearities, ROMs do not require approximating the nonlinear terms through sampling, the reduced order operators can be precomputed. Projection of FOM onto the reduced space yields low-dimensional matrix operators that preserve the polynomial structure of the FOM. This is an advantages because the offline-online computation is separated in contrast to the hyper-reduction methods. For polynomial nonlinearities, tensorial POD preserve the structure of FOM in the reduced space. Recently, tensorial POD methods are applied for PDEs and ODEs with polynomial nonlinearities Benner15 ; Benner18 ; Benner19 ; Gu11 ; Kramer19 .
The naive application of POD/DEIM and tensorial POD may not preserve the geometric structures, like the symplecticness, energy preservation and passivity of Hamiltonian, Lagrangian and port Hamiltonian PDEs. The stability of reduced models over long-time integration and the structure-preserving properties has been recently investigated in the context of Lagrangian systems Lall03 ; Carlberg15 , and for port-Hamiltonian systems Chaturantabu16 . For linear and nonlinear Hamiltonian systems, the symplectic model reduction technique, proper symplectic decomposition (PSD) is constructed for Hamiltonian systems like, linear wave equation, sine-Gordon equation, nonlinear Schrödinger (NLS) equation to ensure long term stability of the reduced model Hesthaven16 ; Peng16 . Recently the AVF method is used as a time integrator to construct reduced order models for Hamiltonian systems like Korteweg-de Vries equation Gong17 and NLS equation Karasozen18 . Reduced order models for the SWEs are constructed in conservative form using POD-DEIM Lozovskiy16 ; Lozovskiy17 , in the -plane by POD-DEIM and tensorial POD Stefanescu13 ; Stefanescu14 , by dynamic mode decomposition (DMD) Bistrian15 ; Bistrian17 , the f-plane using POD Esfahanian09 . In these articles, the preservation of the energy and other conservative quantities in the reduced space are not discussed. We apply the POD/DEIM for the SWE in Poisson form and integrate in time using the AVF method which preserves the Hamiltonian (energy) in the reduced space. Using the quadratic semi-discretized form of the SWE in the f-plane, we also apply two different tensorial POD algorithms Benner15 ; Benner18 ; Benner19
using the matricization of tensor and exploiting the sparse structure of the discretized SWE. The reduced equations are also solved by the linearly implicit Kahan’s method. A comparison with DEIM shows the efficiency of the tensorial POD by exploiting the quadratic structure of the semi-discretized SWE.
In the next section, we introduce 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 in AVF semi-discretized Hamiltonian form, and tensorial POD for the semi-discrete ODE system with quadratic nonlinearities in Section 4. In Section 5, we compare POD, POD/DEIM, and tensorial POD for two dimensional rotating SWE with respect to the accuracy of the reduced solutions, preservation of the energy, enstrophy, vorticity and computational efficiency by the ROMs. The paper ends with some conclusions.
2 Shallow water equation
Many geophysical flows can be written in Hamiltonian form Morrison98 . The rotational SWE in non-canonical Hamiltonian/Poisson form was described in Salmon04 . Later on the Nambu formulation of the SWE Salmon07 , SWE with complete Coriolis force Salmon05a ; Stewart16 , and multi-layer SWE Stewart16 are developed. The two-dimensional rotational SWE is given as Salmon04
where is the (particle) velocity, is the fluid depth, is the potential vorticity and is the Coriolis force, is the gravity constant, is the state vector and . The non-canonical Hamiltonian/Poisson form of the SWE (1) is given by
with the Hamiltonian (energy)
where , and is the functional derivative of with respect to . The functional Jacobian is given by
The Poisson bracket (4) is related to the skew-symmetric Poisson matrix as . Although the matrix in (2) is not skew-symmetric, the skew-symmetry of the Poisson bracket appears after integrations by parts Lynch02 , and the Poisson bracket satisfies the Jacobi identity
Other conserved quantities are the Casimirs of the SWE Salmon04 of the form
where is an arbitrary function of the potential vorticity . The Casimirs are additional constants of motion, they commute with all other state functions
Important special cases are the potential enstrophy
the mass , and the vorticity . Alternatively, inserting the potential vorticity in (1), the SWE can be written in the f-plane as
3 Full order model (FOM)
We discretize the rotational SWE (1) using finite differences in a rectangular space domain using the uniform grid . We take in both space direction the same mesh sizes , . The matrix denotes the first order central finite difference operators and under periodic boundary conditions
The corresponding first order finite difference approximations are , , where
denotes the Kronecker product. The identity matrix of sizeis denoted by . The semi-discrete form of (2) is given as
where denotes element-wise or Hadamard product. The discretized variables are
where denotes the diagonal vector and .
The vorticity is discretized as
and . Discrete form of the conserved quantities are given by:
The AVF method Cohen11 can be viewed as a generalization of the implicit midpoint rule which preserves quadratic Hamiltonians. It preserves higher order polynomial Hamiltonians, which includes the cubic Hamiltonian . Quadratic Casimirs 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 (7). Since the Hamiltonian 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 degree. The AVF method is used with finite element discretization of the rotational SWE Cotter18 ; Cotter19 and for thermal SWE Eldred19 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 with a simplified Jacobian with a quasi-Newton solver is used for the thermal SWE in Eldred19 . In Stewart16 SWEs as a Poisson systems is discretized in space following Salmon04 by the Arakawa-Lamb discretization Arakawa81 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 Stewart16 . A general approach to constructing schemes that conserve energy and potential enstrophy formulates the SWEs Nambu brackets Salmon05 , which is computationally expensive. Also it is not possible to preserve multiple integrals like the enstrophy at the same time Dahlby11a by geometric integrators like the AVF method. The SWE in a different Poisson form was discretized in space by compatible finite elements and solved in time using the AVF method Cotter18 ; Cotter19 .
Semi-discretization of the SWE (5) leads to the following ODE system
It can be rewritten using the Kronecker product as a quadratic ODE
with , , , and is the identity matrix of size . The FOM (8) consists of a linear part , and quadratic parts . For the quadratic ODE systems like (8), Kahan introduced an ”unconventional” discretization Kahan93 given by
The symmetric bilinear forms and are computed by the polarization of the quadratic vector fields and Celledoni15
is polarized in the same form. Kahan’s method is second order and time-reversal and linearly implicit, i.e requires only one Newton iteration per time step Kahan97 :
where denotes the Jacobian of .
A different polarization for polynomial Hamiltonians were introduced in Dahlby11 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 (1) with state-dependent Poisson matrix. Recently in Miyatake18 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 not applicable to the SWE due to state-dependent Poisson matrix in (2).
4 Reduced order model
In this section, we give the construction of two reduced order models for the SWE. The first one preserves the Hamiltonian for the SWE in Poisson form (1) using POD/DEIM. The second one preserves the linear-quadratic structure for the SWE in the -plane representation (5) by applying tensorial POD.
4.1 Hamiltonian preserving reduced order modelling with POD/DEIM
The semi-discretized SWE in Hamiltonian form (6) represents a nonlinear ODE of the form
where the state variable . In reduced order modeling of fluid dynamics problems to avoid the fluctuations, the mean centered snapshots are collected at time instances , in the snapshot matrices , and Stefanescu14 ; Taira17
where , , , and , and denote the mean of the snapshots
We have collected the velocity components and , and the height separately in snapshot matrices. The velocity components and are stacked in column vectors Taira17
The POD basis functions are computed by applying the (thin) singular value decomposition (SVD)
to the snapshot matrices for the velocity components and the height component separately. The columns of the orthonormal matrices and are the left and right singular vectors of the snapshot matrix , respectively, . The diagonal matrix contains the singular values of . The singular vectors for the velocity components and are selected by unstacking (11). The r-POD basis matrix consists of the first columns of . The POD basis minimizes the following least squares error
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 POD basis. The singular values provide quantitative guidance for choosing the size of the POD basis, that accurately represent the given snapshot data. Usually the following relative ”cumulative energy” criteria is used
where is a user-specified tolerance, typically taken to be or greater. For , the state is then approximated in the POD subspace spanned by the first POD basis as
The ROM (14) with the discrete energy , does not necessarily satisfy the invariance of the discrete energy for all time. But the updated ROM
The cost of evaluating scales not only with the dimension of the ROM, also with the dimension the FOM, . The computational cost is reduced by sampling the nonlinearity a subset of the elements of and interpolating, known as hyper-reduction techniques. Several hyper-reduction methods are developed to reduce the computations cost of evaluating the reduced nonlinear terms, empirical interpolation (EIM) Barrault04 , discrete empirical interpolation (DEIM) chaturantabut10nmr , QDEIM Gugercin16a , missing point estimation Astrid08 , best points interpolation Nguyen08 , gappy POD Carlberg13 . We use the DEIM, which is one the most frequently used hyper-reduction method. Rigorous error and stability analysis are in general not amenable for sampling procedures used in hyper-reduction methods.
We apply the DEIM to the ODE form (10) of SWE together with the AVF time integrator, and the ROM (15) is also integrated by AVF. Thus, in order to proceed with the DEIM, the AVF applied nonlinear snapshots are collected on time intervals , , in the snapshot matrices , ,
where we set , , and each corresponds to the nonlinearity in one of the three equations in the semi-discrete system of SWE. Then, we can approximate each in the column space of the snapshot matrices . We first apply POD to the snapshot matrices and find the matrices whose columns are the basis functions spanning the column space of the snapshot matrices . Then, we apply the DEIM algorithm chaturantabut10nmr to find projection matrices so that we have
where is precomputed in the offline stage, and by the term , we need to compute only entry of the nonlinear vector. In addition, the computational complexity for the Jacobian matrix reduces from to . By using DEIM approximation, the full discrete ROM with POD/DEIM reads as
In case of the selection of the number of DEIM basis functions, the ”cumulative energy” criteria (12) is used, as well. But, because of the nature of the nonlinearity, snapshot matrices are more sensitive and larger number of modes are needed for well reflection, we take now 99.99% or greater.
4.2 Tensorial POD
For PDEs and ODEs with polynomial nonlinearities, ROMs do not require approximating the nonlinear function 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 Graham99 ; Holmes12 equation, the quadratic polynomial forms of the FOMs are exploited by constructing reduced models. For polynomial nonlinearities the ROMs preserve the linear-quadratic structure.
The mean-centered ROM for (9) takes the form
with the precomputed reduced order operators
This avoids the approximation of the nonlinear terms by hyper-reduction and allows separation of offline and online computation of FOM and ROM.
By introducing the reduced order solutions and operators, the linear-quadratic reduced equation (16) becomes
We can compute as follows
The nonlinear term above is computed by the Kronecker product as follows
The other quadratic term can be computed in the same way. In Bistrian15 ; Stefanescu14 the semi-discrete form of SWE in the -plane leads also to a linear quadratic ODE system like (8). They construct the reduced tensor as follows
where and . Then, the reduced matrix can be constructed by the mode-1 unfolding of the reduced tensor . It was shown that the tensorial POD (TPOD) has online complexity and DEIM has , where is the number of interpolation points. In Bistrian15 ; Stefanescu14 TPOD was computed using Frobenious dot products, therefore the TPOD was slower than the DEIM. Below we give two new approaches which reduces the computational cost of the TPOD significantly.
The special structure of the matrix which represents the Hessian of right-hand side of (8) will be exploited in computation of TPOD Benner15 . The matrix is sparse due to the local structure of common discretizations like finite differences, finite elements for homogeneous polynomial nonlinearities. For the SWE with quadratic nonlinear terms, the cross terms vanish for in the semi-discretized ODE (8). Therefore the number of nonzero columns of are only. The matrix is an unfolding of the tensor . Reshaping the tensor into a matrix is called matricization, because working with matrices instead of tensors is beneficial. A common matricization of is the so called mode matricization Benner15 . For example is called mode matricization of . For mode- and mode- matricizations, we refer to Benner15 . Using the matricization of the tensors, tensor matrix multiplication is equivalently determined by matrix-matrix products. The main computational burden in TPOD is computation of , which has complexity of order for quadratic nonlinearities. In Benner15 an algorithm is developed to construct the reduced Hessian which avoids the computation of .
Below we give the TPOD algorithm in Benner15 for the quadratic nonlinear term . In (4.2) the computation of will be inefficient due to the dense structure of the reduced matrix , i.e. . Using the -mode (matrix) product the POD modes can be efficiently computed as following Benner15 :
Compute by ,
Compute by ,
Compute by .
Because the matrix is sparse and has nonzero columns, will not result in a dense matrix, i.e. it has not the complexity of . This enables the efficient computation of the POD modes for large quadratic systems. At each step the tensor multiplication is performed by means of matricization.
The reduced linear-quadratic is solved by Kahan’s method with the Jacobian evaluated using the following tensorial form
Recently a new TPOD approach was developed for a more efficient computation of the reduced matrix using the particular structure of Kronecker product Benner18 ; Benner19 . Although the -mode (matrix) product decreases the complexity of evaluating the reduced matrix , still the matrix has to be build for each different polynomial nonlinearities. A more compact form of the evaluation of the reduced matrix is given in MATLAB notation as follows Benner18 ; Benner19 :
where and the complexity of this operation is . Thus, the reduced matrix can be constructed without explicitly defining the matrix . The transpose of the Kronecker products of any given two vectors and can be represented as follows
where vec denotes vectorization of a matrix. Using (19), the matrix in (4.2) can be constructed as follows . In Benner18 ; Benner19 the so called CUR Mahoney09 , an pseudo-skeletal matrix approximation is used to increase computational efficiency. Here we make use of ”MULTIPROD” leva08mmm to reduce the complexity of the reduced nonlinear terms. MULTIPROD use virtual array expansion to perform multiple matrix products. Now, if we reshape the matrix as and compute MULTIPROD of and in and dimensions. MULTIPROD assigns virtually a singleton to the third dimension of which =MULTIPROD(). Thus, we can represent (4.2) as .
5 Numerical results
We consider the SWE on the spatial domain , with data , and with the initial conditions
Periodic boundary conditions Sato18 are prescribed. The final time is set to , and spatial and temporal mesh sizes are taken as and , respectively. All simulations are performed on a Windows 10 machine with Intel Core i7, 2.5 GHz and 8 GB using MATLAB R2014. CPU time is measured in seconds.
The singular values of the snapshot matrices for both methods are plotted in the semi-logarithmic scale in Figure 1, with the index of singular values on the horizontal axis. The singular values decay for both methods slowly, which is the characteristic of the problems with wave phenomena in fluid dynamics.
The errors in Figure 2 for the state variables with respect to the number of POD modes decrease similarly for both methods, in which the slow decay of the singular values in Figure 1 is reflected to slow decreasing behavior of the ROM errors. The number of POD modes in further computations are selected as , according to the relative energy criteria (12) with . Using the same criteria the number of DEIM interpolation points are selected as .
The FOM solutions for the state variables for both methods are similar in Figure 3. The FOM solutions are captured well by the POD/DEIM with the AVF method, and TPOD with the Kahan’s method in Figure 4. The errors for approximate solutions and conserved quantities are the same by POD and TPOD.
The energy, enstrophy and vorticity are well preserved in Figures 5-7. Because the vorticity is a quadratic Casimir, it shows no drift over time, Figure 7. The energy and enstrophy have cubic terms, therefore they show some drifts (for the POD in the AVF method), but the ROM solutions with POD/DEIM and TPOD have bounded oscillations over the time, i.e. they are preserved with some accuracy. The mass is preserved up to machine precision since it is a linear conserved quantity, and it was not shown.
To illustrate the FOM-ROM state errors, we use time averaged relative -error defined by
and for conservation of the energy, enstrophy, and vorticity, we use the following time-averaged relative FOM-ROM error