Davide Torlo
email: davide.torlo@sissa.it

Mario Ricchiuto

Model order reduction strategies for weakly dispersive waves

Keywords: model order reduction, dispersive wave equations, BBM-KdV, Boussinesq, hyper-reduction 2010 MSC: 65M06, 76B15, 65M22

We focus on the numerical modelling of water waves by means of depth averaged models. We consider in particular PDE systems which consist in a nonlinear hyperbolic model plus a linear dispersive perturbation involving an elliptic operator. We propose two strategies to construct reduced order models for these problems, with the main focus being the control of the overhead related to the inversion of the elliptic operators, as well as the robustness with respect to variations of the flow parameters. In a first approach, only a linear reduction strategies is applied only to the elliptic component, while the computations of the nonlinear fluxes are still performed explicitly. This hybrid approach, referred to as pdROM, is compared to a hyper-reduction strategy based on the empirical interpolation method to reduce also the nonlinear fluxes. We evaluate the two approaches on a variety of benchmarks involving a generalized variant of the BBM-KdV model with a variable bottom, and a one-dimensional enhanced weakly dispersive shallow water system. The results show the potential of both approaches in terms of cost reduction, with a clear advantage for the pdROM in terms of robustness, and for the EIMROM in terms of cost reduction.

Introduction

Water waves equations can be modeled with various strategies [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]. Different models provide reasonable approximations in different contexts, according to which type of solution we are interested in. The most complete models consider Euler or Navier-Stokes equations and describe the motion of the water in each point in the 3 dimensional space. A first approximation level is given by averaging the water speed on the vertical direction, reducing by one the dimensions and allowing to display only the water surface as a function of the horizontal variables. The well-known shallow water equations, Boussinesq equations, Green-Naghdi equations, and other models fall in this category. In this work we are interested in studying dispersive models, which are able of reproducing dispersion phenomena. In contrast with other models, as shallow water equations, in these dispersive models shocks are avoided and instead the dispersion fragments the waves into different waves with different wavelength before they break. We will consider two dispersive models which approximate water waves at different levels: the Benjamin-Bona-Mahony (BBM) equation [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF], in a more general form that it is linked to the Korteweg-De Vries (KdV) equation [START_REF] Chen | An EIM-degradation free reduced basis method via over collocation and residual hyper reduction-based error estimation[END_REF][START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF], and an enhanced Boussinesq (EB) system of equations [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF].

The peculiarity of both models is the combination of hyperbolic systems and dispersive terms. Numerical solutions of such problems are often obtained using explicit solver for the hyperbolic part, while the dispersive terms are obtained first solving an elliptic problem, which requires greater computational costs, and then using this solution inside the hyperbolic system [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF][START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF]. In this work, we will try to reduce the computational costs of such problems, focusing on the compression of the elliptic operators, which are responsible of the largest part of the computational time of such methods. In case of parametric problems, where a fast response is necessary or we have a multi-query task, the reduction could lead to strong advantages for the computational costs, without degrading the quality of the solution.

Hyperbolic problems are known for developing shocks and advection dominated solutions. Classically, they are badly reducible as the Kolmogorov n-width decays very slowly for such problems [START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF][START_REF] Strazzullo | POD-Galerkin model order reduction for parametrized nonlinear time dependent optimal flow control: an application to Shallow Water Equations[END_REF][START_REF] Burman | The edge stabilization method for finite elements in CFD[END_REF][START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF][START_REF] Peherstorfer | Stability of discrete empirical interpolation and gappy proper orthogonal decomposition with randomized and deterministic sampling points[END_REF][START_REF] Orszaghova | From the paddle to the beach-A Boussinesq shallow water numerical wave tank based on Madsen and Sørensen's equations[END_REF]. Nevertheless, for dispersive problems there is no shock formations and the wave traveling is often transformed into a oscillation of the whole domain. Hence, also the advection character of the solution is less pronounced. Hence, there is room to attempt a reduction with classical model order reduction (MOR) algorithms coming from the parabolic and elliptic community. In particular, we will use the proper orthogonal decomposition (POD) [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF] to reduce the solution manifold and then we will apply a Galerkin projection to obtain a reduced problem. A further step will consist of interpolating the nonlinear fluxes with the empirical interpolation method (EIM) [START_REF] Bacigaluppi | Implementation and Evaluation of Breaking Detection Criteria for a Hybrid Boussinesq Model[END_REF]. This will provide a second level of reduction that allows more reduction in the computational costs. Few works have already performed a model order reduction directly on shallow water equations, inter alia [START_REF] Stefanescu | Comparison of POD reduced order strategies for the nonlinear 2D shallow water equations[END_REF][START_REF] Strazzullo | Model reduction for parametrized optimal control problems in environmental marine sciences and engineering[END_REF][START_REF] Ricchiuto | Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF][START_REF] Stefanescu | POD/DEIM nonlinear model order reduction of an ADI implicit shallow water equations model[END_REF]. Up to our knowledge, there are no other works trying to reduce the computational costs of dispersive wave equations by the means of model order reduction techniques.

The paper is structured as follows. In Section 3 we introduce the BBM-KdV model, some energy properties and its classical discretization, i.e., the full order model (FOM), and we introduce few tests that we will study along the paper. In Section 3.5 we develop the reduction algorithms and the hyperreduction steps. In Section 4 we test the presented algorithms on all the presented tests, showing the power and limits of these reduction techniques. In Section 5 we introduce an enhanced Boussinesq system of equations, its FOM discretization and some reference tests that we will use. We introduce then its ROM and its hyper-reduction algorithm. Then we apply the reduced algorithms for the presented tests in Section 6, where we compare the two reduced methods and the FOM one. In Section 7 we highlight the major steps of this work and we suggest possible developments and extensions.

Generalities: dispersive wave models and model reduction

We discuss the general form of the dispersive wave model used in this work, which also allows to explain the main underlying idea. Weakly dispersive wave models can be often seen as perturbations of some hyperbolic partial differential system [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF], and can be written in general form

∂ t u + ∂ x F (u) + S(u) -µ 2 X t (∂ t u + ∂ x G(u)) + µ 2 D x u = 0 (1)
where u : R → R d are the unknowns of the system, F : R d → R d is a (nonlinear) flux which describes the hyperbolic operator, S : R d → R d is a source term which might include bathymetry effects, X t : R d → R d is a linear operator that contains second derivative terms, while D x : R d → R d is a linear operator defined by some dispersion terms with third order derivatives. The first three terms in the above equation define a hyperbolic balance law, as e.g. the shallow water equations with bathymetry, and µ 2 is a small parameter multiplying the weakly dispersive regularization [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]. Two specific examples of this general family are used in the paper: the BBM-KdV model [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Chen | An EIM-degradation free reduced basis method via over collocation and residual hyper reduction-based error estimation[END_REF][START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF] and an enhanced Boussinesq model [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF].

The main idea of this work is to exploit the smallness of µ 2 to defined a hybrid model, which enhances the hyperbolic equation with a reduced order approximation of the small terms. To this end, we follow [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Cagniart | Model Order Reduction for Problems with Large Convection Effects[END_REF] and recast the system as

(1 -µ 2 X t )(∂ t u + ∂ x F (u) + S(u)) + µ 2 X x u = 0, (2)
having set

X x u := X t (∂ x F (u) -∂ x G(u) + S(u)) + D x u. (3)
The models share the property of being composed of two parts: an elliptic (linear) operator and a hyperbolic (nonlinear) operator. The two operators are treated very differently in their discretization, where the hyperbolic part is often discretized in an explicit way, while the elliptic operator is discretized implicitly. Moreover, the combination of the two operators allows to obtain stable schemes with little or zero extra numerical viscosity, which is not the case for pure hyperbolic problems. In particular, we can recast the system introducing an auxiliary variable Φ ∈ R d :

   1 -µ 2 X t Φ + µ 2 X x u = 0, ∂ t u + ∂ x F (u) + S(u) = Φ. (4)
The operator (1 -X t ) is elliptic (and invertible), and by classical formal arguments we can show the smallness of Φ simply by noting that (see also [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF][START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF][START_REF] Lannes | Modeling shallow water waves[END_REF])

Φ = -µ 2 1 -µ 2 X t -1 X x u = O(µ 2 ).
Our main idea is to devise an approach where only the term Φ is reduced, possibly benefiting from the combination of the modelling error (controlled by µ 2 ) and model reduction error. The motivation for this is discussed in the next section.

Time and space discretization of the model

As in [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Cagniart | Model Order Reduction for Problems with Large Convection Effects[END_REF] for all models we use a splitting of the elliptic and hyperbolic operators, as in an IMEX time discretization. This leads to the semi discrete prototype

     (1 -X t ) Φ n+1 + X x u n = 0, u n+1 -u n ∆t + ∂ x F (u n ) + S(u n , u n+1 ) = Φ n+1 . (5)
Here, we denote with n the time discretization index, and the µ 2 factor has been included in the operators X t and X x to lighten the expressions. As we can see, all the operations are vectorial except the solution of Φ n+1 , which requires the solution of a linear system. After a further spatial discretization the system can be written in the general form

     M Φ -X t Φ n+1 + X x u n = 0, M u u n+1 -u n ∆t + F(u n ) + S ex (u n ) + S im u n+1 = Φ n+1 , (6)
where, with an abuse of notation, we use the same symbols u ∈ R d×N h and Φ ∈ R N h to refer to the discretized variables, where N h is the number of degrees of freedom for each variable. The discretization can be performed by finite differences (FD) or finite element methods (FEM). Specific choices used for our simulations are discussed later in the text. M u and M Φ are mass matrices, X t and X x are the discretizations of X t and X x respectively. S ex is the discretization of the explicit part of the source function S and S im is the matrix that discretizes the linear implicit source terms and F(u) is a discretization of the (nonlinear) flux ∂ x F (u), which might include some stabilization terms as well. In the following, we will denote [START_REF] Benjamin | Model Equations for Long Waves in Nonlinear Dispersive Systems[END_REF] as full order model (FOM). The time discretization in ( 6) is done through an implicit-explicit (IMEX) scheme, where some of the terms are discretized with implicit Euler (Φ n+1 ) and others with explicit Euler. In practice, we will will use a IMEX Runge-Kutta (RK) schemes with order matching the spatial discretization orders, but all the discussion can be easily derived from the first order discretization of [START_REF] Benjamin | Model Equations for Long Waves in Nonlinear Dispersive Systems[END_REF]. Hence, we will keep discussing it, and only for the specific problems we will go in detail with the IMEX RK discretization. To have stability of the hyperbolic operator, we need to impose some CFL conditions on the timestep of the type ∆tρ(JF (u)) ≤ CFL∆x, where ρ(JF (u)) is the spectral radius of the Jacobian of F in u and CFL is a constant smaller than 1.

Typically, the most computationally intensive part is the solution of the linear systems, derived from the elliptic operator. While in one dimensional cases the linear systems are tridiagonal for simple discretizations and very efficient algorithms like Thomas algorithm can render very fast the solution of such systems, this is not the case in general. Even in this very favourable case, we find that the cost of the solution of these linear systems is still the largest of the problems (between 60% and 90%).This is also due to the fact that the computation of the explicit flux and source terms can be efficiently done in parallel. The main goal of this work is to investigate hybrid methods in which we reduce the cost of the linear system, while keeping the evaluation of the fully nonlinear terms.

Reduction of the general model

In this section, we describe how to perform some reduction techniques on the operators presented in the FOM. As described above, there are two types of operations in the previous discretized model. The first ones are the vector based operations constituted of the fluxes, sources and diffusion terms, of derivatives and of sums of all terms and (sparse) matrix vector multiplications. The second one are the matrix based operations which include essentially the solution of linear systems. Even if in some particular cases efficient ad hoc algorithms can be used (e.g. Thomas algorithm for the tridiagonal matrices), we have in mind more challenging applications, e.g. two-dimensional problems or high order methods, where such algorithms cannot be used. In those cases, we need to recast to sparse linear solver methods, which in general allow at best to obtain costs of the order of O(N h log(N h )) at each time step. These systems require greater computational costs to be solved, with respect to the vector part, and are the ones that mostly weigh on the FOM algorithm.

Our objective here is to perform some model reduction to gain control on the computational costs. To this end, we introduce different approximations of the model, via reduced basis expansions [START_REF] Harris | Array programming with NumPy[END_REF]. We want to find an appropriate set of N RB basis of vectors in R N h , with N RB N h , allowing to express the approximate solutions as

u n ≈ V ûn , ûn ∈ R N RB (7)
with the columns of V ∈ R (d×N h )×N RB containing the basis vectors, and ûn the vector of the reduced basis coefficients.

A classical way of obtaining a reduced basis method is by projecting onto the reduced space. We introduce a test matrix W ∈ R (d×N h )×N RB and then we project the whole equation [START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF], i.e.,

W T M Φ -X t V Φn+1 + W T X x V ûn = 0, ( 8a
)
W T M u V ûn+1 -ûn ∆t + W T F(V ûn ) + W T S ex (V ûn ) + W T S im V ûn+1 = W T V Φn+1 . ( 8b
)
This open the question on which matrix W we should consider. If we choose W = V we obtain a Galerkin projection, while different choices leads to Petrov-Galerkin projections. There are several choices that one can take to obtain W , for example controlling the energy of the system [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. We will comment in each model what we use as test matrix W .

It should be noted that with the reduction (8), if N RB N h , one can get rid of the costs of the solution of the large systems of equations. If we define, for every matrix

A ∈ R N h ×N h , its reduced version as  := W T AV ∈ R N RB ×N RB , (9)
we can reduce the computational costs of ( 8), obtaining

MΦ -Xt Φn+1 + Xx ûn = 0, (10a) Mu ûn+1 -ûn ∆t + W T F(V ûn ) + W T S ex (V ûn ) + Ŝim ûn+1 = Î Φn+1 , (10b)
where I ∈ R N h ×N h is the identity matrix. In this model, the costs of the linear systems are independent of N h , while still some computations (the computations of fluxes and source terms) are still dependent of N h . Nevertheless, their implementation can be easily parallelized and their computational costs do not impact as much as the solver of the full linear system. Still, the computational cost reduction that one can obtain in this context is limited by these terms.

Choice of the reduced space

The underlying technique used in this work is proper orthogonal decomposition (POD) [START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF][START_REF] Harris | Array programming with NumPy[END_REF] or principal component analysis (PCA). It consists in a singular value decomposition (SVD) of the snapshot matrix (a simple juxtaposition of the FOM solutions for different parameters and timesteps) of which we retain the most energetic eigenvectors, ordering the singular values. The retained N RB eigenvectors form the basis V ∈ R N h ×N RB for the reduced order model (ROM).

The number of retained modes/eigenvectors can be chosen either directly, or according to a tolerance on the discarded energy, i.e.,

N RB := arg max N N j=1 σ i Nmax j=1 σ i ≤ 1 -tol P OD , (11)
where N max is the minimum between N h and the number of snapshots considered.

The POD (and the linear ROMs) are based on the superposition ansatz, for which the solutions can be reasonably represented by a linear combination of few modes u(x) ≈ N RB j=1 v j (x)û j . This is not always true, in particular for hyperbolic problems. More sophisticated algorithms exist to obtain nonlinear ROMs [START_REF] Burman | The edge stabilization method for finite elements in CFD[END_REF][START_REF] Tonelli | Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF][START_REF] Strazzullo | POD-Galerkin model order reduction for parametrized nonlinear time dependent optimal flow control: an application to Shallow Water Equations[END_REF][START_REF] Peherstorfer | Stability of discrete empirical interpolation and gappy proper orthogonal decomposition with randomized and deterministic sampling points[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF], but these will not be considered in this work.

Hyper-reduction of nonlinear operators

To further alleviate the computational cost one can introduce some linearization of the nonlinear operator. There are several algorithms for this. The approach selected in this work is the so-called empirical interpolation method (EIM) [START_REF] Bacigaluppi | Implementation and Evaluation of Breaking Detection Criteria for a Hybrid Boussinesq Model[END_REF]. The EIM is a very simple approximate-then-project approach. It is perhaps not the most advanced technique to handle nonlinearities, however it will serve as a reference in terms of cost. Using it with different tolerance will allow to have a fair comparison. The hyper-reduced model will be referred to as EIMROM in the following. We want to mention that several other techniques have been developed in recent years and more stable algorithms are available [START_REF] Wei | Generation of waves in Boussinesq models using a source function method[END_REF][START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF][START_REF] Dingemans | Water wave propagation over uneven bottoms: Linear wave propagation[END_REF][START_REF] Cauquis | Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study[END_REF][START_REF] Duruflé | A numerical study of variable depth KdV equations and generalizations of Camassa-Holm-like equations[END_REF]. Still, the focus of this work is not on the hyper-reduction, hence we will stick to a simple EIM algorithm.

Given a nonlinear flux, the goal of the EIM procedure is to approximate it through an interpolation into few magic points {z j } N EIM j=1 , and some basis functions {ψ j } N EIM j=1 ⊂ R N h . We denote by Z ∈ R N EIM the vector of the magic points, and by Ψ ij := ψ j (x i ) ∈ R N h ×N EIM the matrix of the reduced basis functions evaluated in the magic points, i.e., the slice of V in the magic points. Given a nonlinear flux ϕ{u} : (R N h ) d → R N h , we approximate it as

ϕ EIM {u}(x) := N EIM j=1 ψ j (x)ϕ{u}(z j ), (12)
which, after the projection, reads

φEIM (u) = W T ϕ EIM (u) = N EIM j=1 W T ψ j (x)ϕ{u}(z j ) = W T Ψϕ{u}(Z). ( 13
)
The reduction lies in the computation of the flux ϕ{u} only in few magic points {z j } N EIM j=1 , while we can pre-compute and store the matrix W T Ψ ∈ R N RB ×N EIM in the offline phase.

The choice of the points and of the number of the EIM basis functions is crucial in order to maintain a good accuracy of the whole algorithm. Typically, they are chosen in a greedy fashion that minimizes the ∞-norm of the worst approximated fluxes. All the details can be found in [START_REF] Bacigaluppi | Implementation and Evaluation of Breaking Detection Criteria for a Hybrid Boussinesq Model[END_REF][START_REF] Drmac | A new selection operator for the discrete empirical interpolation method-improved a priori error bound and extensions[END_REF]. In practice we will stop such algorithm when the error gets lower than a certain tolerance set proportionally with the error of the RB approximation.

Remark 1 (Implementation details). To assemble and solve the systems at each stage we will use a Python implementation based on numpy [START_REF] Grimberg | On the stability of projection-based model order reduction for convection-dominated laminar and turbulent flows[END_REF] and on Numba [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF]. These packages allow to perform a multithread computation of the vector assembly (all the nonlinear terms), while the solution of the sparse FOM systems will be performed with Thomas algorithm, since the matrices are tridiagonal. This leads to have the largest cost in the solution of the system O(N h ), while the nonlinear operators can be assembled at a lower cost, depending on the architecture of the machine. We will compare also the sparse solver of scipy to emulate the case of a general sparse solver without the a priori knowledge of the tridiagonality. In that version of the algorithm, the whole computational time, except the fluxes, will be affected by a non precompiled code. For ROM systems, which are not sparse, we will use the solve function of the linear algebra package of numpy. The computations all along this work will be performed on an 8 cores Intel(R) Core(TM) i7-10875H CPU @ 2.30GHz.

η(x, t) h(x, t) b(x) h 0 L a 0 b 0 h(x) (a) Dimensional variables a 0 η(x, t) h 0 h(x, t) b 0 b(x) h 0 L a 0 b 0 h 0 h(x) (b) Dimensionless variables

BBM-KdV model

As a first exercise, we study the BBM-KdV equations, which provide a scalar one way approximation of the water wave equations [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Chen | An EIM-degradation free reduced basis method via over collocation and residual hyper reduction-based error estimation[END_REF][START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]. It is useful to write the scalar partial differential equation in terms of the dimensionless parameters

µ = h 0 L , = a 0 h 0 , β = b 0 h 0 ,
with h 0 a characteristic water depth, L the wave length, a 0 the wave amplitude, and with b 0 a characteristic bathymetry amplitude, see Figure 1 for both dimensional and non-dimensional notations. These parameters are classically used to define the different propagation regimes. More precisely, the dispersive character of the waves is measured by µ, while ε is a measure of their nonlinearity. For ε = O(µ 2 ), the BBM-KdV equation provides a 1-way O(µ 4 + µ 2 ) approximation of the water equations. Its classical form for constant bathymetry reads

1 -α p µ 2 ∂ xx ∂ t η + ∂ x η + 3 2 ε η 2 2 + µ 2 p∂ xxx η = 0, p ≤ 1 6 , (14)
having set α p = 1/6 -p. For different values of the parameter p we obtain different approximations. In particular, the KdV equation is obtained for p = 1 6 , while for p = 0 the equation reduces to the BBM model.

For the purpose of its numerical approximation, we manipulate the model equation following [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF], and recast it as

1 -α p µ 2 ∂ xx (∂ t η + ∂ x F (η)) + µ 2 6 ∂ xxx η = -εµ 2 α p ∂ xxx 3 2 η 2 2 , (15)
having introduced the short notation for nonlinear flux

F (η) := η + 3 2 ε η 2 2 . ( 16
)
Within the same approximation hypotheses, the BBM-KdV equation is equivalent to

1 -α p µ 2 ∂ xx (∂ t η + ∂ x F (η)) + µ 2 6 ∂ xxx η = 0. ( 17
)
which we refer to as the modified BBM-KdV equation. Note that both for ( 15) and ( 17), the limit µ → 0 corresponds to the hyperbolic conservation law

∂ t η + ∂ x F (η) = 0. (18)
Equation ( 17) can be recast into the form of (2) by setting u := η, X x := µ 2 6 ∂ xxx , X t := α p µ 2 ∂ xx and S(u) := 0.

Remark 2 (Dimensional equations). One can easily recover the dimensional form of the equation, and more importantly of the variables involved, from the scaling:

x = Lx, t = L c 0 t, η = h 0 η, b = βh 0 b, h = h 0 h, h = h 0 h (19)
with c 0 = √ gh 0 .

Energy conservation

A property of ( 14) which we will use in the following is the existence of a certain number of additional derived conservation laws (see also [1? ]). In particular, ( 14) is endowed with an energy-energy flux pair also verifying a conservation law. This can be shown quite classically by premultiplying by η and manipulating the resulting expression. Setting Ξ(u) = η 2 /2 + εη 3 /2, this leads to

∂ t E + ∂ x F = 0, E(η) := η 2 2 + µ 2 α p (∂ x η) 2 2 , F(η) := Ξ(η) -µ 2 α p η∂ xt η + µ 2 pη∂ xx η -µ 2 p (∂ x η) 2 2 . (20)
From this relation, we know that the energy E(η) is conserved for [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF]. For the simplified model ( 17), this balance is only valid within the O( µ 2 , µ 4 ) asymptotic error of the model, namely

∂ t E + ∂ x F = O( µ 2 , µ 4 ). ( 21
)
Other balance laws may be derived in a similar spirit, see e.g. [START_REF] Ali | Mechanical Balance Laws for Boussinesq Models of Surface Water Waves[END_REF]. In both cases, it is convenient to consider the scalar product associated to E

η, ρ E := 1 2 Ω ηρ + µ 2 α p ∂ x η∂ x ρ (22)
which is essentially a modified H 1 norm scalar product.

Including the bathymetric effects

Generalizations of the KdV equation and of other dispersive one-wave models are discussed in some detail in [START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF][START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF][START_REF] Israwi | Variable depth kdv equations and generalizations to more nonlinear regimes[END_REF]. For our purposes, these generalizations are useful as they allow to broaden the spectrum of benchmarks. We use here the O(εµ 2 , µ 2 β, µ 4 ) approximation obtained from the variable depth model of [START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF], still referred to here as to the modified BBM-KdV model, and reading

(1 -α∂ xx )(∂ t η + (γ + δη)∂ x η) + ω∂ xxx η + νη = 0, ( 23
)
where the coefficient depend now on the bathymetry profile b(x) as (in dimensionless form)

c(x) = 1 -βb(x), α(x) =α p µ 2 , ω(x) = 1 6 µ 2 c(x) 5 , γ(x) = c(x), δ(x) = 3 2 ε, ν(x) = 3 2 c (x). ( 24
)
The equation can be rewritten up to an O(βµ 2 ) as

(1 -α∂ xx )(∂ t η + (γ + δη)∂ x η + νη) + ω∂ xxx η = 0, (25)
to better match [START_REF]On the Formulation of Mass, Momentum and Energy Conservation in the KdV Equation[END_REF]. We can indeed define S(u) = νη to the previous formulation.

As before, a conservation law for the energy [START_REF] Harris | Array programming with NumPy[END_REF] can only be shown within an O(β, µ 2 , µ 4 ).

Full order model discretization

We discretize (23) using a finite difference (FD) approach. The resulting discrete model will be denoted by full order model (FOM). Classically, we consider a discretization of the spatial domain Ω = ∪ N h i=1 [x i , x i+1 ] into N h equi-spaced intervals, and define ∆x :

= x i+1 -x i . The temporal domain [0, T ] is subdivided into K time steps [0, T ] = ∪ K k=1 [t k-1 , t k ],
. We denote by u k i the numerical approximation of a quantity u in a point x i at time t k , and, with an abuse of notation, we use the same symbol for the continuous quantity u : R + × R → R, and for the discrete one: u ∈ R N h ×K . For convenience, we use the following notation for the classical centered finite difference operators D, D 2 and D 3 :

R N h → R N h : (Du) i = u i+1 -u i-1 2∆x , (D 2 u) i = u i+1 -2u i + u i-1 ∆x 2 , (D 3 u) i = u i+2 -2u i+1 + 2u i-1 -u i-2 2∆x 3 . (26)
We start by defining the IMEX first order Euler discrete approximation

(I -αD 2 )Φ n+1 + ωD 3 η n = 0, η n+1 -η n ∆t + (γ + δη n )Dη n + νη n + J (η n , λ) = Φ n+1 . (27)
Here, the vector-by-vector multiplication is meant component-wise in η n Dη n . The term J (η n , λ) is a high order numerical diffusion term, which depends also on the spectral radius λ, where λ i = γ + δ|η i |. In this work we have chosen to use to this purpose an operator originated from the continuous interior penalty (CIP) approach [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF] which can be written as

J (v, λ) j := d∆x 3 {λ j+1 (D 2 v) j+1 -2λ j (D 2 v) j + λ j-1 (D 2 v) j-1 } (28)
where d ∈ R + is a stabilization parameter, which in practice we have set to 1 in all numerical simulations. Also [START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF] can be reshaped into (6) simply defining

M Φ = M u := I, X t := αD 2 , X x := ωD 3 , F(η) := (γ + δη)Dη + J (η, λ) and S ex (η) := νη.
Starting from [START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF], we build a second order SSPRK(2,2) IMEX approximation. The time-step is computed from the CFL condition: ∆t ≤ CFL∆x/(max i λ n i ), with CFL < 1.

Benchmarks for the modified BBM-KdV equation

In this section we describe some problems that we will use to test the reduced algorithms.

Propagation of a periodic monochromatic wave

The first one starts with a long wave on the whole domain, i.e., a cosine on a periodic domain. This cosine deforms and splits into different waves that travel at different speeds. We use in [START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF] with the dimensionalized coefficients passed through [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]: h 0 = 1, g = 9.81, a 0 = 0.04, u 0 = a 0 cos(x/10) on the domain [0, 20π] and final time is T = 200. The FOM will use N h = 2000 degrees of freedom and time step will be set with a CFL number of 0.2. In Figure 2(a) the initial and the final solution are depicted.

Undular bore propagation

This test emulates a dam break which develops into several waves. The initial condition is 5) and the parameters in [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF]

u 0 = a 0 1 - 1 1 + e -4(x-
(x) = 0.07          x-45 10 if 45 < x ≤ 55, 1 if 55 < x ≤ 75, -x-80 5 if 75 < x ≤ 80, 0 else.
The initial condition is defined by u 0 = a 0 cosh -1 ( x-22 5 ), final time is T = 60, the parameters are h 0 = 1, g = 9.81, a 0 = 0.1, CFL = 0.1 and the domain is [0, 100]. In Figures 2(c) and 2(d) the evolution of the solution in time is shown. Now, we compare the computational times of the solution of the linear systems for the three benchmark problems with respect to other costs in the problem solution. In Table 1 we compare the costs of the linear systems with respect to the whole simulation, while in Table 2 we see the ratio of the cost of the solution of the linear systems over the cost of the nonlinear fluxes. Let us focus on Thomas algorithm first, we notice that boundary conditions make a huge difference. Indeed, periodic boundary conditions imply the solution of Thomas algorithm three times instead of just one and this increases the costs of the system solver.

On the other side, using a general solver for sparse system, like the one provided by scipy, increases the computational costs of the solver even by a factor of around 10. In particular, the distinction between different boundary conditions is less evident in this case. This algorithm is more interesting in perspective of a more general discretization or in case of multidimensional problem, but less efficient.

Projection based reduction

As said before, the computational operation in (27) that we mainly aim at reducing is the solution of the linear system, while the flux and source terms are easily parallelized.

As before, we introduce the reduced variable of dimension N RB N h so that we approximate

η n ≈ V ηn , ηn ∈ R N RB (29)
with the columns of V ∈ R N h ×N RB containing the basis vectors, and ηn the vector of the reduced basis coefficients.

Two projection based model order reduction (MOR) strategies exist [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF], equivalent under certain conditions. The first one is the projection with a test matrix W ∈ R N h ×N RB and it is obtained by projecting the whole equation ( 27), i.e.,

W T V ηn+1 -r(V ηn , Φ n (V ηn )) = 0, (30a) r(η, Φ) = η -∆t ((γ + δη)Dη + J (η, λ) -Φ) . (30b)
This opens the question on which matrix W we should consider. If we choose W = V we obtain a Galerkin projection, while different choices lead to Petrov-Galerkin projections.

An alternative approach, allowing to answer this question, is to construct the reduced space using a minimization problem. Using in particular the scalar product ( 22), we can for example define ηn+1 = arg min

x V x -r(V ηn , Φ n+1 (V ηn ), V x -r(V ηn , Φ n (V ηn ) E . (31)
Using the discrete form of ( 22) (in dimensional variables), we can introduce the scalar product matrix

Θ := ∆x(I + αD T D),
where I ∈ R N h ×N h is the identity matrix and D ∈ R N h ×N h is the (first) derivative matrix, and such that at the discrete level u, v E = u T Θv.

Then, we can write the minimization [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF] as the solution of

2V T ΘV ηn+1 -2V T Θr(V ηn , Φ n+1 (V ηn )) = 0, (32a)
suggesting that minimizing the energy in the projected space requires testing with W = Θ T V in (30a).

Approximation accuracy of the problems

For the benchmarks involving the modified BBM-KdV equation we have applied the POD on 1000 snapshots in time, with fixed parameters. The FOM dimension is here N h ≈ 10 3 , which is a representative order of magnitude. In practice, we have used 2000 points. Concerning the potential for reduction, we need to look at the error decay in terms of singular values. For the tests considered here, the decay has been plotted in Figures 3(a) to 3(c). As we can see, the energy drop is not exceedingly fast, still being exponential from the very beginning. In particular, for compact travelling waves as in the last benchmark, the advection dominated regime slows down the convergence of the eigenvalues. Nevertheless, if we want to keep the error of the order of 10 -2 or 10 -3 the number of eigenmodes needed is in between 50 and 70, which is much below the typical mesh size used in the FOM and can already reduce the computational cost of the linear system (even if full). We remark that the computational reduction is not proportional to the dimension as the original FOM has a sparse linear system to be solved, while the ROM has a full linear system.

Reduction of the linear operator

We consider here the reduction of the linear system which gives Φ in [START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF]. In general this operation can scale as an O(N h 2 ), unless particularly optimized algorithms are implemented. Indeed, in our FOM simulations we will use Thomas algorithm which only scales as O(N h ). When the linear operator is the only one reduced, we still need to assemble the discretization of the nonlinear part of the PDE. For this reason we will refer to this approach as to a partial discrete reduced order model or pdROM.

The pdROM discrete evolution equation can be written as Also here, we can see the connection with [START_REF] Burman | The edge stabilization method for finite elements in CFD[END_REF], where all the terms can be matched. In practice, for N RB sufficiently small, the N RB × N RB matrix A N RB := W T (I -αD 2 ) -1 D 3 V can be pre-computed and stored. Denoting with M N RB := W T V ∈ R N RB ×N RB , the pdROM update can be effectively implemented as

W T V ηn+1 = W T V ηn -W T ∆t ((γ + δV ηn )DV ηn + J (V ηn , λ)) + ∆tW T Φ n+1 , (33a)
W T Φ n+1 = -ωW T (I -αD 2 ) -1 D 3 V (η n + ν ηn ). ( 33b
Φn+1 = -ωA -1 N RB (η n + ν ηn ), (34a) F(η) = W T ((γ + δη)Dη + J (η, λ)) , (34b) ηk+1
= ηn + ∆tM -1 N RB -F(V ηn ) + Φn+1 . (34c)
It is clear that, with this approach, the parameters γ, δ, ω can be easily modified in the online computations. Also α can be decoupled from the reduction process and used to assemble the linear reduced system, but in this context α represent the topography parameters and we do not aim to change it. In the system case we will see how to modify a parameter in the LHS matrix.

Remark 3 (Modification to nondispersive code). It is very important to stress that the computational results, as well as all conclusions concerning cost reduction can be obtained keeping a classical solver for nondispersive nonlinear equations, to which we add the reduced dispersive term, i.e.,

η k+1 = η n -∆t ((γ + δη n )Dη n + νη n + J (η n , λ)) + ∆tβV Φn+1 , (35a) Φn+1 := -(W T (I -αD 2 )V ) -1 W T D 3 η n , (35b)
where

B N RB := W T (I -αD 2 )V ∈ R N RB ×N RB
can be stored and inverted at low computational cost.

Nevertheless, working on two different spaces contemporary does not guarantee the same accuracy as in the previously presented algorithm. We have observed that reducing and projecting only the Φ term leads to much more inaccurate methods. We believe that this is due to spurious modes that the reconstruction of Φ generates and that are not appropriately damped in the full dimensional space, while they are in the pdROM approach. In the test section we will show that both algorithms have comparable computational costs, but (34) provides more stability.

Remark 4 (EIMROM stability). As shown in Section 2.2.2, we can reduce the computational costs of the source and flux terms with a hyper-reduction technique. The extra reduction of the EIM algorithm does not guarantee anymore the energy minimization [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF]. This is noticeable when the dimension of the EIM space N EIM is smaller or close to the dimension of the RB space N RB . In that regime, it is easy to see instabilities and Runge phenomena. In the simulation section, we will observe in which regime we should stay with the dimension N EIM . Other works already shown that the dimension of the EIM space should be larger that the reduced space one to obtain stability of the reduced model [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF][START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF][START_REF]On the Formulation of Mass, Momentum and Energy Conservation in the KdV Equation[END_REF][START_REF] Cauquis | Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study[END_REF]. Some novel algorithms have been proposed in recent papers [START_REF] Cauquis | Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study[END_REF][START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF] based on different techniques as the over-collocation method and variants of the Gappy-POD. We will not investigate in deep this behavior in this work.

Remark 5 (Fully implicit). We want to remark that even for the pdROM simulations, even if the solution is energy stable up to an O(εµ 2 , µ 4 ), oscillations may arise due to the under-resolution in the reduced space. This would not lead to blow ups of the solution, but it might have some unphysical behaviors. In the case of the EIM algorithm the oscillations might also lead to instabilities and blow ups when the EIM space does not properly approximate the nonlinear flux manifold. In order to alleviate this issue, it might be helpful to consider a fully implicit discretization or, even better, an implicit Euler time discretization. On the other side, this would lead to extra computational costs due to nonlinear solvers that might heavily slow down the reduced formulation, both in the pdROM and EIMROM contexts. In particular, the ROM computational costs with respect of the FOM ones, where the benchmark algorithms are explicit for the hyperbolic equation, might be disadvantageous, at least for one dimensional simulations.

Simulations for KdV-BBM

Modus operandi

In this section, we will test the pdROM and EIMROM in different regimes, showing their error and the computational costs with respect to the FOM solutions.

1. We will proceed studying a time-only manifold of solutions, in order to understand which errors our algorithms can reach with different dimensions of reduced basis spaces. In this first phase, we will fix the parameters of the problem, we will compute the FOM solution on a mesh with N h = 2000 points and we will collect 1000 snapshots in time u(t k ) and the respective nonlinear fluxes F (u(t k )). These snapshots will form the training set used to choose the reduced basis space and the EIM space, setting either the dimension of the space or the tolerance. Then, different simulations of the original problem are run with pdROM and EIMROM, recording their errors and their computational times.

We will picture also one significant simulation of the reduced algorithms for a visual comparison for each test. 2. The second stage will consider a solution manifold generated by time and parameters. In practice, we will vary h 0 , on which all the parameters depend in the BBM-KdV problem, see [START_REF] Karczewska | Can simple KdV-type equations be derived for shallow water problem with bottom bathymetry?[END_REF]. In particular, we build the training sample using 10 randomly sampled h 0 ∈ [h min , h max ] = [0.7 h0 , 1.3 h0 ], and for each of these parameters, we collect 1000 snapshots at uniform timesteps. We populate the ROM and EIM spaces using these 10.000 snapshots and their fluxes. Then we test the reduced algorithms for h0 for a larger final time (for every test we will specify which one), checking again how the reduced procedures perform with different tolerances for EIM and POD. Finally, we test the methods for a parameter h 0 = 0.63 • h0 outside the training set and, when sensible, with a different initial condition (for every test we will specify which one), to assess the property of the reduced algorithms in the extrapolation regime.

Propagation of a periodic monochromatic wave

For this test, we use h0 = 1 and final time T = 200. This test is fairly simple and it can be well compressed with few basis functions. Nevertheless, it already shows what can be achieved with and without EIM. We collect samples for different times of one simulation obtaining a time-only reduction. In Figure 4(a) the decay of the relative error is exponential with respect to the dimension N RB of the reduced space (different points correspond to N RB = 10k for k = 1, . . . , 10) and with few basis functions we can reach a very small error. On the other side, this error could be worsen or even the simulation might not run until the final time if we add the EIM algorithm with not enough basis functions. The error for those methods stagnates around the tolerance of EIM and if N EIM N RB the method becomes unstable.

The computational time for the KdV-BBM pdROM compared to the FOM using Thomas method is between 20% and 30% in Figure 4(d). The hyper-reduction allows to further reduce the in between 4% and 10% of the original computational time. In Figure 4(g) we see the simulations run with pdROM with different dimensions of the reduced space. Already with N RB = 20 we obtain an accurate solution. In Figure 4(j) we see that a minimum amount of EIM basis is necessary not to have unstable solutions and to be accurate.

For the second phase, we introduce the reduction in the parameter space. For Figures 4(b), 4(e), 4(h) and 4(k), we test the obtained reduced spaces and the respective algorithms onto the parameter h0 = 1 and for final time T = 250. The error decay in Figure 4(b) is similar to the one of the time-only reduced case. The largest differences stays in the differences of the number of basis functions needed to reach a certain tolerance. In the example in Figure 4(h), with N RB = 27 we reach a tolerance of 0.3. With timeonly reduction N RB = 10 are enough to reach the same tolerance, though not reaching the same accuracy in the simulation see Figure 4(a). The same holds for the EIM algorithm. This is noticeable in the computational times, that increase for all algorithms given certain thresholds, see Figures 4(d) and 4(e), this is a measure of the increased difficulty of the problem. In Figure 4(h) a representative simulation is shown, where we see that with N RB = 27 basis function we get an already accurate approximation with pdROM, while N EIM = 72 is still not enough to obtain a good accuracy for EIMROM. Increasing N EIM = 119 in Figure 4(k), we can obtain an accurate solution also for EIMROM. This is accordance with the observations made in [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF][START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF][START_REF]On the Formulation of Mass, Momentum and Energy Conservation in the KdV Equation[END_REF][START_REF] Cauquis | Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study[END_REF]. In this case, the computational time of EIMROM is 3 times faster than the pdROM, and decreasing tol EIM we can achieve even smaller errors. A particular attention must be observed in choosing the tolerance of EIM small enough with respect to tol P OD , as one can see in Figure 4(b) that the EIMROM simulations can explode if tol EIM is not small enough.

In Figures 4(c), 4(f), 4(i) and 4(l) we test the algorithms with h 0 = 0.63 / ∈ [h min , h max ], a 0 = 0.05 instead of 0.04 and initial conditions u 0 (x) = -1.2a 0 cos( 410 π(x -2)) + a0 10 cos( 8 10 πx). Similar conclusion can be drawn also in this case. The main difference is a larger error, due to the higher wave amplitude of the solution.

In Figure 5 we can see the relative error of pdROM N RB = 34 and EIMROM with N EIM = 93 varying a 0 and h 0 in the initial conditions and the problem settings. We plot the error as a function of ε and µ, in order to understand how the nonlinearity and the dispersion affect the solutions. We know that as h 0 gets smaller or a 0 gets bigger (ε larger), the nonlinearity of the problem becomes more pronounced, hence leading to more oscillations and steeper gradients. In the strongly nonlinear regime (ε = a 0 /h 0 10 -1 ) the hypothesis of having small ε are not valid anymore and the error increases rapidly. On the other side, the dispersion coefficient µ = h 2 0 /L 2 does not affect much the error for a fixed ε, in Figure 5. Overall we can see that, as long as nonlinearity does not become predominant, the pdROM provides lower errors over the parameter space showing potential for a more robust approximation outside the training set.

Undular bore propagation

This is a very complex test for model reduction. We start by choosing h0 = 1, a 0 = 0.04 and final time T = 20. In order to take into account the non homogeneous Dirichlet boundary condition at the left boundary, we subtract from the method the residual equation computed on a lift solution u lif t that verify the Dirichlet boundary conditions. In this way, the resulting scheme can be applied with homogeneous Dirichlet boundary condition on the left.

This tests is way more challenging than the previous one, in particular because flat zones and oscillating ones coexist and move around the domain. In particular, the EIM method easily fails in capturing this behavior without falling into Gibbs phenomena in the flat area. Indeed, in Figure 6(a) we see the error decay for only-time reduction both for EIMROM and pdROM and the former become very oscillatory when too many POD bases are used. The pdROM uses between 10% and 20% of the FOM computational time, while the EIMROM, in the few situations where it does not explode, uses around 5% of the FOM time, see Figure 6(d). In Figure 6(g) we see for different POD basis functions how we approach the exact solution, while in Figure 6(j) we fix the number of POD basis functions very low and we change the EIM ones. In this stable situation as the EIM bases increase we slowly converge towards the pdROM solution, even if the right part of the solution struggles with obtaining the right behavior.

When exploring the parameter domain in the training set h 0 ∈ [0.7, 1.3], we observe a similar situation. The decay of the error on a parameter inside the training domain computed at time T = 15 for pdROM decays exponentially in Figure 6(b), while EIMROM is affected by oscillations and struggles with obtaining reasonable results. The computational time of the pdROM is between the 15% and the 25% of the FOM computational time as shown in Figure 6(e), while the EIMROM, when the simulation is reasonable for few basis functions, needs only around the 5% of the FOM computational time. In Figure 6(h) we see for a tolerance of the POD of 3 • 10 -2 we need N RB = 33 and the pdROM is quite close to the FOM solution, while EIMROM with same tolerance and N EIM = 52 start oscillating around the original position of the dam. In Figure 6(k) the situation worsen for the EIMROM which becomes very oscillatory for N RB = 62 and N EIM = 77 and this explains why the error is so large.

Taking the parameter out of the training set, considering h 0 = 0.63 and a 0 = 0.03, the EIMROM becomes even more oscillatory and its error often explodes, but also the pdROM suffer of the larger flat zone in the part of the domain, leading to instabilities. The error decay is indeed affected by this and when too many POD bases are considered the error grows, see Figure 6(c) and the simulation in Figure 6(l). Computational times stays similar to the previous case Figure 6(f), and, if properly choosing the number of basis functions, both pdROM and EIMROM gives decent solutions, though the EIMROM oscillates around the flat area, see Figure 6(i).

Clearly this test case is pushing over the limit of the method as the advection dominated part is leading the computations and the contrast between flat and oscillating region leads very easily to Gibbs In Figure 7 we see how the error of the pdROM with N RB = 46 varies with respect to a 0 and h 0 . We can see a strong gradient when h 0 and a 0 increase in the area of unstable solutions. Again, the main factor here is the relevance of nonlinearity. As long as the waves remain dispersive and we are far from very shallow bores, we observe an error plateau which seems to indicate a certain robustness with respect to the variation of parameters. It must be remarked however, that outside this area, the error grows relatively fast also for pdROM. Some improved treatment of nonlinearity is necessary to handle this kind of problems when getting close the the shallow/highly nonlinear range.

Solitary waves interacting with a submerged bar

For this test, we use h0 = 1, final time T = 60 and a 0 = 0.1. This test has similar results to the Propagation of a periodic monochromatic wave. Its advection dominated character is more pronounced than in the Propagation of a periodic monochromatic wave, but the error decays quickly enough as we can see in Figure 8(a). The pdROM solution with N RB = 40 basis function is indistinguishable from the FOM one in Figure 8(g) and, adding the EIM approximation, we can see different results. In Figure 8(j) for the EIMROM simulation with N RB = 42 and N EIM = 65 we still do not have enough accuracy and the solution oscillates, with N EIM = 75 the solution is much more stable, but we can still see differences in the peak of the wave, finally for N EIM = 118 observe a very accurate solution. Similarly to before, the computational costs scale between 15% and 30% of FOM for pdROM and between 5% and 25% for EIMROM, see Figure 8(d).

In the second phase, we consider also here 10 simulations for randomly chosen h 0 ∈ [0.7 h0 , 1.3 h0 ] and we use them to compute the reduced space with the POD. The error decay with respect to the POD tolerance in Figure 8(b) is similar to the first phase, but the dimension of the spaces are larger as we notice in the computational time plot in Figure 8(e) up to a 40% of the FOM computational time for pdROM. In Figures 8(h) and 8(k) we see accurate pdROM simulations with N RB = 72 for a POD tolerance of 0.01 and a computational time of around 25% of the FOM time, while we need N EIM = 259 to get very accurate EIMROM solutions for computational costs around 20% of the FOM ones. Again, the EIM space should be very large in order to obtain accurate and stable results [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF][START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF][START_REF]On the Formulation of Mass, Momentum and Energy Conservation in the KdV Equation[END_REF][START_REF] Cauquis | Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study[END_REF].

When testing on a parameter outside the training domain, i.e., a 0 = 0.15 and h 0 = 0.63, we have that the exact solution is more oscillatory, see Figures 8(i) and 8(l), and to catch properly the solution we need more basis functions in both spaces. Indeed, the errors are larger in Figure 8(c), while the computational times stays similar.

In Figure 9 we see how the error varies for pdROM with N RB = 72 and EIMROM with N EIM = 227 for different parameters h 0 and a 0 , plotting the error with respect to ε and µ. As before we see that nonlinearity really plays a major role, with the shallow water/hyperbolic limit being poorly handled by both approaches. Away from this limit (left part of the figures) the EIMROM error is larger than the pdROM one.

We remark that the computational cost of the linear system is the largest one in the FOM and, while doing a pdROM, we obtain already a great reduction. In Appendix Appendix B we compare the pdROM costs with the reduction of only the linear system for Φ (evolving the full FOM) as suggested in Remark 3. Essentially, the cost of the two methods are identical. Moreover, changing the linear solver into something more generic results in an even stronger reduction of the computational costs. We compare the costs of the Thomas algorithm and of a generic sparse solver in Appendix Appendix B.

An extension to a Boussinesq system

We propose in this section a possible extension to well known dispersive shallow water model by Madsen and Sørensen [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF]. In the dimensional form the model equations read [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF]

∂ t η + ∂ x q = 0, ∂ t q -T t [∂ t q] + ∂ x (qu) + gh∂ x η -T x [η] = 0, with T t [•] := B 1 h2 ∂ xx [•] + 1 3 h∂ x h∂ x [•] and T x [•] := gB 2 h2 2∂ x h + h∂ x ∂ xx [•], (36)
where, with the notation of Figure 1,η(x,t) is the water level variation from the reference value, q(x, t) is the discharge, b 0 b(x) is the bathymetry profile, h 0 is the reference value for water height, h(x

) = h 0 -b 0 b(x)
is the water height at rest, h(x, t) = η(x, t) + h(x) is the water height, u(x, t) is the depth averaged speed and it is linked to q by q(x, t) = h(x, t)u(x, t). System (36) provides a O(εµ 2 , µ 4 ) approximation of the full nonlinear wave equations. It contains the nonlinear terms belonging also to shallow water equations, and they are the only nonlinear terms of the model. Similarly to Boussinesq equation and to the previous scalar model there are other linear dispersion terms which scale with constants B 1 and B 2 . The model constants B 1 and B 2 are set to B 2 = 1/15 and B 1 = B 2 + 1/3 to optimize the linear characteristics of the model (phase and shoaling) with respect to the full Euler equations [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF].

As before, we manipulate the equations for the purpose of their numerical solution. First we recast the model in dimensionless form:

∂ t η + ∂ x q = 0, ∂ t q -µ 2 T t [∂ t q] + ε∂ x (qu) + h∂ x η -µ 2 T x [η] = 0, with T t [•] := B 1 h2 ∂ xx [•] + 1 3 h∂ x h∂ x [•] and T x [•] := B 2 h2 2∂ x h + h∂ x ∂ xx [•]. ( 37
)
We then isolate the nonlinear hyperbolic operator and write the system as

∂ t η + ∂ x q = 0, (1 -µ 2 T t )[∂ t q + ε∂ x (qu) + h∂ x η] + µ 2 T t [ε∂ x (qu) + h∂ x η] -µ 2 T x [η] = 0. ( 38
)
In this case we do not neglect any of the small order terms to keep the original model, often used in literature (see e.g. [START_REF] Benjamin | Model Equations for Long Waves in Nonlinear Dispersive Systems[END_REF][START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF][START_REF] Argaud | Stabilization of (G) EIM in presence of measurement noise: application to nuclear reactor physics[END_REF]) and references therein. The final model can be written as a system of 3 equations, given by

     ∂ t η + ∂ x q = 0, (1 -µ 2 T t )[ψ] = µ 2 {T x [η] -T t [ε∂ x (qu) + h∂ x η]} , ∂ t q + ε∂ x (qu) + h∂ x η = ψ. (39)
Also for this model, we show the analogy with the model defined in [START_REF]On the Formulation of Mass, Momentum and Energy Conservation in the KdV Equation[END_REF]. Indeed, if we consider

u := (η, q) T , ∂ x F (u) := (∂ x q, ε∂ x (qu) + h∂ x η) T , S(u) = 0, X x (u) := 0, -µ 2 T x [η] -T t [ε∂ x (qu) + h∂ x η] T , X t (u) := 0, µ 2 T t T ,
one obtains [START_REF] Ricchiuto | Upwind residual discretization of enhanced Boussinesq equations for wave propagation over complex bathymetries[END_REF]. For this model, one cannot show an exact conservation for a total energy. For other dispersive nonlinear systems one can show the conservation of an energy of the form E = gh 2 2 + hu 2 2 + h 3 (ux) 2 6 [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF]. Differently from the BBM-KdV equation, the relation between this energy and the main unknowns is nonquadratic, and it is not straightforward to use it to build the reduced basis and the test reduced space as the one studied before. Based on regularity requirements, we will instead investigate a formal extension of the method proposed in the scalar case.

Wave generation and boundary conditions

We will consider two families of solutions: solitary waves, monochromatic periodic waves. Exact solitary waves for (36) can be obtained assuming η = η(x -Ct), and q = Cη, which allows to obtain an ODE from the second equation in [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF]. The solution is uniquely defined given amplitude and depth at rest (a 0 , h 0 ). In particular, setting C 2 0 = gh 0 , the celerity C verifies the consistency condition

C 2 C 2 0 = a 2 0 2h 2 0 1 + a0 3h0 a0 h0 -log(1 + a0 h0 ) (40)
Details can be found in [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF][START_REF] Ohlberger | Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing[END_REF]. The ODE is integrated here with a built-in solver in scipy.

To generate monochromatic waves, we have included an internal wave generator. Different definitions exist in literature [START_REF] Taddei | A registration method for model order reduction: data compression and geometry reduction[END_REF][START_REF] Mojgani | Lagrangian basis method for dimensionality reduction of convection dominated nonlinear flows[END_REF][START_REF] Torlo | Model reduction for advection dominated hyperbolic problems in an ALE framework: Offline and online phases[END_REF][START_REF] Walkley | A numerical method for extended Boussinesq shallow-water wave equations[END_REF]. Following [START_REF] Torlo | Model reduction for advection dominated hyperbolic problems in an ALE framework: Offline and online phases[END_REF][START_REF] Walkley | A numerical method for extended Boussinesq shallow-water wave equations[END_REF], we add to the η equation a compact forcing term with periodic variation in time:

∂ t (η + h iwg ) + ∂ x q = 0, (41)
where the internal generator is defined by

h iwg (x, t) := A iwg f iwg (x) sin(ωt), (42)
with ω = 2π/T the frequency of the signal and T the period, f iwg the spatial damping function

f iwg (x) = Γ iwg e -(x-xiwg) 2 /d 2 iwg . (43)
We refer to [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF][START_REF] Walkley | A numerical method for extended Boussinesq shallow-water wave equations[END_REF] for the tuning of the parameters and we will use the following ones

Γ iwg = 0.185 g/h 0 T, d 2 iwg = gα 2 iwg h 0 T 2 80 , T = 2.525, α iwg = 4, (44)
and x iwg the center of the generator and α iwg ∈ [0, 4] is a parameter controlling the amplitude of the wave and β iwg is just a rescaling factor. Also A iwg is a rescaling factor that will be changed along the simulations to test the problem for different parameters.

Finally, outflow conditions are imposed through sponge layers. Following [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF] these terms are written as diffusion operators which will be discretized implicitly and added to the mass matrix. At the continuous level they have the form

S(v) = -ν∂ xx v, ν(x) :=    n 1 1-e n 2 x-X s1 X s2 -X s1 1-e , X s1 ≤ x ≤ X s2 , 0, else. (45)
Here n 1 and n 2 are two coefficients that should be tuned with the problem, we will use the values n 1 = 10 -3 and n 2 = 10 as suggested in [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF]. The points X s1 and X s2 are fixed close to the boundary and X s2 will coincide with the boundary itself. Also here, both these linear source terms can be added in the formulation to match (2) defining S im accordingly.

Full order discrete model

The approach used to discretize [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF] is similar to the one used in the scalar case. We combine a linear continuous Galerkin method (see [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF] and references therein) with an SSPRK(2,2) time integrator [START_REF] Ghavamian | POD-DEIM model order reduction for strain-softening viscoplasticity[END_REF]. To introduce a little dissipation we again make use of a continuous interior penalty operator (CIP) [START_REF] Brocchini | A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics[END_REF]. The resulting discrete equations can be written as

M (η + h iwg ) (s) - s-1 r=0 ρ s r (η + h iwg ) (r) + ∆tSη (s) = -∆t s-1 r=0 θ s r N η (η (r) , q (r) ) + J (η (r) , λ (r) ) (46a) (M -T t )ψ = - s-1 r=0 θ s r T t M -1 N q (η (r) , q (r) ) + T x η (r) (46b) M q (s) - s-1 r=0 ρ s r q (r) + ∆tSq (s) = ∆tMψ -∆t s-1 r=0 θ s r N q (η (r) , q (r) ) + J (q (r) , λ (r) ) (46c)
where the index in bracket notation is to denote the Runge-Kutta i.e., {η (s) , q (s) } s=1,2 , and ρ s r , θ s r are the Runge-Kutta coefficients in the Shu-Osher formulation as in the previous sections. The matrix M ∈ R N ×N is the usual P 1 mass matrix M ij = ∆x(1 + 3δ ij )/6, while the nonlinear operators N η and N q are given by

N η N q j = ∆xD q q 2 /h j + g 6 0 (2h j + h j+1 )(η j+1 -η j ) -(2h j + h j-1 )(η j -η j-1 )
.

The CIP stabilization J is defined exactly as in [START_REF] Lannes | The Water Waves Problem: Mathematical Analysis and Asymptotics[END_REF], with λ j = u j + gh j . The operators T t and T x are the finite element approximation of the corresponding T operators appearing in [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF], while S is the linear diffusion operator obtained from the discretization of the sponge layer terms. These are reported in Appendix A for completeness. Also for this discrete model we can see the matching with (6), with few differences in the discretization.

The mass matrices can be rewritten as M u = M Φ := M 0 0 M , the dispersion terms and the fluxes can be defined as

X t := 0 0 0 T t , X x u := 0, -µ 2 T x η -T t M -1 N q (η, q) T , F(u) := N η + J (η, λ) N q + J (q, λ) (47)
and the source terms can be written as

S im = S 0 0 S , S ex (u) = ∂ t h iwg 0 . (48)
The main difference with the KdV-BBM model is that here the explicit dispersive terms are nonlinear as they contain a flux term. For more details on the discretization of the spatial operators we refer to [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF].

Benchmarks for weakly dispersive free surface waves

Solitary waves interacting with a submerged bar

We consider the simulation of solitary waves propagating over a trapezoidal bar, defined as in Section 3.4.3, but with maximal height equal to b 0 = 0.2 and points of change of slope equal to {11, 17, 19, 22}. We choose periodic boundary conditions on the domain Ω = [-20, 30] and we center the initial solitary wave at the initial condition in x = 5. Even though this is a traveling solution, the shape of the soliton is smooth and large enough not to encounter problems in the reduction due to the advection character of the solution. We use the following parameters a 0 = 0.2, h 0 = 1, g = 9.81, CFL = 0.5, T = 18. Simulations at various times are displayed in Figures 10(a

Periodic waves over a submerged bar

We consider monochromatic waves propagating in the domain Ω = [0, 35] on a trapezoidal submerged bar of maximum height 0.3. An internal wave generator is positioned at x iwg = 10. Sponge layers are set on both ends of the domain. The wave parameters a 0 and h 0 will be changed along the simulations. For the first test we use a 0 = 0.027, as in [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF] and h 0 = 0.5 (original was h 0 = 0.4 that we will leave for the final simulations). As before, g = 9.81 and CFL = 0.5, while we set T = 40. FOM simulations for different times are depicted in Figures 10(b) and 10(d). This test is particularly complicated to be reduced, as one can heuristically see from the simulations. Indeed, the gradients get steep after the trapezoid and travels all along the right domain.

Here, we compare as for the BBM-KdV benchmarks the computational cost of the system solution with respect to the whole simulation. As before, we observe that the cost of the system is a large part of the simulation. In Table 1 we list the time percentage devoted to the solution of the linear system, while in Table 2 we show the ratio between the solver time and the flux evaluation. For the very optimized Thomas algorithm in numba we observe that the linear solver is not so expensive as in the previous cases, in particular for the non periodic boundary condition case. On the other side, for a generic sparse solver the computational costs of the linear solver are huge with respect to the other operations and it is noticeable comparing it with the cost of the explicit flux evaluation.

Reduction of the enhanced Boussinesq model

For the reduction of the FOM we proceed analogously to the scalar case. The main difference is that, in absence of a clear energy statement allowing a linear reduction, each variable has its own reduced basis associated to the minimization of the L 2 norm as in classical Galerkin projection methods. We discuss the reduction main steps hereafter.

Projection on reduced spaces and linear reduction

We consider two reduced spaces for the two variables: V η ∈ R N h ×N RB η and V q ∈ R N h ×N RB q . The spaces are obtained again with the POD algorithm, using the same tolerance on both variables. Then, we project the equations ( 46) onto these spaces, obtaining

W η,T (M + ∆tS)V η η(s) =W η,T M s-1 r=0 ρ s r (V η η(r) + h (r) iwg ) -h (s) iwg -∆t s-1 r=0 θ s r W η,T N η V η η(r) , V q q(r) + J (V η η(r) , λ (r) ) W q,T (M -T t )V q ψ = - s-1 r=0
θ s r W q,T T t M -1 N q (V η η(r) , V q q(r) ) + T x V η η(r)

W q,T (M + ∆tS)V q q(s) =W q,T MV q s-1 r=0 ρ s r q(r) + ∆tW q,T M q V q ψ -∆t s-1 r=0 θ s r W q,T N q (V η η(r) , V q q(r) ) + J (V q q(r) , λ (r) ) [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF] for every stage s of SSPRK2. The reduced variables are denoted with the • symbol. As already remarked several times, we have no energy allowing to define a minimization problem leading to an efficient linear reduction strategy. Hence we will make use here of classical L 2 projection arguments, and consider Θ to be the identity, and W = V .

Similarly to the scalar case, we can precompute the projected operators and use them in the reduced simulation to save computational time. In this first stage, we obtain the pdROM algorithm. Let us define Mη := W η,T MV η , Mq := W q,T MV q , Tt := W q,T T t V q , Tx := W q,T T x V η , Ŝη := W η,T SV η , Ŝq := W q,T SV q , recalling that h iwg = µ iwg (t, a 0 )f iwg , we also define fiwg := W η,T f iwg . We obtain substituting in [START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF] ( Mη + ∆t Ŝη )η (s) = Mη

( Mq + ∆t Ŝq )q (s) = Mq

s-1 r=0 ρ s r q(r) + ∆t Mq ψ -∆t s-1 r=0 θ s r N q,(r) , (50c)
where the reduced nonlinear fluxes are defined as N η,(r) :=W η,T N η V η η(r) , V q q(r) + J (V η η(r) , λ (r) ) , N q,(r) :=W q,T N q (V η η(r) , V q q(r) ) + J (V q q(r) , λ (r) ) .

Notice that there has been a further projection inside the matrix multiplication W q,T T t M -1 N q ≈ W q,T T t V q V q,T M -1 W q W q,T N q = Tt Mq,-1 N q , (52) in order to have an efficient implementation of the operations without the need of too many operators. Hence, all the matricial operations are reduced and only the nonlinear fluxes have computational costs that scale as an O(N h ). The solution of the reduced system can be done after the simple assembly of the reduced mass and sponge matrices, with dimension N RB . As before, this method will be denoted with pdROM in the benchmarks section and all the terms can be matched with [START_REF] Burman | The edge stabilization method for finite elements in CFD[END_REF] following the same definitions of the FOM discretization. Remark 6 (A more efficient implementation). All of the above manipulations retain the structure of a perturbation of the nonlinear shallow water solver via the use of the auxiliary variable ψ. This is very interesting from the practical point of view for several reasons as in perspective it could be used to enhance an existing shallow water solver. It also allows to more easily embed certain wave breaking closures (see e.g. [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models[END_REF][START_REF] Karczewska | Can simple KdV-type equations be derived for shallow water problem with bottom bathymetry?[END_REF]). However, this reduced model can be more efficiently implemented using a global form that does not isolate the shallow water equations. In particular, instead of (50c) and (50b) we could use the following

( Mq -Tt + ∆t Ŝq )q (s) -( Mq -Tt ) s-1 r=0 ρ s r q(r) + s-1 r=0 θ s r ∆t N q,(r) + Tx η(r) = 0. ( 53
)
Avoiding the explicit evaluation of ψ allows to avoid several projections and reduce the online cost.

Hyper-reduction

The final reduction that can be applied is the EIM on the nonlinear fluxes N η , N q . The same procedures presented in the scalar case can be applied in this case, obtaining an extra level of approximation which allows to get rid of all the computations that scale as an O(N h ). Same caveats hold for this problem: the tuning of the tolerance of the EIM algorithm must be done accordingly to the POD one otherwise instabilities and Runge phenomena might appear, see [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF][START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF][START_REF]On the Formulation of Mass, Momentum and Energy Conservation in the KdV Equation[END_REF][START_REF] Cauquis | Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study[END_REF].

Simulations for enhanced Boussinesq

In this section we study the simulations for the EB model. We follow the modus operandi of Section 4.1 We will show results mainly for the variable η, omitting the ones for q for brevity. In the second stage we will let not only h 0 vary to form the training set, but also a 0 , we will specify their range in the different tests. To make a comparison with a reasonable solver, being the problem more complicated, we will use a sparse scipy solver for the FOM, where the tridiagonal structure of the matrices is not exploited a priori.

Solitary waves over a submerged bar

For this test we set a 0 = 0.2 and h0 = 1. For the first phase we just compress the solutions for a simulation until T = 25 with the POD. The problem is relatively simple and the wave does not show much dispersion. We can see with different choices of N RB η = N RB q the behavior of the pdROM approximated solutions in Figure 11(a). Adding also the EIM to the algorithm we obtain expected results as we can see both in the error behavior in Figure 11(b), and in computational times of Figure 11(c), where we can achieve around 4% to 10% of FOM computational time for pdROM and between 1% and 8% for EIMROM. In this case, the system of the FOM requires more computational time than the scalar case one, as it is composed of more terms. Now, we consider 10 snapshots with randomly chosen parameters h 0 ∈ [0.8, 1.2] and a 0 ∈ [0.16, 0.24]. This training set is used to compute the POD and EIM basis functions. Simulating pdROM and EIMROM for a 0 = 0.2 and h 0 = 1, with relatively few basis functions the error is already low, see Figure 12(a) and the computational costs stay in the same range between 1% and 10% of FOM computational times, see Figure 12(b). In Figure 12(c) we see that with only N RB η = 20, N RB q = 20 for 0.03 POD tolerance the pdROM already approximates very well the solution and adding the EIM with N EIM η = 47 and N EIM q = 60 for 0.0003 EIM tolerance results in an accurate approximation of the FOM in just 7% of the FOM computational time.

Considering h 0 = 0.76 and a 0 = 0.252 slightly outside the training set we can already see larger errors (one order of magnitude) in Figure 13(a) and comparable computational times in Figure 13(b). In Figure 13(c) we see that the FOM is a bit more oscillating than in the previous case and that the pdROM and EIMROM, for the same parameters used above, struggle more with approximating accurately the solution, still being not too far from the FOM solution.

For N RB η = 28, N RB q = 28 relative to a POD tolerance of 0.01 and N EIM η = 44, N EIM q = 50 for EIM tolerance of 0.001, we see in Figure 14 how the parameters a 0 and h 0 influence the relative error of the pdROM and EIMROM approximations. The structure is similar to the scalar case, where increasing the nonlinearity of the problem the oscillations rise and the error follows. Again, we see that the EIMROM do not obtain such an accurate solution as pdROM for problems close to the training set, though using much less computational time (around 3 times faster than pdROM).

Monochromatic waves on a submerged bar

The last test is very challenging as its solutions are very close to show steep gradients and transport phenomena. Moreover, the water level is often close to the bathymetry and, if not well represented, the approximation can result into invalid values. Fortunately, the steep gradients are all close to the end of the trapezoid and with not so many basis functions it is possible to approximate them. We take a 0 = 0.027, h0 = 0.5 and final time T = 40. Taken the POD over the time evolution of this FOM, we approximate in Figure 15(a) with pdROM the same evolution for different N RB η = N RB q of the same problem. We see that with 40 basis functions the pdROM approximation is close to the FOM one, but only for 60 basis functions we do not observe any large deformation. In Figure 15(b) we see more precisely the error decay for different N RB and different EIM tolerances. The error is quite large, and the decay with the EIM approximation is very slow. It must be noticed that in many of these EIMROM simulations, the water height reaches negative values and stop at an earlier time, in those cases the error is not reported or it is larger than the scale. We see for an EIM tolerances of 10 -4, 3 • 10 -5 and 10 -6 in Figure 15(d) how the error for the first one is very large, but also for the other two cases the error does not decrease a lot in particular in the area after the trapezoid, where the steepest part of the solutions are located. The computational costs for this problem range between 4% and 30% for pdROM and between 1% and 30% for EIMROM as shown in Figure 15(c).

Introducing more parameters in the training set, with h 0 ∈ [0.4, 0.6] and a 0 ∈ [0.0216, 0.0324], leads to a much more rich manifold of solutions. This implies that the eigenvalue decay of the POD gets slower. For the test parameter h0 and a 0 = 0.027, we observe in Figure 16(b) that the computational costs for pdROM and EIMROM are very similar and they are around the 5% of the FOM computational time to obtain an error of 10 -1 , around 15% of the time for an error of 10 -2 and around 30% for an error of 10 -3 . For very accurate results with error lower than 10 -4 a computational time of more of the 30% of the FOM computational time may be required. The pdROM simulation in Figure 16(a) shows the accurate result for tol P OD = 0.001, N RB η = 99, N RB η = 92 and tol EIM = 0.00003, N EIM η = 213 N EIM q = 291. The number of EIM basis functions is large in order to obtain a good approximation of the solution and almost not gaining anything in terms of computational time with respect to the pdROM. This number must be large in order to avoid the simulation to oscillate so much that it hits negative values and it crashes the simulation. Again, this behavior has been observed in many other works [START_REF] Peherstorfer | Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling[END_REF][START_REF] Yano | Discontinuous Galerkin reduced basis empirical quadrature procedure for model reduction of parametrized nonlinear conservation laws[END_REF][START_REF] Filippini | A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up[END_REF][START_REF]On the Formulation of Mass, Momentum and Energy Conservation in the KdV Equation[END_REF][START_REF] Cauquis | Lax-Wendroff Schemes with Polynomial Extrapolation and Simplified Lax-Wendroff Schemes for Dispersive Waves: A Comparative Study[END_REF] and we do not aim at finding better hyper reduction algorithm in this work.

For the parameter outside the training set, we chose a 0 = 0.027 and h 0 = 0.4, which correspond to the test case in [START_REF] Reiss | The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena[END_REF] which was validated with experimental data. This simulation is quite challenging also for pdROM and when it is not enough resolved it can have oscillations in h going below 0. This happens for example for tol P OD = 5 • 10 -2 , as shown from the incomplete curves in Figure 17(a). Increasing the dimension of the reduced space we overcome this issue. Computational costs are similar to the ones discussed above. In Figure 17(c) we plot the simulation for tol P OD = 0.005, N RB η = 99, N RB η = 92 and tol EIM = 0.00003, N EIM η = 213 N EIM q = 291. Here it is the EIMROM simulation with the largest tolerance that does not crashes along the simulation. The computational time reduction due to the EIM at this stage is almost negligible with respect to the pdROM and the approximation is very similar. Both simulations oscillate more than the original FOM after the trapezoid at this resolution, though using only 15% of the computational time of the FOM. To reach errors of the order of 10 -1 we need computational costs of around the 25% of the FOM ones as one can observe comparing Figure 17 In Figure 18 we see the strong influence that h 0 and a 0 have on these simulations. In this case, it is crucial the level of h 0 . When this level is too low, we fall in another regime and all the hypotheses made on the dispersive equations are not valid anymore. Moreover, the simulations risk to hit negative values and break in this regime. Above this regime, for pdROM we can simply say that as the oscillations and the nonlinearity increase, the error increases. On the other side, the dispersion term is actually helping in reducing the error in the reduced space. For EIMROM we observe more areas where the algorithm struggles with obtaining good results, also for not too small h 0 values.

Conclusion

In this work we have proposed some strategies for model order reduction for dispersive waves equations. Applications have been shown to a BBM-KdV type model, and to the enhanced Boussinesq equations of Madsen and Sørensen. Both models contain hyperbolic and dispersive terms that can be decoupled in the numerical solver. The dispersive terms can be obtained solving an elliptic and linear problem, which is well suited to be reduced with standard projection based reduction techniques after the choice of a reduced basis space. We apply the proper orthogonal decomposition on some snapshots obtained for different times and parameters to obtain a reduced space. When the gradients of the solutions are not so steep and the nonlinear character is not the dominant one, we can obtain large reduction with the pdROM in computational times, up to 20 times less, and still obtain a reliable solution. Introducing a second level of approximation with the empirical interpolation method, this factor can increase up to 100, but stability issues may arise. Hence, it is not always obvious how to choose how many interpolation functions are needed. There are some clear bounds that cannot be overtaken when using these simple reduction techniques: the nonlinearity cannot be too pronounced and the water level must be far enough from zero. Nevertheless, the results are encouraging in particular knowing that already for one dimensional problems we can can obtain good reduction.

This idea can be easily applied to a pre-existent hyperbolic algorithm, as a shallow water code, with a cheap additional term which takes into account of the dispersive effect. This can, in fact, be used without the burden of having to solve a large linear system. In the future, we aim to apply this algorithm to more complicated two dimensional problems, where we expect larger reduction in computational times (as the FOM will become more expensive), and to selectively switch on and off the dispersive term into a shallow water code, in order to allow to pass between different regimes, according to the solution shape.

To our knowledge this is the first work in this direction, and the initial results are quite promising. Several future improvements and developments are foreseen both involving improved reduction strategies, and more complex wave dynamics and models. Notable the extension to multiple space dimensions, as well as breaking waves are under investigation.

Figure 1 :

 1 Figure 1: Water waves: notation

  are h 0 = 1, g = 9.81, a 0 = 0.04, domain [0, 20π], CFL = 0.1 and final time T = 20. Initial and final simulations are shown in Figure 2(b).

Figure 2 :

 2 Figure 2: Initial and final solution for tests.

Figure 3 :

 3 Figure 3: Singular values decay

Figure 4 :Figure 5 :

 45 Figure 4: Propagation of a periodic monochromatic wave. pdROM and EIM reduction: simulation, errors and computational time, varying POD and EIM dimensions

  ROM errors time-parameter reduction on parameter and η 0 in the training set ROM computational time timeparameter reduction on parameter and η 0 in the training set (legend above) ROM computational time timeparameter reduction on parameter and η 0 outside the training set (legend above)

Figure 7 :

 7 Figure 7: Undular bore propagation. Error plot varying the parameters a 0 and h 0 for pdROM N RB = 46

  ROM errors time-parameter reduction on parameter and η 0 in the training set ROM errors time-parameter reduction on parameter and η 0 outside the training set EIM ROM computational times time-only reduction (legend above)

  ROM computational time timeparameter reduction on parameter and η 0 in the training set (legend above)

  ROM computational time timeparameter reduction on parameter and η 0 outside the training set (legend above) ROM solutions time-parameter reduction with tol P OD = 0.01, N RB = 72 and tol EIM = 0.03, N EIM = 152 on parameter and η 0 in the training set

Figure 9 :

 9 Figure 9: Solitary waves interacting with a submerged bar. Error plot varying the parameters a 0 and h 0 for pdROM N RB = 72 and EIMROM with N EIM = 227

  ) and 10(c).

22

 22

  Periodic waves over a submerged bar

( d )Figure 10 :

 d10 Figure 10: Full order model solutions for EB system tests.

r

  Tt Mq,-1 N q,(r) + Tx η(r) ,

Figure 11 :

 11 Figure 11: Solitary waves over a submerged bar. pdROM and EIM (time-only) reduction: simulation, errors and computational time, varying POD and EIM dimensions

Figure 12 :

 12 Figure 12: Solitary waves over a submerged bar. pdROM and EIM (time-parameter) reduction: simulation, errors and computational time, varying POD and EIM dimensions. Parameter inside the training domain.

Figure 13 :

 13 Figure 13: Solitary waves over a submerged bar. pdROM and EIM (time-parameter) reduction: simulation, errors and computational time, varying POD and EIM dimensions. Parameter outside the training domain.

Figure 15 :

 15 Figure 15: Monochromatic waves on a submerged bar. pdROM and EIM (time-only) reduction: simulation, errors and computational time, varying POD and EIM dimensions

Figure 16 :

 16 Figure 16: Monochromatic waves on a submerged bar. pdROM and EIM (time-parameter) reduction: simulation, errors and computational time, varying POD and EIM dimensions. Parameter inside the training domain.

  0.01 EIM tol = 0.005 EIM tol = 0.001 EIM tol = 0.0001 EIM tol = 3e-05 EIM tol = 1e-05 Errors (a) ROM errors time-parameter reduction on parameter outside the training set

Figure 17 :

 17 Figure 17: Monochromatic waves on a submerged bar. pdROM and EIM (time-parameter) reduction: simulation, errors and computational time, varying POD and EIM dimensions. Parameter outside the training domain.

  (a) and Figure 17(b) for tol P OD = 2 • 10 -3 .

Table 1 :

 1 Estimate percentage of computational time for solving the linear systems for the three benchmark problems: comparison of the implementation of Thomas algorithm all in numba framework and sparse solver by scipy with fluxes computed in numba

		numba and Thomas	scipy and numba
	N h	3.4.1	3.4.2	3.4.3	3.4.1	3.4.2	3.4.3
	2000 86.4% 63.4% 85.2% 86.6% 84.4% 86.6%
	4000 84.1% 62.4% 83.9% 86.2% 83.5% 88.8%
	8000 86.1% 61.5% 85.0% 94.1% 85.3% 91.0%

Table 2 :

 2 Estimate of ratio of computational time of the linear system over computational time of the nonlinear fluxes for the three benchmark problems: comparison of the implementation of Thomas algorithm all in numba framework and sparse solver by scipy with fluxes computed in numba

		numba and Thomas	scipy and numba
	N h	3.4.1 3.4.2 3.4.3 3.4.1 3.4.2 3.4.3
	2000	11.2	3.5	11.3	72.4	52.1	71.8
	4000	11.8	3.7	11.6	73.4	54.9	74.3
	8000	11.2	3.7	12.8	78.4	56.1	75.8

3.4.3. Solitary waves interacting with a submerged bar

This test consists of a soliton evolving on a ramp bathymetry defined by b

Table 3 :

 3 Estimate percentage of computational time for solving the linear systems for the three benchmark problems: comparison of the implementation of Thomas algorithm all in numba framework and sparse solver by scipy with fluxes computed in numba

		numba and Thomas	numba and scipy
	N h	5.3.1	5.3.2	5.3.1	5.3.2
	2000 75.8%	39.1%	74.3%	77.7%
	4000 72.1%	32.8%	87.4%	80.8%
	8000 69.1%	27.7%	94.5%	81.1%

Table 4 :

 4 Estimate of ratio of computational time of the linear system over computational time of the nonlinear fluxes for the three benchmark problems: comparison of the implementation of Thomas algorithm all in numba framework and sparse solver by scipy with fluxes computed in numba

		numba and Thomas	numba and scipy
	N h	5.3.1	5.3.2	5.3.1	5.3.2
	2000	6.8	2.25	39.8	31.5
	4000	6.7	2.18	41.1	31.5
	8000	6.9	2.15	39.0	32.1

  Solitary waves over a submerged bar. Error plot varying the parameters a 0 and h 0 for pdROM N RB

								Reduced solutions at final time			q
	0.0											0.06
	0.1											0.04		10 1
	0.5 0.4 0.3 0.2	16	18	20	22	24	26	28	30 ROM N=10 FOM ROM N=20 ROM N=40 ROM N=60 ROM N=100	32	0.06 0.04 0.02 0.00 0.02 Error	16 10 4 10 3 10 2 10 1 10 0	18	20 ROM EIM tol 0.0003 22 24 20 40 60 Dimension reduced basis 26 28 80 100 ROM EIM tol 0.0001 ROM EIM tol 3e-05 ROM EIM tol 1e-05 ROM EIM tol 3e-06 ROM EIM tol 1e-06 pdROM	30	120	32	140		ROM Time/FOM Time	10 1	20	40	60 Dimension reduced basis 80 100 ROM EIM tol 0.0003 120 140 ROM EIM tol 0.0001 ROM EIM tol 3e-05 ROM EIM tol 1e-05 ROM EIM tol 3e-06 ROM EIM tol 1e-06 pdROM
	(a) ROM solutions time-only reduction	(b) EIM ROM errors time-only reduc-	(c) EIM ROM computational times
												tion												time-only reduction
	0.06 0.08						FOM pdROM 50 50 EIM 71 125 ROM 50 50		0.04 0.05						FOM pdROM 50 50 EIM 100 161 ROM 50 50	0.04 0.05	FOM pdROM 50 50 EIM 170 231 ROM 50 50
	0.04												0.03												0.03
	0.06 0.04 0.02 0.00 0.02												0.01 0.00 0.01 0.02												0.01 0.00 0.01 0.02
	0.08												0.02												0.02
	16	18	20	22		24	26	28	30	32			16	18	20	22	24	26		28	30	32		16	18	20	22	24	26	28	30	32
						Log10 of relative error PDROM 28 28, POD tol 0.01									Log10 of relative error EIM 44 50 ROM 28 28, EIM tol 0.001, POD tol 0.01
		1.4												Training		0.055		1.4									Training	0.055
																	0.435										0.435
		1.2															0.926	1.2									0.926
		1.0															1.416	1.0									1.416
	h0																		h0
		0.8															1.906	0.8									1.906
																	2.397										2.397
		0.6															2.887	0.6									2.887
		0.10 0.4		0.15		0.20		0.25 asurf	0.30	0.35	0.40	3.377	0.10 0.4		0.15	0.20	0.25 asurf	0.30	0.35	0.40	3.377
	Figure 14: η = 31
	and N RB	q = 30 and EIMROM with N EIM	η = 40 and N EIM	q = 51

Acknowledgments D. T. has been funded by an Inria Postdoc in Team Cardamom and by a SISSA Mathematical Fellowship. M. R. is a permanent member of Inria Team Cardamom.

Declarations of interest: none.

We report here the expressions of the finite element approximations of some of the operators arising in the Madsen and Sørensen model. The first is T t ∈ R N ×N , representing the matrix discretization of the terms T t . It can be split into two matrices. The first one is simply defined as a tridiagonal matrix with entries for the jth row B h2 j [1 -2 1], while the second one has the following entries for the jth row

The sum of the two defines the dispersion matrix

. The operator can be approximated by a matrix T x with the multiplication T x η. It can be conveniently written with two matrices T x,1 and T x,2 such that T x η = T x,1 w + T x,2 w, where T x,1 has the following 3 entries on the j-th row in the 3 main diagonals

and where T x,2 has the following 3 entries on the j-th row in the 3 main diagonals Finally, the discretization of this sponge layer term leads to an operator S which is a tridiagonal matrix with entries for the jth row

Appendix B. Computational costs of different solvers

We compare here the computational costs of different sparse linear solvers and different strategies for the reduction. We will refer to the KdV-BBM problem [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF], similar results can be drawn from the enhanced Boussinesq system [START_REF] Torlo | Model reduction for advection dominated hyperbolic problems in an ALE framework: Offline and online phases[END_REF].

First of all, let us compare different sparse linear solvers for the FOM. We have used Thomas algorithm in the previously proposed examples. This linear solver very optimized and tailored towards the problem we are solving. As soon as the discretization gets more involved or as soon as we move to more dimensions, we would lose the tridiagonal structure of the system matrices and it would become impossible to apply Thomas algorithm. In order to have a fairer comparison, in Figure 19(a) we show the ROM computational costs over the FOM computational costs when the FOM uses Thomas algorithm implemented in numba and when the FOM uses the sparse solver of scipy. We observe a huge difference between the two solvers of a factor between 4 and 6. This would make the pdROM computational costs even more appealing for a general sparse solver implementation in the FOM, passing from 20% of the computational time to only 4%! Secondly, we want to prove again that the main computational cost of the FOM is given by the linear system and that the advantage obtained with the pdROM is focused on this step of the algorithm. Indeed, we can compare the performance of the pdROM with the scheme given by the reduction only of the Φ equation, as suggested in Remark 3. In that algorithm, Φ is computed by projecting the RHS of its equation onto the reduced basis space, then solving the reduced system and, finally, reconstructing the full Φ. Its computational costs are comparable to the one of the pdROM. In Figure 19(b) we plot the ratio of these costs for some N RB for all the previous tests. We can see that they are very close to one and the fluctuations around that are probably due to the low precision in measuring the times. We remind that the ROM where only Φ is reduced has larger errors than the pdROM. This is maybe due to a mismatch of the spaces along the computations that allows the propagation of spurious modes. We are still investigating the phenomenon. Hence, in general, we would recommend a fully projected approach.