# The Simple Finite Volume Lax-Wendroff Weighted Essentially Nonoscillatory Schemes for Shallow Water Equations with Bottom Topography

Authors: Changna Lu; Luoyan Xie; Hongwei Yang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

A Lax-Wendroff-type procedure with the high order finite volume simple weighted essentially nonoscillatory (SWENO) scheme is proposed to simulate the one-dimensional (1D) and two-dimensional (2D) shallow water equations with topography influence in source terms. The system of shallow water equations is discretized using the simple WENO scheme in space and Lax-Wendroff scheme in time. The idea of Lax-Wendroff time discretization can avoid part of characteristic decomposition and calculation of nonlinear weights. The type of simple WENO was first developed by Zhu and Qiu in 2016, which is more simple than classical WENO fashion. In order to maintain good, high resolution and nonoscillation for both continuous and discontinuous flow and suit problems with discontinuous bottom topography, we use the same idea of SWENO reconstruction for flux to treat the source term in prebalanced shallow water equations. A range of numerical examples are performed; as a result, comparing with classical WENO reconstruction and Runge-Kutta time discretization, the simple Lax-Wendroff WENO schemes can obtain the same accuracy order and escape nonphysical oscillation adjacent strong shock, while bringing less absolute truncation error and costing less CPU time for most problems. These conclusions agree with that of finite difference Lax-Wendroff WENO scheme for shallow water equations, while finite volume method has more flexible mesh structure compared to finite difference method.


---

## Body

## 1. Introduction

In this paper, the simple finite volume WENO (weighted essentially nonoscillatory) scheme with Lax-Wendroff-type time discretization is proposed to numerically solve the system of shallow water equations with bottom topography influence in source terms.It has been an important work to search for the solutions of nonlinear differential equations due to their rich mathematical structures and features [1–4] as well as important applications in fluid dynamics and plasma physics [5–9]. The solutions of differential equations can be simply divided into two kinds, exact solutions and numerical solutions. For the exact solutions, a great number of methods have been introduced, like Painlevé test, Darboux Transformation, bilinear method, symmetry method [10], and so on [11–19], while for the forced nonlinear differential equations, such as shallow water equations with topography influence in source terms, it is very difficult to get the exact solutions and we have to derive the numerical solutions and analyze the topography forced effect.WENO is a procedure of spatial discretization for partial differential equation (PDE); in other words, it is a numerical method to discretize the derivative terms in space. The WENO scheme is an important high order accuracy numerical method, especially to simulate strong shocks and contact discontinuous and complicated smooth solutions, as they can keep high order and nonoscillatory characteristic for both continuous and discontinuous solutions. WENO scheme is popularly used in many communities [20–23] in recent years; it has attracted much attention in CFD (computational fluid dynamics); especially for shallow water flows, finite volume WENO scheme made important contributions as more flexibility of mesh shape. Xing and Shu [24] introduced a finite volume WENO wave propagation scheme on rectangular mesh. Lu et al. [25] investigate the performance of finite volume WENO scheme for the system of shallow water equations on unstructured triangle meshes; simulation of a tidal bore on Qiantang river is performed.In recent years, WENO scheme has been generalized for hyperbolic conservation laws [22, 26, 27] after the first WENO scheme was originally derived in [28] for the third-order finite volume frame based on ENO type schemes [29, 30], such as CWENO and hybrid WENO [31] arising from various applications. In 2016, a successful type of WENO [32] was proposed to approximate hyperbolic conservation laws. The new WENO reconstruction is a convex combination of two linear polynomials with a fourth-degree polynomial using the same five-point big stencil as in the classic WENO fashion; the linear weights are constants, which are positive, and their summation equals one so they are more simple and easily extended to multidimensions in engineering applications. We call the new type of WENO as simple WENO (SWENO).For time-dependent problems, there are mainly two different ways to approximate derivative with respect to time [33–35]. One common approach is well known ODE solver, such as Euler’ method, Runge-Kutta method, Adams’ multilevel method, and Radau scheme. These approaches have been popularly used for the advantages of good stability and simplicity in idea and codes. The disadvantage is that the highest accuracy order for total variation diminishing (TVD) Runge-Kutta method is fourth order. Another way is via the Lax-Wendroff-type procedure; the idea is converting all the time derivatives into spatial derivatives using PDE and Taylor expansion with respect to time, then discretizing the spatial derivatives using numerical methods. The advantages are smaller stencil and frequent use of the original PDE. The disadvantages are complicated in formulation and coding. Both approaches can produce any order accuracy in time theoretically. The finite volume Lax-Wendroff time discretization formula is a little simpler than that of the finite difference frame. Qiu and Shu [34] derived a scheme consisting of Lax-Wendroff-type time discretization and the finite difference WENO scheme in compressible gas dynamics; the scheme has important contributions to reducing CPU cost and maintaining the nonoscillatory characteristic. Lu and Qiu [36] also developed Lax-Wendroff with finite difference WENO scheme for shallow water equations; similar results were obtained.In this paper, we developed the fifth-order SWENO finite volume scheme with the third-order Lax-Wendroff-type time discretization to simulate shallow water flows. In the scheme we follow the ideas in [32, 37] about simple finite volume WENO scheme and the algebraic prebalance method for the flux and the source terms [36, 38] in shallow water equations. An outline of the paper is as follows: in Section 2 we describe details of the discretization with finite volume SWENO scheme and Lax-Wendroff-type time discretization for shallow water equations. In Section 3, a range of 1D and 2D numerical results show that the scheme has accuracy, resolution, efficiency, and robustness. Finally, Section 4 contains conclusions.


## 2. Description of Numerical Model

### 2.1. The Shallow Water Equations

The 2D system of shallow water equations with conservative form is(1)Ut+FUx+GUy=SU,with conservational vector (2)U=D,Du,DvT,fluxes(3)FU=Du,Du2+gD22,DuvT,GU=Dv,Duv,Dv2+gD22T,and source terms(4)SU=0,-gDbx,-gDbyT,where D is the total water depth; u and v are the depth-averaged water velocities in the x- and y-directions, respectively; t is the time; g is the gravity constant; b is the topography function above; bx and by are derivatives with respect to x and y; they are used to describe the bed slopes stresses. Friction on the bottom, wind stress, and Coriolis term could also be considered as parts of source terms.We rewrite (1) in quasilinear form as(5)Ut+AUUx+BUUy=SU,where (6)AU=∂F∂U=010c2-u22u0-uvvu,BU=∂G∂U=001-uvvuc2-v202v.The eigenvalues of Jacobian matrixes A and B are(7)λA=u-c,u,u+c;λB=v-c,v,v+c.The normalized right and left eigenvectors of Jacobian matrixes A and B are as follows:(8)RA=101u-c0u+cv1v,LA=u+c2c-12c0-v01-u-c2c12c0;RB=101u1uv-c0v+c,LB=v+c2c0-12c-u10-v-c2c012c.We can easily check that RALA=I, RBLB=I, which is an important characteristic for hyperbolic conservation law. Neglecting variables in y-direction, we have the 1D conservative unsteady shallow water equations: (9)Dt+Dux=0,Dut+Du2+12gD2x=-gDbx.Similar way is that we can get Jacobian matrix and other information for 1D case.


### 2.2. The Fifth-Order Finite Volume SWENO Schemes

The fifth-order finite volume SWENO scheme will be used to numerically solve the spatial derivatives of shallow water equations [32, 37]. In this section, the basic procedure of SWENO will be described as short overview for the following 1D scalar hyperbolic conservation law:(10)ut+fux=0.DenoteIi=[xi-1/2,xi+1/2] as the ith cell, a control volume for finite volume method, centered on xi; we assume that the grid points are uniform with Δx=Δxi=xi+1/2-xi-1/2. For a finite volume scheme, we evolve the proximate average value u¯i of u¯(xi,t)=1/Δx∫Iiu(x,t)dx at mesh cell Ii in time. Using integral average, the conservation law (10) is rewritten into (11)ddtu¯xi,t=-1Δxfuxi+1/2,t-fuxi-1/2,t.f(u(xi+1/2,t)) is approximated by numerical flux f^i+1/2=f^(ui+1/2L,ui+1/2R). The popular Lax-Friedrich flux is used for its simpleness and monotone, which is given by (12)f^ui+1/2L,ui+1/2R=12fui+1/2L+fui+1/2R-αui+1/2R-ui+1/2L,where α=maxu|f′(u)|; ui+1/2L,R are approximations to u(xi+1/2,t) of left side and right side at xi+1/2 considering the discontinuous solution. We call this process reconstruction; here the reconstructions of interface values depend on the cell averages by the fifth-order SWENO reconstruction procedure.Onlyui+1/2L will be introduced in detail as follows. For the fifth-order SWENO scheme, three approximated point values of ui+1/2L should be constructed. Choose the first big stencil {Ii-2,Ii-1,Ii,Ii+1,Ii+2}; a fourth-degree polynomial p1(x) can be obtained by requiring the same mean value on the five cells. The first approximated point value p1(xi+1/2) for ui+1/2L is(13)ui+1/2L1=-120u¯i-2+920u¯i-1+4760u¯i-1360u¯i+1+130u¯i+2.Choose the second stencil {Ii-1,Ii}; another linear polynomial p2(x) can be obtained with the same idea. The second approximation p2(xi+1/2) is(14)ui+1/2L2=-12u¯i-1+32u¯i.Choose the third stencil {Ii,Ii+1}; a linear polynomial p3(x) can be obtained and the third approximation p3(xi+1/2) is(15)ui+1/2L3=12u¯i+12u¯i+1.Convex combination coefficients are the nonlinear weights. Before that we have to compute the smoothness indicators of the three above polynomials first. For the fifth-order SWENO reconstructions, the convex combination of all these three-point values is(16)ui+1/2L=ω11d1ui+1/2L1-d2d1ui+1/2L2-d3d1ui+1/2L3+ω2ui+1/2L2+ω3ui+1/2L3.The nonlinear weights ωm  (m=1,2,3) are defined by (17)ωm=αm∑s=13αs,αs=ds1+τε+βs,τ=β1-β2+β1-β322,and here ε is meaningless, just to avoid the denominator to be zero, so ε is a positive constant and typically ε=10-6, ds are the linear weights, and “smoothness indicators” βs are used to measure the smoothness of the three polynomials in cell Ii. The smaller the value of smooth indicator is, the smoother the polynomial is. We use the same smooth indicator formula for k-degree polynomial p(x) in [39](18)β=∑l=1k∫xi-1/2xi+1/2Δxi2l-1∂lpx∂xl2dx.The explicit expression of the three smoothness indicators βs is given by(19)β1=1144u¯i-2-8u¯i-1+8u¯i+1-u¯i+22+115600-11u¯i-2+174u¯i-1-326u¯i+174u¯i+1-11u¯i+22+7812880-u¯i-2+2u¯i-1-2u¯i+1+u¯i+22+14214611310400u¯i-2-4u¯i-1+6u¯i-4u¯i+1+u¯i+22;β2=u¯i-1-u¯i2;β3=u¯i-u¯i+12.For the fifth-order SWENO scheme on uniform mesh, the important idea is that the linear weights are your chosen positive constants; the only condition is that the summation is one. So the procedure of linear weights in classical WENO is avoided, and the negative weights are prevented. In this paper, we use the positive linear weights ds as (20)d1=0.98,d2=0.01,d3=0.01.On uniform mesh, the reconstruction ofui+1/2R is mirror symmetry on big stencil {Ii-1,Ii,Ii+1,Ii+2,Ii+3} to that of ui+1/2L. For shallow water equations, the reconstruction is operated in local eigenspace for more robust, so the local characteristic decomposition should be used, in which properties of eigenspace and Roe’s averages are needed.For 2D finite volume case, we cannot simply proceed in a dimension-by-dimension fashion as in finite difference method. The spatial surface integration on rectangular cellIi,j=[xi-1/2,xi+1/2]×[yj-1/2,yj+1/2] can be converted to line integration by Gauss theorem. In order to get the fifth order, we can use three-point Gaussian quadrature for line integration of every rectangular cell like in Figure 1. For example, to calculate numerical flux in x-direction, three Gaussian point values on every vertical line should be reconstructed from cell average values in y-direction first, then left and right Gaussian point values should be reconstructed by the SWENO scheme, then the numerical fluxes are calculated, and last the cell average value can be approximated on next time level.Figure 1
A sketch of SWENO reconstruction on a rectangular cell.


### 2.3. Prebalance of the Flux and Bottom Topography

We adopt the prebalance of the flux gradients and bottom topography in source term by Rogers et al. [38] using algebraic technique. Taking 1D shallow water equations (9), for example, the prebalanced formulation is(21)U~t+F~Ux=S~U,with (22)U~=ζ,DuT,F~=Du,Du2+gζ2+2hζ2T,S~U=0,-gζbxT.The main algebraic technique is to separate U into(23)U=Ueq+U~,with (24)Ueq=h,0T,U~=U-Ueq=ζ,huT,in which ζ is free surface elevation, the height of water from the still level h+b to free surface elevation; h is still water height, the height of water from bottom to still water level; D=ζ+h is the total water depth; see Figure 2 [36].Figure 2
A sketch of variables in shallow water equations.It is easy to verify that modified system (21) and system (9) have the same Jacobian matrix. Then use the finite volume SWENO discretization for (21).The specific value of the still water level depends on examples; we will describe the chosen datum below. It is perfectly reasonable to choose a fixed horizontal datum elsewhere and derive the prebalanced equations using a stage-discharge approach. Mostly, we choose the still water level between the maximum of bottom function and the minimum of free surface.


### 2.4. The Lax-Wendroff-Type Discretization for 1D Shallow Water Equations

We focus on the procedure of Lax-Wendroff-type time discretization scheme for 1D prebalanced shallow water equation (21) [33, 34]; for simplicity we still describe (21) as(25)Ut+FUx=SU.In this process, SWENO finite volume method will be used. Denote time level as tn+1=tn+Δt, Δt is time step, and U¯in is the approximate average value of U on Ii at the nth time level. The important idea in Lax-Wendroff-type time discretization procedure is to replace time derivatives by the spatial derivatives using Taylor expansion and PDE (25) [34]. By a temporal Taylor expansion with respect to time t for conservative vector U we obtain(26)Ux,t+Δt=Ux,t+ΔtUt+Δt22Utt+Δt36Uttt+Δt424Utttt+⋯,and the sliding average is(27)U¯x,t+Δt=U¯x,t+ΔtUt¯+Δt22Utt¯+Δt36Uttt¯+Δt424Utttt¯+⋯,where Ut¯ means mean value of time derivative Ut in x-direction, which will be given in detail as follows. In this paper we would like to obtain the third-order accuracy in time, so we need to approximate those terms on the right hand side in (27) until the third time derivatives and neglect other terms behind the third time derivatives. Theoretically, we can obtain any accuracy order in time and just need to approximate those terms until the according time derivatives.Subject to shallow water equation (25), by Lax-Wendroff type dime discretization, the time derivatives in (26) are turned to the spatial derivatives as follows: (28)Ut=-FUx+SU;Utt=-F′UUtx+SUt=F′U2Ux-F′USUx+S′U-F′UUx+SU;Utttt=-F′′UUt2+F′UUttx+SUtt=-3F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2SUUx+F′USUtx+-F′UUx+SU-S′′UF′UUx-S′UF′′UUx+S′′USU+S′U2-S′UF′U-F′′UUx2-F′UUxx+S′UUx.The abbreviated form is used like F′′(U)(Ut)2, which means F′′(U)UtUt; the details will be described later. By the way, the above vector source term is easy to be generalized to any dimensions; indeed, as the source term S(U)=[0,-gζbx]T, we only need to handle time derivative as scalar case in practice. For simplicity, we just keep time derivative of source term temporally as follows.So (26) can be written as(29)Ux,t+Δt=Ux,t-ΔtFU-Δt2F′U2Ux-F′USU+Δt263F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2SUUx+F′USUtx+ΔtSU+Δt2SUt+Δt26SUtt+⋯.Then the sliding average (27) with third order in time will be approximated by(30)U¯in+1=U¯in-ΔtΔxFUi+1/2-FUi-1/2+ΔtSU¯+Δt2SUt¯+Δt26SUtt¯,with(31)FU=FU-Δt2F′U2Ux-F′USU+Δt263F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2SUUx+F′USUt.F(U) is approximated by numerical flux F^, such as F(Ui+1/2)≈F^i+1/2=F^(Ui+1/2L,Ui+1/2R). With this form, the finite volume scheme has more flexibility in terms of mesh structure than that of finite difference method [36].Using the fifth-order SWENO described in Section2.2, we reconstruct point values at interface Ui+1/2L,R; meantime, we calculate point values of (Ux)i+1/2L,R and (Uxx)i+1/2L,R. We have to point out that only one order lower accuracy approximation for Ux is needed compared with that of U, as the extra Δt factor multiplied by the second term in flux F in (31). In the same way, we need the third accuracy order of approximation for Uxx, which is in the third term in (31). Indeed, for linear polynomial p2(x),p3(x), the first derivative is constant, the second derivative is zero, the nonlinear weights do not need to calculate for reconstruction of  Ux,Uxx, and only the same linear weight is used, which makes the process easier. For example,(32)Uxi+1/2L=p1′xi+1/2.For the fourth-degree polynomial p1(x), the computation is a little more complicated.In reconstruction ofF(U), we also need Jacobian matrix(33)F′U=AU=01c2-u22u,F′′U=A′U=00g+2u2h-2uh00-2uh2h,and, for simplicity, let F′′(U)=(F1,F2), where F1, F2 are the two 2×2 matrixes; F′′(U)(U′)2=F1U′,F2U′U′ is a 2×1 vector.In order to fit the scheme to problems with discontinuous or big slop in bottom topography, the bottom topography effect is included in numerical flux. In other words, we use SWENO reconstruction for bottom function and then use Gaussian quadrature for the source termsS¯i=[0,-gζ¯bx¯]iT, in which ζ¯ is known. This method can be generalized to varying bottom topography with time. In St¯i=[0,-gζt¯bx¯]iT and Stt¯i=[0,-gζtt¯bx¯]iT, ζt and ζtt have been calculated as part in flux F.For simplicity, we use SWENO5-LW3 to denote the scheme with the fifth-order finite volume SWENO scheme in space and the third-order Lax-Wendroff-type method in time and WENO5-RK3 to denote classical WENO scheme and the third-order Runge-Kutta methods for comparison.


### 2.5. The Lax-Wendroff-Type Discretization for 2D Shallow Water Equations

In 2D case, we need more information; this section aims to display the frame and idea. For 2D shallow water equation(34)Ut+FUx+GUy=SU,the sliding average of Taylor expansion is(35)U¯x,y,t+Δt=U¯x,y,t+ΔtUt¯x,y,t+Δt22Utt¯x,y,t+Δt36Uttt¯x,y,t+Δt424Utttt¯x,y,t+⋯.Using shallow water system (34), by Lax-Wendroff type dime discretization, the time derivatives in (35) are turned to the spatial derivatives as follows. For simplicity we keep time derivative of source term temporally.(36)Ut=-FUx-GUy+SU;Utt=-F′UUtx-G′UUty+SUt=F′U2Ux-F′USUx+G′U2Uy-G′USUy+SUt;Uttt=-F′′UUt2+F′UUttx-G′′UUt2+G′UUtty+SUtt=-3F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2S′UUx+F′USUtx-3G′′UG′U2Ux2-3G′′UG′USUUx+G′′USU2+G′U3Uxx-G′U2S′UUx+G′USUty+SUtt.Plugging this to the Taylor expansion (35), then the sliding average with the third order in time on a rectangular cell Iij will be approximated by(37)U¯i,jn+1=U¯i,jn-ΔtΔxΔy∫yj-1/2yj+1/2FUi+1/2,∗-FUi-1/2,∗dy-ΔtΔxΔy∫xi-1/2xi+1/2GU∗,j+1/2-GU∗,j-1/2dx+ΔtSU¯i,j+Δt2SUt¯i,j+Δt26SUtt¯i,j,with(38)FU=FU-Δt2F′U2Ux-F′USU+Δt263F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2S′UUx+F′USUt;GU=GU-Δt2G′U2Uy-G′USU+Δt263G′′UG′U2Ux2-3G′′UG′USUUx+G′′USU2+G′U3Uxx-G′U2S′UUx+G′USUt.F(U) and G(U) are approximated by numerical fluxes F^ and G^, for example, F(Ui+1/2,∗)≈F^i+1/2,∗=F^(Ui+1/2,∗L,Ui+1/2,∗R). In the fifth-order finite volume SWENO scheme, on rectangle control volume Ii,j,(39)-ΔtΔxΔy∫yj-1/2yj+1/2FUi+1/2,∗-FUi-1/2,∗dy-ΔtΔxΔy∫xi-1/2xi+1/2GU∗,j+1/2-GU∗,j-1/2dx.Three-point Gaussian quadrature formula is used for the above line integration. The numerical flux vectors F^(Ui±1/2,GL,Ui±1/2,GR) and G^(UG,j±1/2L,UG,j±1/2R) at Gaussian points (xG,yj±1/2), (xi±1/2,yG) are obtained by the regular fifth-order SWENO finite volume reconstruction and the Lax-Friedrich numerical flux as described in Section 2.2, where xG={xi-3/5Δx/2,xi,xi+3/5Δx/2} are the three Gaussian points. Meantime, reconstruction of Ux, Uxx and Uy, Uyy in the second and third term (38) for Gaussian points on cell interface will be calculated in this procession of SWENO reconstruction. That is, in this procession, we calculate UG,j+1/2, (Ux)G,j+1/2, (Uxx)G,j+1/2,…  in subroutine of SWENO reconstruction.In the reconstruction, the following Jacobin matrixes are needed,F′(U)=A(U), G′(U)=B(U), which have been mentioned above.(40)F′′U=A′U=000000000g+2u2D-2uD0-2uD2D00002uvD-vD-uD-vD01D-uD1D0,G′′U=B′U=0000000002uvD-vD-uD-vD01D-uD1D0g+2v2D0-2vD000-2vD02D,F′′′U=A′′U=A1,A2,A3,G′′′U=B′′U=B1,B2,B3,with(41)A1U=000000000-6u2D24uD204uD2-2D20000-6uvD22vD22uD22vD20-1D22uD2-1D20;A2U=0000000004uD2-2D20-2D2000002vD20-1D2000-1D200;A3U=0000000000000000002uD2-1D20-1D200000;B1U=000000000-6uvD22vD22uD22vD20-1D22uD2-1D20-6v2D204vD20004vD20-2D2;B2U=0000000002vD20-1D2000-1D200000000000;B3U=0000000002uD2-1D20-1D2000004vD20-2D2000-2D200.For source term, we can separate it into two parts; then the derivatives can be dealt with as scalar case as in 1D. (42)S¯i,j=0,-gζ¯bx¯,-gζ¯by¯i,jT=0,-gζ¯bx¯,0i,jT+0,0,-gζ¯by¯i,jT.Nine points of Gaussian surface integral are used for cell average source term.


## 2.1. The Shallow Water Equations

The 2D system of shallow water equations with conservative form is(1)Ut+FUx+GUy=SU,with conservational vector (2)U=D,Du,DvT,fluxes(3)FU=Du,Du2+gD22,DuvT,GU=Dv,Duv,Dv2+gD22T,and source terms(4)SU=0,-gDbx,-gDbyT,where D is the total water depth; u and v are the depth-averaged water velocities in the x- and y-directions, respectively; t is the time; g is the gravity constant; b is the topography function above; bx and by are derivatives with respect to x and y; they are used to describe the bed slopes stresses. Friction on the bottom, wind stress, and Coriolis term could also be considered as parts of source terms.We rewrite (1) in quasilinear form as(5)Ut+AUUx+BUUy=SU,where (6)AU=∂F∂U=010c2-u22u0-uvvu,BU=∂G∂U=001-uvvuc2-v202v.The eigenvalues of Jacobian matrixes A and B are(7)λA=u-c,u,u+c;λB=v-c,v,v+c.The normalized right and left eigenvectors of Jacobian matrixes A and B are as follows:(8)RA=101u-c0u+cv1v,LA=u+c2c-12c0-v01-u-c2c12c0;RB=101u1uv-c0v+c,LB=v+c2c0-12c-u10-v-c2c012c.We can easily check that RALA=I, RBLB=I, which is an important characteristic for hyperbolic conservation law. Neglecting variables in y-direction, we have the 1D conservative unsteady shallow water equations: (9)Dt+Dux=0,Dut+Du2+12gD2x=-gDbx.Similar way is that we can get Jacobian matrix and other information for 1D case.


## 2.2. The Fifth-Order Finite Volume SWENO Schemes

The fifth-order finite volume SWENO scheme will be used to numerically solve the spatial derivatives of shallow water equations [32, 37]. In this section, the basic procedure of SWENO will be described as short overview for the following 1D scalar hyperbolic conservation law:(10)ut+fux=0.DenoteIi=[xi-1/2,xi+1/2] as the ith cell, a control volume for finite volume method, centered on xi; we assume that the grid points are uniform with Δx=Δxi=xi+1/2-xi-1/2. For a finite volume scheme, we evolve the proximate average value u¯i of u¯(xi,t)=1/Δx∫Iiu(x,t)dx at mesh cell Ii in time. Using integral average, the conservation law (10) is rewritten into (11)ddtu¯xi,t=-1Δxfuxi+1/2,t-fuxi-1/2,t.f(u(xi+1/2,t)) is approximated by numerical flux f^i+1/2=f^(ui+1/2L,ui+1/2R). The popular Lax-Friedrich flux is used for its simpleness and monotone, which is given by (12)f^ui+1/2L,ui+1/2R=12fui+1/2L+fui+1/2R-αui+1/2R-ui+1/2L,where α=maxu|f′(u)|; ui+1/2L,R are approximations to u(xi+1/2,t) of left side and right side at xi+1/2 considering the discontinuous solution. We call this process reconstruction; here the reconstructions of interface values depend on the cell averages by the fifth-order SWENO reconstruction procedure.Onlyui+1/2L will be introduced in detail as follows. For the fifth-order SWENO scheme, three approximated point values of ui+1/2L should be constructed. Choose the first big stencil {Ii-2,Ii-1,Ii,Ii+1,Ii+2}; a fourth-degree polynomial p1(x) can be obtained by requiring the same mean value on the five cells. The first approximated point value p1(xi+1/2) for ui+1/2L is(13)ui+1/2L1=-120u¯i-2+920u¯i-1+4760u¯i-1360u¯i+1+130u¯i+2.Choose the second stencil {Ii-1,Ii}; another linear polynomial p2(x) can be obtained with the same idea. The second approximation p2(xi+1/2) is(14)ui+1/2L2=-12u¯i-1+32u¯i.Choose the third stencil {Ii,Ii+1}; a linear polynomial p3(x) can be obtained and the third approximation p3(xi+1/2) is(15)ui+1/2L3=12u¯i+12u¯i+1.Convex combination coefficients are the nonlinear weights. Before that we have to compute the smoothness indicators of the three above polynomials first. For the fifth-order SWENO reconstructions, the convex combination of all these three-point values is(16)ui+1/2L=ω11d1ui+1/2L1-d2d1ui+1/2L2-d3d1ui+1/2L3+ω2ui+1/2L2+ω3ui+1/2L3.The nonlinear weights ωm  (m=1,2,3) are defined by (17)ωm=αm∑s=13αs,αs=ds1+τε+βs,τ=β1-β2+β1-β322,and here ε is meaningless, just to avoid the denominator to be zero, so ε is a positive constant and typically ε=10-6, ds are the linear weights, and “smoothness indicators” βs are used to measure the smoothness of the three polynomials in cell Ii. The smaller the value of smooth indicator is, the smoother the polynomial is. We use the same smooth indicator formula for k-degree polynomial p(x) in [39](18)β=∑l=1k∫xi-1/2xi+1/2Δxi2l-1∂lpx∂xl2dx.The explicit expression of the three smoothness indicators βs is given by(19)β1=1144u¯i-2-8u¯i-1+8u¯i+1-u¯i+22+115600-11u¯i-2+174u¯i-1-326u¯i+174u¯i+1-11u¯i+22+7812880-u¯i-2+2u¯i-1-2u¯i+1+u¯i+22+14214611310400u¯i-2-4u¯i-1+6u¯i-4u¯i+1+u¯i+22;β2=u¯i-1-u¯i2;β3=u¯i-u¯i+12.For the fifth-order SWENO scheme on uniform mesh, the important idea is that the linear weights are your chosen positive constants; the only condition is that the summation is one. So the procedure of linear weights in classical WENO is avoided, and the negative weights are prevented. In this paper, we use the positive linear weights ds as (20)d1=0.98,d2=0.01,d3=0.01.On uniform mesh, the reconstruction ofui+1/2R is mirror symmetry on big stencil {Ii-1,Ii,Ii+1,Ii+2,Ii+3} to that of ui+1/2L. For shallow water equations, the reconstruction is operated in local eigenspace for more robust, so the local characteristic decomposition should be used, in which properties of eigenspace and Roe’s averages are needed.For 2D finite volume case, we cannot simply proceed in a dimension-by-dimension fashion as in finite difference method. The spatial surface integration on rectangular cellIi,j=[xi-1/2,xi+1/2]×[yj-1/2,yj+1/2] can be converted to line integration by Gauss theorem. In order to get the fifth order, we can use three-point Gaussian quadrature for line integration of every rectangular cell like in Figure 1. For example, to calculate numerical flux in x-direction, three Gaussian point values on every vertical line should be reconstructed from cell average values in y-direction first, then left and right Gaussian point values should be reconstructed by the SWENO scheme, then the numerical fluxes are calculated, and last the cell average value can be approximated on next time level.Figure 1
A sketch of SWENO reconstruction on a rectangular cell.


## 2.3. Prebalance of the Flux and Bottom Topography

We adopt the prebalance of the flux gradients and bottom topography in source term by Rogers et al. [38] using algebraic technique. Taking 1D shallow water equations (9), for example, the prebalanced formulation is(21)U~t+F~Ux=S~U,with (22)U~=ζ,DuT,F~=Du,Du2+gζ2+2hζ2T,S~U=0,-gζbxT.The main algebraic technique is to separate U into(23)U=Ueq+U~,with (24)Ueq=h,0T,U~=U-Ueq=ζ,huT,in which ζ is free surface elevation, the height of water from the still level h+b to free surface elevation; h is still water height, the height of water from bottom to still water level; D=ζ+h is the total water depth; see Figure 2 [36].Figure 2
A sketch of variables in shallow water equations.It is easy to verify that modified system (21) and system (9) have the same Jacobian matrix. Then use the finite volume SWENO discretization for (21).The specific value of the still water level depends on examples; we will describe the chosen datum below. It is perfectly reasonable to choose a fixed horizontal datum elsewhere and derive the prebalanced equations using a stage-discharge approach. Mostly, we choose the still water level between the maximum of bottom function and the minimum of free surface.


## 2.4. The Lax-Wendroff-Type Discretization for 1D Shallow Water Equations

We focus on the procedure of Lax-Wendroff-type time discretization scheme for 1D prebalanced shallow water equation (21) [33, 34]; for simplicity we still describe (21) as(25)Ut+FUx=SU.In this process, SWENO finite volume method will be used. Denote time level as tn+1=tn+Δt, Δt is time step, and U¯in is the approximate average value of U on Ii at the nth time level. The important idea in Lax-Wendroff-type time discretization procedure is to replace time derivatives by the spatial derivatives using Taylor expansion and PDE (25) [34]. By a temporal Taylor expansion with respect to time t for conservative vector U we obtain(26)Ux,t+Δt=Ux,t+ΔtUt+Δt22Utt+Δt36Uttt+Δt424Utttt+⋯,and the sliding average is(27)U¯x,t+Δt=U¯x,t+ΔtUt¯+Δt22Utt¯+Δt36Uttt¯+Δt424Utttt¯+⋯,where Ut¯ means mean value of time derivative Ut in x-direction, which will be given in detail as follows. In this paper we would like to obtain the third-order accuracy in time, so we need to approximate those terms on the right hand side in (27) until the third time derivatives and neglect other terms behind the third time derivatives. Theoretically, we can obtain any accuracy order in time and just need to approximate those terms until the according time derivatives.Subject to shallow water equation (25), by Lax-Wendroff type dime discretization, the time derivatives in (26) are turned to the spatial derivatives as follows: (28)Ut=-FUx+SU;Utt=-F′UUtx+SUt=F′U2Ux-F′USUx+S′U-F′UUx+SU;Utttt=-F′′UUt2+F′UUttx+SUtt=-3F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2SUUx+F′USUtx+-F′UUx+SU-S′′UF′UUx-S′UF′′UUx+S′′USU+S′U2-S′UF′U-F′′UUx2-F′UUxx+S′UUx.The abbreviated form is used like F′′(U)(Ut)2, which means F′′(U)UtUt; the details will be described later. By the way, the above vector source term is easy to be generalized to any dimensions; indeed, as the source term S(U)=[0,-gζbx]T, we only need to handle time derivative as scalar case in practice. For simplicity, we just keep time derivative of source term temporally as follows.So (26) can be written as(29)Ux,t+Δt=Ux,t-ΔtFU-Δt2F′U2Ux-F′USU+Δt263F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2SUUx+F′USUtx+ΔtSU+Δt2SUt+Δt26SUtt+⋯.Then the sliding average (27) with third order in time will be approximated by(30)U¯in+1=U¯in-ΔtΔxFUi+1/2-FUi-1/2+ΔtSU¯+Δt2SUt¯+Δt26SUtt¯,with(31)FU=FU-Δt2F′U2Ux-F′USU+Δt263F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2SUUx+F′USUt.F(U) is approximated by numerical flux F^, such as F(Ui+1/2)≈F^i+1/2=F^(Ui+1/2L,Ui+1/2R). With this form, the finite volume scheme has more flexibility in terms of mesh structure than that of finite difference method [36].Using the fifth-order SWENO described in Section2.2, we reconstruct point values at interface Ui+1/2L,R; meantime, we calculate point values of (Ux)i+1/2L,R and (Uxx)i+1/2L,R. We have to point out that only one order lower accuracy approximation for Ux is needed compared with that of U, as the extra Δt factor multiplied by the second term in flux F in (31). In the same way, we need the third accuracy order of approximation for Uxx, which is in the third term in (31). Indeed, for linear polynomial p2(x),p3(x), the first derivative is constant, the second derivative is zero, the nonlinear weights do not need to calculate for reconstruction of  Ux,Uxx, and only the same linear weight is used, which makes the process easier. For example,(32)Uxi+1/2L=p1′xi+1/2.For the fourth-degree polynomial p1(x), the computation is a little more complicated.In reconstruction ofF(U), we also need Jacobian matrix(33)F′U=AU=01c2-u22u,F′′U=A′U=00g+2u2h-2uh00-2uh2h,and, for simplicity, let F′′(U)=(F1,F2), where F1, F2 are the two 2×2 matrixes; F′′(U)(U′)2=F1U′,F2U′U′ is a 2×1 vector.In order to fit the scheme to problems with discontinuous or big slop in bottom topography, the bottom topography effect is included in numerical flux. In other words, we use SWENO reconstruction for bottom function and then use Gaussian quadrature for the source termsS¯i=[0,-gζ¯bx¯]iT, in which ζ¯ is known. This method can be generalized to varying bottom topography with time. In St¯i=[0,-gζt¯bx¯]iT and Stt¯i=[0,-gζtt¯bx¯]iT, ζt and ζtt have been calculated as part in flux F.For simplicity, we use SWENO5-LW3 to denote the scheme with the fifth-order finite volume SWENO scheme in space and the third-order Lax-Wendroff-type method in time and WENO5-RK3 to denote classical WENO scheme and the third-order Runge-Kutta methods for comparison.


## 2.5. The Lax-Wendroff-Type Discretization for 2D Shallow Water Equations

In 2D case, we need more information; this section aims to display the frame and idea. For 2D shallow water equation(34)Ut+FUx+GUy=SU,the sliding average of Taylor expansion is(35)U¯x,y,t+Δt=U¯x,y,t+ΔtUt¯x,y,t+Δt22Utt¯x,y,t+Δt36Uttt¯x,y,t+Δt424Utttt¯x,y,t+⋯.Using shallow water system (34), by Lax-Wendroff type dime discretization, the time derivatives in (35) are turned to the spatial derivatives as follows. For simplicity we keep time derivative of source term temporally.(36)Ut=-FUx-GUy+SU;Utt=-F′UUtx-G′UUty+SUt=F′U2Ux-F′USUx+G′U2Uy-G′USUy+SUt;Uttt=-F′′UUt2+F′UUttx-G′′UUt2+G′UUtty+SUtt=-3F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2S′UUx+F′USUtx-3G′′UG′U2Ux2-3G′′UG′USUUx+G′′USU2+G′U3Uxx-G′U2S′UUx+G′USUty+SUtt.Plugging this to the Taylor expansion (35), then the sliding average with the third order in time on a rectangular cell Iij will be approximated by(37)U¯i,jn+1=U¯i,jn-ΔtΔxΔy∫yj-1/2yj+1/2FUi+1/2,∗-FUi-1/2,∗dy-ΔtΔxΔy∫xi-1/2xi+1/2GU∗,j+1/2-GU∗,j-1/2dx+ΔtSU¯i,j+Δt2SUt¯i,j+Δt26SUtt¯i,j,with(38)FU=FU-Δt2F′U2Ux-F′USU+Δt263F′′UF′U2Ux2-3F′′UF′USUUx+F′′USU2+F′U3Uxx-F′U2S′UUx+F′USUt;GU=GU-Δt2G′U2Uy-G′USU+Δt263G′′UG′U2Ux2-3G′′UG′USUUx+G′′USU2+G′U3Uxx-G′U2S′UUx+G′USUt.F(U) and G(U) are approximated by numerical fluxes F^ and G^, for example, F(Ui+1/2,∗)≈F^i+1/2,∗=F^(Ui+1/2,∗L,Ui+1/2,∗R). In the fifth-order finite volume SWENO scheme, on rectangle control volume Ii,j,(39)-ΔtΔxΔy∫yj-1/2yj+1/2FUi+1/2,∗-FUi-1/2,∗dy-ΔtΔxΔy∫xi-1/2xi+1/2GU∗,j+1/2-GU∗,j-1/2dx.Three-point Gaussian quadrature formula is used for the above line integration. The numerical flux vectors F^(Ui±1/2,GL,Ui±1/2,GR) and G^(UG,j±1/2L,UG,j±1/2R) at Gaussian points (xG,yj±1/2), (xi±1/2,yG) are obtained by the regular fifth-order SWENO finite volume reconstruction and the Lax-Friedrich numerical flux as described in Section 2.2, where xG={xi-3/5Δx/2,xi,xi+3/5Δx/2} are the three Gaussian points. Meantime, reconstruction of Ux, Uxx and Uy, Uyy in the second and third term (38) for Gaussian points on cell interface will be calculated in this procession of SWENO reconstruction. That is, in this procession, we calculate UG,j+1/2, (Ux)G,j+1/2, (Uxx)G,j+1/2,…  in subroutine of SWENO reconstruction.In the reconstruction, the following Jacobin matrixes are needed,F′(U)=A(U), G′(U)=B(U), which have been mentioned above.(40)F′′U=A′U=000000000g+2u2D-2uD0-2uD2D00002uvD-vD-uD-vD01D-uD1D0,G′′U=B′U=0000000002uvD-vD-uD-vD01D-uD1D0g+2v2D0-2vD000-2vD02D,F′′′U=A′′U=A1,A2,A3,G′′′U=B′′U=B1,B2,B3,with(41)A1U=000000000-6u2D24uD204uD2-2D20000-6uvD22vD22uD22vD20-1D22uD2-1D20;A2U=0000000004uD2-2D20-2D2000002vD20-1D2000-1D200;A3U=0000000000000000002uD2-1D20-1D200000;B1U=000000000-6uvD22vD22uD22vD20-1D22uD2-1D20-6v2D204vD20004vD20-2D2;B2U=0000000002vD20-1D2000-1D200000000000;B3U=0000000002uD2-1D20-1D2000004vD20-2D2000-2D200.For source term, we can separate it into two parts; then the derivatives can be dealt with as scalar case as in 1D. (42)S¯i,j=0,-gζ¯bx¯,-gζ¯by¯i,jT=0,-gζ¯bx¯,0i,jT+0,0,-gζ¯by¯i,jT.Nine points of Gaussian surface integral are used for cell average source term.


## 3. Numerical Results

In this section, we present numerical results obtained by SWENO5-LW3 scheme and WENO5-RK3 scheme for a number of 1D and 2D examples for shallow water equations. We shall demonstrate the fifth-order convergence of the scheme in space and compare the numerical results by SWENO5-LW3 and WENO5-RK3, addressing CPU time and absolute truncation errors inL1 and L∞ norm. The parameter gravitation constant g=9.812m/s2. CFL condition number is 0.4 for SWENO5-LW3 and 0.6 for WENO-RK3. Uniform meshes are used.Example 1 (test forC-property).
In this test case [40, 41], the bottom function is (43)bx=5exp⁡-25x-52,x∈0,10. The initial conditions are (44)Dx,0=10-bx,ux,0=0.We choose the still water level at 10, which means still water height h(x,t)=10-b(x), and initial ζ(x,0)=0. The ending time is taken as t=0.5s, and continuous boundary condition is used. The surface should remain flat anytime.
In Figure3, we give the L1 and L∞ numerical truncation errors as a function of cost CPU time for water height D when cell numbers are 200, 400, 800, and 1500, and log scales are used; we can see that the two lines are crossed; we cannot say which method is absolutely better than the other one, but the accuracy is good and reaches around 10-14 with double precision support for each scheme, which means both methods can keep C-property very well. These conclusions are fit to numerical solutions of discharge Du.Figure 3
Numerical errors and CPU time for water depth in the 1DC-property test in Example 1. (a) L1 errors. (b) L∞ errors.
(a) (b)Example 2 (test for the accuracy order).
In scheme of the fifth order in space and the third order in time, in order to check accuracy order in space, the modified time stepΔtn=CFL·(Δx)5/3/maxi(|uin|+gDin) should be used to replace regular adaptive time step Δtn=CFL·Δx/maxi(|uin|+gDin). Of course, the smooth solution of the test is required.
We use the same test with smooth solution over a sinusoidal hump as Jiang and Shu [41] to check the order accuracy of the 1D scheme. In domain [0,1], the bottom function is (45)bx=sin22πx, and the initial data are given by (46)Dx,0=5+ecos⁡2πx,Dux,0=sin⁡cos2πx.
The still water level is taken at 5; in other words, still water heighth(x,t)=5-b(x); ζ(x,0) is then set uniquely. We compute until time t=0.1s; when the solution is still smooth, no shocks appear; periodic boundary conditions are adopted. Since the exact solution is unknown, numerical solution on 25600 cells is taken as reference.
In Table2, the L∞ errors and convergence orders by SWENO5-LW3 of water depth D and the discharge Du are listed when mesh is refined. L∞ errors and orders by WENO5-RK3 are also displayed in Table 2 for comparison. It can be seen that both schemes lead to the same fifth order of convergence, and SWENO5-LW3 produces more accuracy solutions than WENO5-RK3, the errors being an order of magnitude smaller. In Table 1, the L1 numerical errors and orders are shown, which have similar characteristics to L∞ measurement.
In Figures4 and 5, we give the L∞ and L1 errors against cost CPU time for water height D and discharge Du when cell numbers are 25∗2k-1,k=1,2,…,6. Log scale for both numerical errors and CPU time is used. In the two figures, the line of SWENO5-LW3 is always at the bottom relative to the line of WENO5-RK3, which means costing the same CPU time; SWENO5-LW3 can get higher accuracy. To get the same error, SWENO5-LW3 can save CPU time. For example, when cell number N=800, the cost CPU time of SWENO5-LW3 is 114.9 s, while 148.8 s for WENO5-RK3 and SWENO5-LW3 can save around 13% CPU time. To get the same error, SWENO5-LW3 can save more CPU time than 13%.Table 1
L 1 errors and numerical orders of accuracy for the one-dimensional order test in Example 2, SWENO5-LW3 and WENO5-RK3 using N cells.
N SWENO5-LW3 WENO5-RK3 D D u D D u L 1 error Order L 1 error Order L 1 error Order L 1 error Order 25 9.37 E - 03 0.00 5.95 E - 02 0.00 1.22 E - 02 0.00 1.14 E - 01 0.00 50 1.80 E - 03 2.38 1.80 E - 02 1.73 1.98 E - 03 2.62 1.99 E - 02 2.52 100 2.27 E - 04 2.98 1.93 E - 03 3.22 3.10 E - 04 2.67 2.68 E - 03 2.89 200 1.44 E - 05 3.98 1.23 E - 04 3.97 2.23 E - 05 3.80 1.92 E - 04 3.81 400 5.60 E - 07 4.69 4.79 E - 06 4.68 8.78 E - 07 4.67 7.56 E - 06 4.66 800 1.83 E - 08 4.94 1.56 E - 07 4.94 2.82 E - 08 4.96 2.43 E - 07 4.96Table 2
L ∞ errors and numerical orders of accuracy for the one-dimensional order test in Example 2, SWENO5-LW3 and WENO5-RK3 using N cells.
N SWENO5-LW3 WENO5-RK3 D D u D D u L ∞ error Order L ∞ error Order L ∞ error Order L ∞ error Order 25 4.57 E - 02 0.00 3.90 E - 01 0.00 6.38 E - 02 0.00 4.78 E - 01 0.00 50 1.32 E - 02 1.80 1.27 E - 01 1.62 1.76 E - 02 1.86 1.72 E - 01 1.47 100 3.48 E - 03 1.92 3.04 E - 02 2.06 4.62 E - 03 1.93 4.04 E - 02 2.09 200 3.54 E - 04 3.30 3.08 E - 03 3.30 5.35 E - 04 3.11 4.67 E - 03 3.11 400 1.83 E - 05 4.28 1.54 E - 04 4.32 2.88 E - 05 4.22 2.43 E - 04 4.26 800 6.38 E - 07 4.84 5.40 E - 06 4.84 1.01 E - 06 4.83 8.56 E - 06 4.83Figure 4
L ∞ errors and CPU time for the 1D order test in Example 2. (a) Water depth D. (b) discharge Du.
(a) (b)Figure 5
L 1 errors and CPU time for the 1D order test in Example 2. (a) Water depth D. (b) discharge Du.
(a) (b)Example 3 (dam breaking on a flatbed).
The stability of shock capturing must be tested and verified; the dam break problem is a proper test for shallow water flows [40] subjected to continuous boundary condition, and bottom topography is flat b(x)=0, so the source term disappeared. The initial conditions are taken as(47)Dx,0=D1x<0D2x≥0,Dux,0=0,with D1=1m and D2=0.1m. The computation is performed on [-1,1] until t=0.1s.
The numerical solutions using the SWENO5-LW3 and WENO5-RK3 on a mesh with 200 cells are shown in Figure6, which demonstrate good capturing of shock [42], and no spurious oscillations occurred for both schemes. For CPU time comparison, the value of CPU time is the total CPU time of running 200 times to reduce error; the total CPU time is 5.55 s for SWENO5-LW3 and 6.79 s for WENO5-RK3 scheme. The conclusion is that the cost CPU time of WENO5-RK3 is around 12% more than that of SWENO5-LW3 for this test.Figure 6
Dam break on a flatbed by SWENO5-LW3 and WENO5-RK3 in Example3. (a) Water depth D; (b) discharge Du.
(a) (b)Example 4 (dam breaking over a rectangular bump).
This test [41] is a generalization of above problem; discontinuous bottom topography is added to check the performance of scheme for great slope stress, or even bx does not exist. On the computational domain [0,1500], the Riemann bottom topography is given by (48)bx=8x-750≤150080otherwise, and the initial data are given by (49)Dx,0=20-bxx≤75015-bxotherwise,Dux,0=0. We take still water height h(x,t)=15-b(x), and still water level is 15. We compute up to t=15s and t=60s with continuous boundary condition. For the final time t=15s, the flow has not propagated to the rectangle bump, so the results are similar to above test. For the final time t=60s, we have to point out that the surface level is smooth, while the total water height D is discontinuous at the exact discontinuous points of bottom topography. The numerical results with 500 cells and reference results with 5000 cells are shown in Figure 7.
The following CPU time is run 20 time; the cost CPU time is 8.34 s for SWENO5-LW3 and 9.39 s for WENO-RK3; also around 13% CPU time is saved by SWENO5-LW3.Figure 7
The surface level of dam break on a discontinuous bed by SWENO5-LW3 and WENO5-RK3 in Example4. (a) t=15s; (b) t=60s.
(a) (b)Example 5 (a 2D accuracy test).
To check the accuracy order for 2D system of shallow water equations, we use the generalization of Example for the 1D accuracy test [41]. On the unit square [0,1]×[0,1], the bottom topography is (50)bx,y=sin⁡2πx+cos⁡2πy, and the initial data are given by (51)Dx,y,0=10+esin⁡2πxcos⁡2πy,Dux,y,0=sin⁡cos⁡2πxsin⁡2πy,Dvx,y,0=cos⁡2πxcos⁡sin⁡2πy.
We take the still water heighth(x,t)=10-b(x,y); 4-periodic boundary conditions are subjected in the problem. The time step is modified copying 1D case. Program is performed until time t=0.05s when the solution is still smooth. Numerical solution on 1600×1600 cells is taken as reference solution.
Similar numericalL1 and L∞ errors and orders are generated by the SWENO5-LW3 and WENO5-RK3; both schemes get the fifth convergence orders. Because the results are similar to that in [36], we omit the data. In Figure 8, we give the L∞ and L1 errors as functions of cost CPU time for D. One more time, the line of SWENO5-LW3 is at the bottom relative to the line of WENO5-RK3. When cell number N=400×400, the cost CPU time is 10242 s for SWENO5-LW3, while being 13465 s for WENO5-RK3; SWENO5-LW3 can save around 13% CPU time.Figure 8
Numerical error and CPU time of water heightD for the 2D order test in Example 5. (a) L1 errors. (b) L∞ errors.
(a) (b)Example 6 (a small perturbation of 2D steady-state water).
A classical problem [38] is used to check the capability of capturing the perturbation of the stationary state for 2D shallow water equations. The function for bottom topography is (52)bx,y=0.8exp⁡-5x-0.92-50y-0.52,x,y∈0,2×0,1. The initial conditions are (53)Dx,y,0=1-bx,y+0.010.05≤x≤0.151-bx,yotherwise,ux,y,0=vx,y,0=0.
The upper and lower boundaries are the boundaries of solid wall; the left and right boundaries are the extrapolation boundaries.
We take the still water level at 10. Figure9 shows the right-going disturbance as it propagates the hump; the 30 contours of the surface level D+b are showed at different time by SWENO5-LW3 on a mesh of 200×100 cells. The results by WENO5-RK3 are similar; we ignored these figures. These results demonstrate that the simulation of complex small features in water flow has high resolution in 2D case. On the same mesh, for final time 0.48 s, the cost CPU time for SWENO5-LW3 is 22.80 s separately, while for WENO5-RK3 it is 25.11 s. So the cost CPU time for WENO5-RK3 is around 11% more than that of SWENO5-LW3.Figure 9
A small perturbation of 2D steady-state water in Example6. 30 contours of water level at different time t=0.12s, 0.24s, 0.36s, and 0.48s.


## 4. Concluding Remarks

In the previous sections, we have studied the Lax-Wendroff-type SWENO scheme for numerical solution of nonhomogeneous hyperbolic conservation laws for shallow water equations with source terms. The scheme can simulate discontinuous flows correctly for prebalanced shallow water equations, while the scheme has more flexible mesh structure than finite difference method, less procedure of local feature decomposition, and simpler idea than classical WENO. We compare the numerical solutions by both SWENO5-LW3 and WENO5-RK3 addressing CPU time, numerical truncation error, and efficiency. The numerical solutions show that the SWENO5-LW3 scheme we established in previous section can simulate shallow water flow and catch the stronger shock. The SWENO5-LW3 scheme can save around 12% CPU costs while generating similar numerical error to that by WENO5-RK3. These conclusions agree with that of the finite difference WENO with Lax-Wendroff-type and Runge-Kutta time discretization for shallow water equations.



---

*Source: 2652367-2018-02-22.xml*