Moving mesh version of wave propagation algorithm based on augmented Riemann solver

Applied Mathematics and Computation 375 (2020) 125087

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Moving mesh version of wave propagation algorithm based on augmented Riemann solver Mina Bagherpoorfard a,∗, Ali R. Soheili b,c a

Department of Mathematics, Fasa Branch, Islamic Azad University, Fasa, Iran Department of Applied Mathematics, Ferdowsi University of Mashhad, Mashhad, Iran c The Center of Excellence on Modelling and Control Systems, Iran b

a r t i c l e

i n f o

Article history: Received 14 June 2019 Revised 1 November 2019 Accepted 25 January 2020

Keywords: Finite volume methods Moving mesh methods Riemann solvers Shallow water equations Wave propagation algorithm

a b s t r a c t Shallow water equations are a nonlinear system of hyperbolic conservation laws. This system admits discontinuous solutions. In addition, a discontinuous bathymetry can lead to the discontinuous solutions. The accuracy of numerical methods depends on the accuracy of the numerical solutions around these discontinuities. In order to increase the accuracy, the adaptive mesh methods are more cost-effective than the methods based on the global refinement, In terms of computational cost. Moving mesh method is an adaptive mesh method. In this method, an increase of the number of mesh points is not needed for increasing the accuracy and efficiency of numerical solutions. Highly accurate solutions only can be obtained by concentrating points in areas where more accuracy is needed and decreasing the number of points in smooth areas. In this paper, LeVeque wave propagation algorithm is improved by the moving mesh method, for shallow water equations with variable bathymetry. Moreover, augmented Riemann solver is used for solving Riemann problems in each time step. Finally, moving mesh version of wave propagation algorithm based on augmented Riemann solver is proposed. The numerical solutions based on the moving mesh method and fixed mesh method have a significant difference and application of a moving mesh method yields highly accurate solutions around shocks and discontinuities. © 2020 Elsevier Inc. All rights reserved.

1. Introduction Shallow water equations were first introduced by Saint Venant in 1871 for modelling of flow in a channel. These equations are based on conservation of mass and momentum laws. Moreover, shallow water equations can be estimated via 3D depth integration of Navier-Stokes equations and based on such hypothesis as fluid incompressibility and hydrostatic pressure regardless of vertical acceleration of particles [51,64,65]. These equations are utilized for those flows in which flow wavelength is greater than flow depth, such as tsunami [58]. Furthermore, a wide spectrum of hydraulic engineering problems and coastal phenomena, including tide flows, rivers [11,32], floods [12,21,83], dam break [2,81], propagation of tsunami waves [50,58,71,83], reservoirs, and channel flows [66,67], can be efficiently simulated by these set of equations. Shallow water equations are a nonlinear system of hyperbolic conservation laws. Hyperbolic characteristics of these equations make their numerical and analytical solutions more complicated. Special and simple analytical solutions are accessible ∗

Corresponding author. E-mail addresses: [email protected] (M. Bagherpoorfard), [email protected] (A.R. Soheili).

https://doi.org/10.1016/j.amc.2020.125087 0 096-30 03/© 2020 Elsevier Inc. All rights reserved.

2

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

in papers of other sciences [6]. Delstre et al. (2013) collected a set of accessible analytical solutions [22]. However, analytical solutions in real simulations are inaccessible, and consequently, strong, efficient, and trustable numerical schemes are of great importance. Some of the finite difference and finite element methods can be observed in [8,10,16,61]. Among the literature of shallow water equations, different types of finite volume methods can be observed. In [19,68], high order finite volume methods (5 or 6) have been analyzed in order to acquire numerical solutions with high accuracy. Felcman and Havle have obtained numerical results using a vijayasundaram finite volume scheme [26]. In some of the works, staggered finite volume discretization is used instead of collocated discretization for nonlinear hyperbolic conservation laws [9,33,55,77,79]. In this method, different unknowns of the system are approximated on staggered meshes. See [23,34,36]. Godunov-type methods are significant types of finite volume methods [31]. These methods are successful in challenge with most of the hyperbolic conservation complications [3,4,30,59,63,78,80,87]. Other numerical methods can be extracted from Kärnä et al. [47], Kesserwani and Liang [49], Vater et al. [82], Xing and Shu [85]. An efficient numerical scheme for shallow water equations should be well-balanced. In most of the numerical methods, it is important for the presented scheme to be in a steady state [7,14,15,18,20,29,48,60,69,70]. Also, simulation of real problems with shallow water equations imposes the source terms on the equations system without which, description of these phenomena will doom to failure. Most of the source terms are related to topography and friction force. Source terms create some complications for numerical approximations; therefore, expansion of numerical method which can encounter with these complications have always been taken into consideration [18,45,53,55,86]. George has introduced an augmented Riemann solver for Levegue wave propagation algorithm [54] of shallow water equations which is also applicable with source term [30]. This solver includes appropriate characteristics of f-wave method, Roe and HLLE Riemann solvers. Furthermore, some tricks are presented for specific conditions, such as strong rarefaction waves, dry state problems, resonance problems, and coastal phenomena. Usually, in real solutions of shallow water equations shock or discontinuity can be observed in areas of the physical domain due to hyperbolic characteristics. Moreover, a discontinuous bed leads to discontinuity of solutions. Numerical analysis of these areas require a very refined grid in order to reach an acceptable accuracy. Using global grid refinement is not considered to be efficient. Therefore, using efficient adaptive mesh methods is necessary for increasing accuracy and decreasing the computational costs. One of the main groups of adaptive mesh methods is moving mesh method. In moving mesh methods, a monitor function, which is used for identification of areas with rapid changes of solutions, makes the points more concentrated around these areas, and therefore, density of points can be reduced around smooth areas of solutions. For time-dependent differential equations, moving mesh methods are divided into two static [46,62,75,76] and dynamic groups [28,37,41,43,44,74]. In dynamic groups, a time dependent mesh equation, including velocity of mesh points, is used for continuous movement of mesh points [13,17,38–40,42]. This moving mesh partial differential equation (MMPDE) is coupled with physical partial differential equations, therefore, physical solutions of the problem and geometrical location of the gird point can be obtained simultaneously. However, in static methods, numerical solutions of partial differential equations and distribution of new grid points are two separate processes. An initial mesh and approximated solutions of physical equations are accessible in each time step. Afterwards, a new mesh is built upon the structure of numerical solutions, and finally, approximate solutions on the new mesh are obtained by interpolating on the old mesh. In this paper, we equip augmented Riemann solver with moving mesh method. We expect that our new method will be more efficient for shallow water equations containing shock or discontinuity in their solutions. 2. Shallow water equations and Riemann problem Shallow water equations are a nonlinear system of hyperbolic conservation laws

∂h ∂ + (hu ) = 0, ∂t ∂ x ∂ ∂  2 1 2 (hu ) + hu + gh = −ghbx , ∂t ∂x 2

(2.1a) (2.1b)

where h(x, t) stands for fluid depth, u(x, t) for horizontal velocity of the fluid, b(x) for bottom elevation, z(x, t ) = h(x, t ) + b(x ) for moving free surface, and g for gravity. If we put

 

q(x, t ) =

 

h q = 1 , hu q2





hu f (q ) = = hu2 + 12 gh2 and





0 ψ (x, t ) = , −ghbx





q2 q22 , 1 + gq21 q1 2

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

3

system (2.1) can be written as the following conservation law:

qt + f (q )x = ψ (q, x )

x ∈ R, t > 0.

Jacobian matrix of f(q) is



f  (q ) =

0

1

−( qq21 )2 + gq1

2q2 q1

that includes eigenvectors



r 1 (q ) =

u−

1

 ,



 ,

gh

r 2 (q ) =

and corresponding eigenvalues

λ1 (q ) = u −



λ2 (q ) = u +

gh

u+



1

 gh

,

gh.

Weak solutions of hyperbolic conservation laws which are created in physical problems should be applied around discontinuities and shocks in Rankine–Hugoniot condition [56]

s[[q]] = [[ f (q )]],

(2.2)

in which s is the discontinuity propagation speed, and [[.]] is the jump in the desired quantity. Uniqueness of the weak solutions is not approved and more entropy conditions are required for validation of physical solutions [53,78]. Generally, a Riemann problem is the hyperbolic conservation law along with the constant piecewise initial conditions as following:

qt + f (q )x = ψ (q, x )

q(x, 0 ) =

ql qr

; x ∈ R, t > 0,

f or x < 0, f or x > 0.

(2.3a)

(2.3b)

In nonlinear Riemann problems, the solution includes a set of transitions in each characteristic family which connect ql to qr . These transitions separate the areas with constant solutions. Moreover, they may be shock, rarefaction or contact discontinuities. Both characteristic families of shallow water equations are genuinely nonlinear and Riemann problem for these equations includes two transitions. Both transitions are shock wave or rarefaction wave [56]. These two waves separate the single constant middle state q∗ from ql and qr . q∗ in term of the type of characteristic in each family is measurable based on the Hugoniot locus or integral curves [56]. For example, if Riemann problem includes two shocks in both characteristic families, q∗ should be on Hugoniot locus passing through ql and through qr , then

⎧

⎨u∗ = ul − (h∗ − hl ) g 1∗ + 1 , hr

2 h g 1 ⎩u∗ = ur + (h∗ − hr ) + h1r . 2 h∗

By removing u∗ , a single nonlinear equation will be achieved for h∗ . Numerical methods which are dependent on solving Riemann problems are so costly in comparison with other methods, similar to Godunov-type finite volume methods. In these methods, one Riemann problem should be solved in each cell interface in order to update the solutions in each time step. Therefore, an exact or approximate Riemann solver which produces efficient enough solutions is of great importance. There are a wide spectrum of approximate Riemann solvers which have low costs in comparison with exact Riemann solvers. Moreover, these solvers produce efficiently accurate approximate solutions, such as Roe [72] and HLLE [24,25,35] approximate Riemann solvers. 3. The wave propagation algorithm Godunov type methods is a kind of finite volume methods for conservation laws which yield consistent and stable approximations for numerical fluxes. In these methods, in each time step, the Riemann problem with the following initial data should be solved in each cell interface xi− 1 :



q(x, t n ) = where

Qin ≈

1 x

2

Qin−1 Qin

 Ci

f or x < xi− 1 , 2 f or x > xi− 1 ,

q(x, t n )dx,

2

(3.1)

4

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

is the average value of exact solution on grid cell Ci = [xi− 1 , xi+ 1 ] with the length x = xi+ 1 − xi− 1 . Wave propagation algo2

2

2

2

rithm is a Godunov-type finite volume method introduced by Levegue [54,56]. This algorithm is applicable for all hyperbolic systems

qt + A(q, x )qx = 0,

q ∈ Rm ,

A(q, x ) ∈ Rm×n .

In this method, the structure of Riemann solution in each cell interface is used for updating numerical solutions. Solutions p p of Riemann problem (3.1) in xi− 1 produce a collection of waves w 1 with speeds s 1 , p = 1, . . . , Mw , so that i− 2

2

Qin − Qin−1 =

Mw  p=1

i− 2

wip− 1 .

(3.2)

2

High resolution wave propagation formula is

Qin+1 = Qin −

 t   t  + n ˜n 1 , A Qin− 1 + A− Qin+ 1 + F˜i+ 1 + F i− 2 2 2 2 x x

in which, fluctuations A± Q n

A+ Qin− 1 = 2

A− Qin− 1 = 2

i− 12

 { p:s

p i− 1 2

>0}

 { p:s p

i−

<0} 1

are defined as

sip− 1 wip− 1 ,

(3.4a)

sip− 1 wip− 1 .

(3.4b)

2

2

2

2

2

and, the modified term F˜ n

i− 12

is defined based on the first order waves [52,56]:



n F˜i− 1 = 2



Mw 1 t p |sip− 1 | 1 − |s | w˜ ip− 1 , 2 x i− 12 2 2 p=1

˜p 1 w i− 2

in which

(3.3)

= φw

p i− 12

is the limited version of wi−1/2 , and φ is a TVD high resolution limiter such as MC limiter or p

minmod limiter [56,73]. If this algorithm is used for conservation law problems, the following condition should be satisfied in order to preserve the conservative numerical solution:

f (Qi ) − f (Qi−1 ) = A− Qin− 1 + A+ Qin− 1 . 2

(3.5)

2

In order to satisfy numerical conservation condition (3.5), instead of factorization (3.2), flux jump can be split in to a set of vectors:

f (Qi ) − f (Qi−1 ) =

Mw  p=1

Z

p , i− 12

Zip− 1 .

(3.6)

2

p = 1, . . . , Mw , are called f-waves [5], and s

p i− 12

are corresponding wave propagation speeds. In this method, the fluc-

tuations are defined based on the f-waves, as (3.4). If for some of the p, Z

p i− 12

are non-zero and static, these waves should

be used in computations of fluctuations to maintain conservative numerical solutions. There is no guarantee for positivity of middle depths q∗l and q∗r , since jumps in solutions are not used in this factorization. For non-homogeneous conservation laws,f-wave method is still applicable. The reader can refer to [5,56] for further details. Also, the modified term F˜ n 1 for i− 2

f-wave method is n F˜i− 1 = 2





  Mw 1 t p sgn sip− 1 1− |s | Z˜ip− 1 , 2 x i− 12 2 2 p=1

p in which Z˜

i− 12

= φZ

p i− 12

is a limited wave.

4. Augmented Riemann solver George has introduced a comprehensive Riemann solver for wave propagation algorithm [30] that includes strengths of f-wave method, and standard Riemann solvers, such as Roe and HLLE. Moreover, it can produce physically valid numerical solutions for other special states, such as resonant problem, dry-state problem, and strong rarefaction.

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

5

4.1. Structure of augmented Riemann solver Augmented Riemann solver is based on the splitting the vector of initial Riemann data

⎡

⎤

hi − hi−1

⎢ hu − hu ⎥  4 i i−1 ⎢ ⎥ α p 1 wp 1 ⎢ ⎥= ⎣φ (Qi ) − φ (Qi−1 )⎦ p=1 i− 2 i− 2

(4.1)

bi − hi−1 in which φ (q ) = hu2 + 12 gh2 , bi is the approximate average value of bottom elevation on ci , and w a set of waves with corresponding speeds s

⎡ ⎤

⎡

1

w1i− 1 2

1

⎤

p . i− 12

p i− 12

∈ R4 , p = 1, . . . , 4, are

These waves are defined as

⎡ ⎤ 0

⎢s1 ⎥ ⎢ s2 ⎥ ⎢0 ⎥ ⎢ ε⎥ ⎢ ε ⎥ ⎢ ⎥ = ⎢ ⎥, w2i− 1 = ⎢ , w1i− 1 = ⎢ ⎥, ⎥ 2 2 2 2 2 ⎣sε ⎦ ⎣ ( sε ) ⎦ ⎣1⎦ 0

(4.2)

0

0

⎡

w4i− 1 2

⎤  1 2 s+ − (Qi−1 , Qi )/ − λ λ (Qi−1 , Qi ) ⎢ ⎥ 0 ⎢ ⎥ =⎢ ⎥, + ⎣ ⎦ s− (Qi−1 , Qi ) 1

and

s1i− 1 = s1ε , 2

s3 + s2ε s1i− 1 = ε , 2 2

s2i− 1 = s2ε , 2

In the third component of w4

i− 12

s4ε = 0.

(4.3)

,

¯ s+ − (Qi−1 , Qi ) = −gh (Qi−1 , Qi )

1 λ2 (Q λ i−1 , Qi ) , λ1 λ2 (Qi−1 , Qi )

in which

hi−1 + hi h¯ (Qi−1 , Qi ) = . 2 1 2 1 λ2 (Q 1 2 If we consider λ1 (q )λ2 (q ) = u2 − gh, average values λ i−1 , Qi ) and λ λ (Qi−1 , Qi ), as approximations of λ (q)λ (q), are defined as

λ1 λ2 (Qi−1 , Qi ) = u¯ 2 − gh¯ , 1 λ2 (Q ¯ λ i−1 , Qi ) = ui−1 ui − gh,

where

u¯ =

ui−1 + ui . 2

In order to protect conservative numerical solutions, in definition of Z are used:



Zip− 1 = 02×1 2

I2×2



02×1 ,

p i− 12

∈ R2 , two middle components of α

p p w i− 12 i− 12

∈ R4

(4.4)

in which 0 is the zero matrix, and I is the identity matrix [30]. 4.2. Some properties of augmented Riemann solver Shock preservation. In shallow water equations, if a shock with propagation speed sp exists in pth characteristic family, then p ˆ p (Qi−1 , Qi ). Therefore, augmented Riemann solver, such as the Roe solver, produces based on Lax entropy condition, sε = λ the real solutions [72].

6

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

Fig. 1. Strong rarefaction in the first characteristic family, s1ε = λ1 (Ql ) << λ1 (q∗ ).

Maintaining depth positivity. For shallow water equations, when there is a variable bed, α4 = bi − bi−1 = 0 which yields two effective middle states h∗l and h∗r . Despite f-wave method which has no guarantee for positivity of h∗l and h∗r , splitting (4.1) provides an opportunity for effective middle state control. For example, in a subsonic problem wich s1ε < 0 < s2ε , using the following bound for (W4 )1 maintaines positivity of effective middle states:

s2ε − s1ε h∗ ≤ (w4 )1 ≤ −1 1 sε bi − bi−1

i f bi − bi−1 > 0,

s2ε − s1ε h∗ ≥ (w4 )1 ≤ −1 s2ε bi − bi−1

i f bi − bi−1 < 0,

where h∗ is HLLE middle state depth. Such a similar bound exists for supersonic problems. Increasing accuracy for the strong rarefaction. If there is a strong rarefaction in kth characteristic family, the selection

⎡

1

⎤

⎢ λk (q∗ ) ⎥ ⎢ ⎥ ⎥, ⎣(λk (q∗ ))2 ⎦

Wi3− 1 = ⎢ 2

0 creates a Riemann approximation with higher accuracy. For example, If the strong rarefaction is in the first characteristic family

s1ε = λ1 (Qi−1 ) << λ1 (q∗ ), then, based on 1-Riemann invariant for middle state we have

  λ1 (q∗ ) = ui−1 + 2 ghi−1 − 3 gh∗ ,

where λ1 (q∗ ) corresponds to the leading edge of the rarefaction (see Fig. 1). Einfeldt speeds correction for dry-state problems. When we have an initial dry-state, solution of real Riemann problem merely includes one rarefaction wave. Leading edge of this rarefaction wave, or wet-dry interface, moves with the following speeds:

 λ1 (q∗− ) = ui−1 + 2 ghi−1 i f hi = 0,

(4.5a)

 λ2 (q∗+ ) = ui − 2 ghi

(4.5b)

i f hi − 1 = 0.

then, the real single wave in a family expands more rapidly than Einfeldt speed in another family. This issue is not valid from physical point of view, since in a real Riemann problem, the wave speed (shock or rarefaction) in a family does not exceed the wave speed in the other family [78] (see Fig. 2). The new Einfeldt speeds in dry-state problem leads to

s˜1ε = λ2 (q∗ ), s˜1ε = s1ε ,

s˜2ε = s2ε if hi−1 = 0,

s˜2ε = λ1 (q∗ ) if hi = 0.

that the inaccuracy is solved via this definition.

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

7

Fig. 2. Insufficient speed estimates for near-dry-state Riemann problems. hr ≈ 0, s2ε = λˆ2 < λ1 (q∗ )(le f t ) and hl ≈ 0, s1ε = λˆ1 > λ2 (q∗ ) (right).

Fig. 3. Near-dry-state:

rarestrength rarestrength rarestrength ≈ 1, weak rarefaction: ≈ 0, and strong rarefaction: α1 < < α2 , with 0 < α 1 < α 2 < 1. s2ε − s1ε s2ε − s1ε s2ε − s1ε

We use the following criteria in our computations for Riemann problems which include a refraction wave in one of their characteristic families. This criteria is used to determine initial dry-state or strong refraction wave. Experimental coefficients α 1 and α 2 are selected, so that 0 < α 1 < α 2 < 1, α 1 is close to zero, and α 2 is close to 1. We put

αst =

strength of the rarefaction , s2ε − s1ε

where the strength of the rarefaction in the first and the second family are defined as

  λ1 (q∗ ) − λ1 (Qi−1 ) = 3( ghi−1 − gh∗ ),   λ2 (qi ) − λ2 (q∗ ) = 3( ghi − gh∗ ),

respectively. If α st > α 1 , there is a Riemann problem near to dry-state, if α st < α 1 , there is a weak refraction wave, if α 1 < α st < α 2 , there is a strong refraction wave. This is indicated in Fig. 3. For more details about augmented Riemann solver, such as well balanced properties, challenge with sonic problem, and etc, the [30] can be studied. 5. The moving mesh method Suppose that [a,b] is the physical domain with variable x, and [0,1] is the computational domain for a computational variable ξ . A coordinate transformation between these domains is

c = [0, 1] → p = [0, b], x = x ( ξ , t ), ξ ∈ [0, 1], x(0, t ) = a, x(1, t ) = b. Let ξ j = Nj , j = 0, 1, . . . , N, is a uniform mesh on the computational domain. Then, a = x0 < x1 < · · · < xN = b, is a corresponding mesh on the physical domain where x j = x j (ξ , t ), j = 0, 1, . . . , N, is obtained based on the equidistribution principle(EP) [39], so that some measure of solution error is equal on each subinterval. We introduce a monitor function M(x, t) which indicates some error measurement in the approximate solutions of original PDE. The mesh in x satisfies in the EP for

8

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

all values of time t that could be written down as

 x(ξ ,t ) a

where

θ (t ) =

M (x˜, t )dx˜ = ξ θ (t ), 

b a

(5.1)

M (x˜, t )dx˜.

Twice differentiations from (5.1) results in the equdistribution principle

∂ ∂ [M ( x ( ξ , t ), t ) x ( ξ , t )] = 0 , ∂ξ ∂ξ

(5.2)

which is called moving mesh equation. 5.1. Static moving mesh method In static methods [46,62,75,76], the process of obtaining numerical solutions in each time step depends on two separate sections: a redistribution algorithm for mesh points, which is an iterative process, and a numerical algorithm for solving the original PDE, which should be an efficient method compatible with physical problem. One of the most important characteristics of static methods is updating the solutions in an iterative process. In each iteration, the numerical solutions are interpolated from the old mesh into the new mesh. For example, in [76], a conservative interpolation has been presented for conservation laws. Its main idea is preservation of mass conservation law for numerical solutions in each step of mesh redistribution. We have used this interpolation in our computations. The new mesh should be satisfied in the equdistribution principle (5.2) in each iteration [75]. By considering artificial time τ , Eq. (5.2) changes to

xτ = (Mxξ )ξ ,

(5.3)

with the boundary conditions x(0, τ ) = a and x(1, τ ) = b. By discretization of Eq. (5.3), the following moving mesh equation can be obtained:

x˜j = x j + in which ξ =

 τ  M ( u 1 )(x j+1 − x j ) − M (u 1 )(x j − x j−1 ) , 1 ≤ j ≤ N − 1, j+ j− 2 2 ξ 2

1 N,

(5.4)

is the step length in c , and u j+ 1 is the average value of PDE solution on [x j , x j+1 ]. 2

The following algorithm can be used in static moving mesh method. In dynamic moving mesh methods, a moving mesh partial differential equation (MMPDE) [17,42] is required. This MMPDE can be coupled with original PDE and then physical solutions and location of grid points can be obtained simultaneously [39]. For further details about dynamic moving mesh method, MMPDEs and their theories, refer to [38–40]. 5.2. Moving mesh version of wave propagation algorithm based on augmented Riemann solver Augmented Riemann solver, which was introduced in Section 4, has appropriate properties. In most cases, shallow water equations include shock or discontinuity in areas of the physical domain due to some of the inherent hyperbolic characteristics. For these problems, the moving mesh techniques are more efficient than using a fixed uniform mesh. In this paper, we have improved augmented Riemann solver for Leveque wave propagation algorithm by moving mesh with the aim of increasing accuracy and efficiency of the numerical solutions. For this end, Algorithm 1 along with the following items are used: Algorithm 1 Static Moving Mesh. [0]

Step 1: Define a uniform mesh on domain c . Obtain the initial mesh x j = x j on p using moving mesh equation (5.3), and then compute u

[0] j+ 12

values based on the initial data.

Step 2: Do the following for ν = 0, 1, · · · [ν ]

a) Move the mesh points x j b) Interpolating the

[ν +1] u 1 j+ 2

[ν +1]

to x j

, using equidistribution principle (5.3).

values on the new mesh.

c) Repeat updating procedures (a) and (b) in a fixed number of repetitions or until x[ν +1] − x[ν ] < ε . [ν +1]

+1 Step 3: Obtain the solution of the original PDE unj+1 on x j

Step 4: If tn+1 < T , put

[0] u 1 j+ 2

=

un+11 , j+ 2

[0] xj

=

[ν +1] xj ,

, with an appropriate numerical algorithm.

and return to the step 2.

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

9

Fig. 4. Mesh trajectories for problem 6.1, dam break over a step bed, for t ∈ [0, 0.4].

-In the second stage of step 2, the interpolation introduced in [75,76] is used in order to update numerical solutions on the new mesh. -In the third step, high resolution wave propagation algorithm (3.3) is applied for numerical solutions of shallow water equations. -Augmented decomposition (4.1), and (4.4) are utilized for preparing the f-waves. In the next section, numerical results of multiple test problems are presented based on this technique. 6. Numerical results and discussion In the previous sections, a moving mesh version of the Riemann solver was presented for shallow water equations with the source terms. For shallow water equations with the shock or discontinuity in their solutions, or with the static discontinuous bathymetry, the moving mesh methods are expected to be more successful than the fixed mesh scheme. In Sections 6.1–6.3, multiple test problems are analyzed with the aim of more accurate and more efficient presentation of moving mesh method. Also, the CPU time is compared between the moving mesh method and fixed mesh method at the level of accuracy. In Section 6.4, a test problem is presented to show that the new technique preserves the well-balanced property. The semi arclength monitor function has been used in our computations:



M = 1 + α1 |hx | + α2 |ux | + α3 |hux | + α4 |zx |



2

(6.1)

6.1. Dam break over a step bed This test is a Riemann problem on a channel with the length of 12. The initial conditions are



h(x, 0 ) =

5 1

if x ≤ 6, , if x > 6

u(x, 0 ) = 0,

with a step bed



b( x ) =

0 1

if x ≤ 6, , if x ≥ 6

and final time is T = 0.4 [1]. We have considered a mesh with NR = 500 as the reference mesh and its corresponding numerical solutions as the reference solutions. Then, the numerical solutions with n = 100 have been compared with the reference solutions. In (6.1), α1 = α2 = α3 = 10, α4 = 0 have been used. Moreover, step 2 of Algorithm 1 is repeated 10 times for redistribution of the mesh points and updating the solutions on the new mesh. In Fig. 4, mesh trajectories have been drawn on the physical domain of the problem for t ∈ [0, 0.4]. According to this figure, concentration of the mesh points is around the rightward going shock, therefore, in T = 0.4, they are more concentrated around x = 8.3. Moreover, in all the times, a dense area of mesh points is around the static discontinuity bed. Numerical solutions z(x, 0.4) and u(x, 0.4) based on the moving mesh scheme and fixed mesh scheme have been indicated in Figs. 5, 6. In the areas with smooth solutions, accuracy of both schemes are the same. However, in the areas with the rapid changes in the solutions, the moving mesh technique leads to more accurate numerical solutions. 8 percent of the grid points have been concentrated in neighborhood of radius 0.2 around the shock and 10 percent of them have been located in neighborhoods of radius 0.3 around the discontinuity. To obtain this level of accuracy by the fixed mesh method, a grid

10

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

Fig. 5. Numerical solutions z(x, 0.4 ) = h (x, 0.4 ) + b(x ) using two methods, moving mesh and fixed mesh with n = 100, and comparison them with reference solution.

Fig. 6. Numerical solutions u(x, 0.4) using two methods, moving mesh and fixed mesh with n = 100, and comparison them with reference solution. Table 1 Comparison CPU time for moving mesh and fixed mesh methods. Test problem

CPU time (sec) moving mesh

CPU time (sec) fixed mesh

Dam Break over a Step Vacuum over Discontinuous Topography Flow over Non-Horizontal Bed

127.49 130.94 33.44

267.49 317.48 96.39

with at least n = 240 must be used. The comparison of the CPU times of these two methods, with the same accuracy, has been presented in Table 1. The results indicate that the moving mesh is more effective in terms of the computational cost. CF L = 0.4 has been applied in these computations. 6.2. Vacuum over discontinuous topography This problem includes incipient cavitation on a discontinuous bed [9,30]. The initial conditions are

h(x, 0 ) = 10,



hu(x, 0 ) =

−350 350

if x < o.w.

50 , 3 ,

with the final time T = 0.25. The topography includes two discontinuities as the following:



b( x ) =

1 0

if 25 < x < 12.5, 3 . o.w.

For this problem, NR = 800, α1 = α4 = 10, α2 = α3 = 0, and in the numerical algorithm, ν = 10, have been applied. Then, the numerical solutions have been computed based on the moving and fixed mesh schemes with n = 150 and CF L = 0.8. In Fig. 7, the reference solution h(x, 0.25) has been presented. This solution includes two discontinuities around x = 25 3 and x = 12.5; furthermore, relatively rapid variations in solution are observed around x = 1.5. In order to demonstrate

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

11

Fig. 7. Reference solution z(x, 0.25 ) = h (x, 0.25 ) + b(x ) with NR = 800.

Fig. 8. Closed figure of numerical solutions z(x, 0.25 ) = h (x, 0.25 ) + b(x ) for x ∈ [8, 9] using two methods, moving mesh and fixed mesh with n = 150, and comparison them with reference solution.

Fig. 9. Closed figure of numerical solutions z(x, 0.25 ) = h (x, 0.25 ) + b(x ) for x ∈ [12, 13] using two methods, moving mesh and fixed mesh with n = 150, and comparison them with reference solution.

efficiency of moving mesh method, closed figures around these three areas have been plotted in Figs. 8–10. Mesh trajectories for t ∈ [0, 0.25] can be seen in Fig. 11. These mesh trajectories indicate higher concentration of points around the static discontinuities and those areas with rapid varieties in solution. This behavior is efficiently consistent with the physical problem. The neighborhoods of radius 0.2 around both main discontinuities contain approximately 2 percent of the grid points. A uniform grid with at least n = 375 can provide such accuracy around the discontinuities but it increases the CPU time, significantly (see Table 1).

12

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

Fig. 10. Closed figure of numerical solutions z(x, 0.25 ) = h (x, 0.25 ) + b(x ) for x ∈ [1.2, 3] using two methods, moving mesh and fixed mesh with n = 150, and comparison them with reference solution.

Fig. 11. Mesh trajectories for problem 6.2, vacuum over discontinuous topography, for t ∈ [0, 0.25].

6.3. Flow over non-horizontal bed The third test problem is defined on a non-horizontal bed



b( x ) =

0 0.1(x − 10 ) 1

if 0 ≤ x ≤ 10, if 10 < x ≤ 20, if 20 ≤ x ≤ 30

,

for x ∈ [0, 30] with the initial conditions

u(x, 0 ) = 0, z(x, 0 ) = h(x, 0 ) + b(x ) =



4 2

if 0 ≤ x < 5, , if 5 ≤ x ≤ 30

for more details see [70]. Reference solutions with NR = 600, and numerical solutions h(x, t) and u(x, t) based on the moving and fixed mesh techniques for n = 120, have been obtained at T = 1. Diagrams of numerical solutions h(x, 1) and u(x, 1), and also mesh trajectories diagram for t ∈ [0, 1], have been plotted in Figs. 12–14, respectively. After x = 12.5, numerical solutions h(x, 1) and u(x, 1) are constant. Therefore, distribution of mesh points is nearly uniform. In Fig. 14, a high density area is located around x = 11 which corresponds to rapid variations in u(x, 1) and h(x, 1) in this area. High concentration of points in this area causes high accuracy of numerical solutions. As can be seen in Figs. 12, 13, the numerical solutions based on the moving mesh method are much better than the approximations which are done with fixed mesh, and the difference between these two types of solutions is significant. Furthermore, there are the solutions with middle variations in [0,2] that the numerical solutions with moving mesh are still more successful. Eighteen points equivalent to 15 percent of grid points are placed in the neighborhood of radius 1 around the shock. The fixed mesh method with n = 270 produces the same results around the shock but, according to Table 1, imposes more computational cost. In computations of this problem α1 = α2 = 15, ν = 5, and CF L = 0.5.

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

Fig. 12. Numerical solutions h(x, 1) using two methods, moving mesh and fixed mesh with n = 120, and comparison them with reference solution.

Fig. 13. Numerical solutions u(x, 1) using two methods, moving mesh and fixed mesh with n = 150, and comparison them with reference solution.

Fig. 14. Mesh trajectories for problem 6.3, flow over non-horizontal bed, for t ∈ [0, 1].

6.4. Steady flow over a bump The last test problem is selected from the transcritical flow without shock class. The bump is defined as



b( x ) =

0.2 − 0.05(x − 10 )2 0

if 25 < x < 12.5, 3 , o.w.

on a rectangular channel with the length of 25. The initial and boundary conditions are

h(x, 0 ) = 0.66 − b(x ), hu(x, 0 ) = 1.53, 0  x  25, hu(0, t ) = 1.53,

h(25, t ) = 0.66, 0  t  200.

13

14

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

Fig. 15. Numerical solutions h (x, 200 ) + b(x ) using two methods, moving mesh and fixed mesh with n = 100, and comparison them with the analytical solution.

Fig. 16. Mesh trajectories for problem 6.4, steady flow over a bump, for t ∈ [0, 200].

The analytical solution of h(x, 200 ) + b(x ) and also the numerical solutions with n = 100 and CF L = 0.9 by using the fixed mesh and moving mesh methods have been plotted in Fig. 15. The analytical solution has been computed using the software SWASHES [22]. Both types of numerical solutions are good approximations to this problem. The solutions from both methods are very close to each other and have almost the same accuracy (in some areas, the accuracy of the moving mesh method is still slightly better than the fixed mesh method). In Fig. 16, the mesh trajectories for t ∈ [0.200] with α2 = α4 = 10 and ν = 5 have been presented. We have used this technique for other problems including trans-critical flow with shock [34,70,84], dry bed simulation with discontinuous bathymetry [27], and several other test problems with successful results. However, the results are not included here to make the paper more concise.

7. Conclusion We have presented a moving wave propagation algorithm based on augmented Riemann solver for shallow water equation. The CLAWPACK software package has been developed for hyperbolic conservation laws, and GEOCLAW software package [57], as its subset, has been presented for modelling shallow water equations with a specific inclination toward tsunami problems. Two main components of GEOCLAW are wave propagation algorithm and augmented Riemann solver which have been introduced in Section 3. Since shock and discontinuity exist in most of the shallow water problems, using a moving mesh methods is more efficient than using a fixed mesh method. In this paper, we have improved the augmented Riemann solver using a moving mesh method. Our numerical results indicated that in problems with shock or discontinuity, or discontinuous bed, the numerical solutions based on the moving mesh method and fixed mesh method have a significant difference. Application of a moving mesh method yields highly accurate solutions around shocks and discontinuities, while the number of mesh points are the same in both methods.

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

15

We believe that moving mesh method can improve GEOCLAW software for problems which include shocks. Moreover, this technique can be so efficient for tsunami problems near costs with broken waves which cause discontinuities in the solutions. It can also predict coastal phenomena resulted from tsunami in one or two dimensions.

References [1] F. Alcrudo, F. Benkhaldoun, Exact solutions to the Riemann problem of the shallow water equations with a bottom step, Comput. Fluids 30 (6) (2001) 643–671. [2] F. Alcrudo, E. Gil, The malpasset dam break case study, in: Proceedings of the Fourth Concerted Action on Dambreak Modelling Workshop, Spain, 1999, pp. 95–109. [3] N.A. Alias, Q. Liang, G. Kesserwani, A Godunov-type scheme for modelling 1d channel flow with varying width and topography, Comput. Fluids 46 (1) (2011) 88–93. [4] E. Audusse, F. Bouchut, M.O. Bristeau, R. Klein, B.T. Perthame, A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows, SIAM J. Sci. Comput. 25 (2004) 2050–2065. [5] D. Bale, R.J. Leveque, S. Mitran, J.A. Rossmanith, A wave-propagation method for conservation laws and balance laws with spatially varying flux functions, SIAM J. Sci. Comput. 24 (2002) 955–978. [6] R. Bernetti, V.A. Titarev, E.F. Toro, Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry, J. Comput. Phys. 227 (6) (2008) 3212–3243. [7] C. Berthon, F. Foucher, Efficient well-balanced hydrostatic upwind schemes for shallow-water equations, J. Comput. Phys. 231 (15) (2012) 4993–5015. [8] R.P. Bonet, N. Norberto, M. Storti, S. Idelsohn, An explicit scheme for the numerical solution of the shallow water equation, Mecánica Comput. XIX (1) (20 0 0) 37–42. [9] F. Bouchut, Non-linear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources, Birkhüser Basel (2004). [10] O.V. Bulatov, T.G.E. Elizarova, Regularized shallow water equations for numerical simulation of flows with a moving shoreline, Comput. Math. Math. Phys. 56 (4) (2016) 661–679. [11] J. Burguete, P. Garcia-Navarro, Implicit schemes with large time step for non-linear equations: application to river flow hydraulics, Int. J. Numer. Methods Fluids 46 (6) (2004) 607–636. [12] V. Caleffi, A. Valiani, A. Zanni, Finite volume method for simulating extreme flood events in natural channels, J. Hydraul. Res. 41 (2) (2003) 167–177. [13] W.M. Cao, W.Z. Huang, R.D. Russell, An r-adaptive finite element method based upon moving mesh PDEs, J. Comput. Phys. 149 (1999) 221–244. [14] M. Castro, J.M. Gallardo, J.A. López-GarcÍa, C. Parés, Well-balanced high order extensions of Godunov’s method for semilinear balance laws, SIAM J. Numer. Anal. 46 (2) (2008) 1012–1039. [15] M. Castro, A. Pardo, C. Parés, E. Toro, On some fast well-balanced first order solvers for nonconservative systems, Math. Comput. 79 (271) (2010) 1427–1472. [16] D.M. Causon, D.M. Ingram, C.G. Mingham, G. Yang, R.V. Pearson, Calculation of shallow water flows using a cartesian cut cell approach, Adv. Water Resour. 23 (5) (20 0 0) 545–562. [17] H.D. Ceniceros, T.Y. Hou, An efficient dynamically adaptive mesh for potentially singular solutions, J. Comput. Phys. 172 (2001) 609–639. [18] A. Chertock, S. Cui, A. Kurganov, T. Wu, Well-balanced positivity preserving central-upwind scheme for the shallow water system with friction terms, Int. J. Numer. Methods Fluids 78 (6) (2015) 355–383. [19] S. Clain, S. Diot, R. Loub‘ere, A high-order finite volume method for hyperbolic systems: multi-dimensional optimal order detection (MOOD), J. Comput. Phys. 230 (10) (2011) 4028–4050. [20] D.H. Cuong, M.D. Thanh, A well-balanced van leer-type numerical scheme for shallow water equations with variable topography, Adv. Comput. Math. 43 (5) (2017) 1197–1225. [21] O. Delestre, S. Cordier, F. James, F. Darboux, Simulation of Rain-water Overland-flow, in: Proceedings of the Twelfth International Conference on Hyperbolic Problems, volume 67, American Mathematical Society, 2009, pp. 537–546. [22] O. Delestre, et al., SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies, Int. J. Numer. Methods Fluids 72 (3) (2013) 269–300. [23] D. Doyen, P.H. Gunawan, An explicit staggered finite volume scheme for the shallow water equations, in: Finite Volumes for Complex Applications VII-Methods and Theoretical Aspects, 2014, pp. 227–235. [24] B. Einfeldt, On Godunov-type methods for gas dynamics, SIAM J. Numer. Anal. 25 (1988) 294–318. [25] B. Einfelde, C.D. Munz, P.L. Roe, B. Sjogreen, On Godunov-type methods for near low densities, J. Comput. Phys. 92 (1991) 273–295. [26] J. Felcman, O. Havle, On a numerical flux for the shallow water equations, Appl. Math. Comput. 217 (11) (2011) 5160–5170. [27] J.M. Figueiredo, S. Clain, Second-order finite volume MOOD method for the shallow water with dry/wet interface, in: Proceedings of the SYMCOMP, 2015, pp. 191–205. [28] Q. Gao, S. Zhang, Moving mesh strategies of adaptive methods for solving non-linear partial differential equations, Algorithms 9 (4) (2016) 86. [29] D.L. George, Numerical Approximation of the Nonlinear Shallow Water Equations: a Godunov-Type Scheme, Master’s thesis, University of Washington, 2004. [30] D.L. George, Augmented Riemann solvers for the shallow water equations over variable topography with steady states and inundation, J. Comput. Phys. 227 (6) (2008) 3089–3113. [31] S.K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Mat. Sb. 47 (1959) 271–306. [32] N. Goutal, F. Maurel, Finite volume solver for 1d shallow-water equations applied to an actual river, Int. J. Numer. Methods Fluids 38 (1) (2002) 1–19. [33] J.M. Greenberg, A.Y. Leroux, R. Baraille, A. Noussair, Analysis and approximation of conservation laws with source terms, SIAM J. Numer. Anal. 34 (1997) 1980–2007. [34] P.H. Gunawan, X. Lhébrard, Hydrostatic relaxation scheme for the 1d shallow water-exner equations in bedload transport, Comput. Fluids 121 (2015) 44–50. [35] A. Harten, P.D. Lax, B. van Leer, On upstream differencing and Godunov-tpye schemes for hyperbolic conservation laws, SIAM Rev. 25 (1983) 235–261. [36] R. Herbin, W. Kheriji, J.C. Latch, On some implicit and semi-implicit staggered schemes for the shallow water and euler equations, ESAIM Math. Model. Numer. Anal. 48 (6) (2014) 1807–1857. [37] W. Huang, Practical aspects of formulation and solution of moving mesh partial differential equations, J. Comput. Phys. 171 (2) (2001) 753–775. [38] W. Huang, Y. Ren, R.D. Russell, Moving mesh methods based upon moving mesh partial differential equations, J. Comput. Phys. 113 (1994) 279–290. [39] W. Huang, Y. Ren, R.D. Russell, Moving mesh partial differential equation based on the equidistribution principle, SIAM J. Numer. Anal. 31 (1994) 709–730. [40] W. Huang, R.D. Russell, Analysis of moving mesh partial differential equations with spatial smoothing, SIAM J. Numer. Anal. 34 (3) (1997) 1106–1126. [41] W. Huang, R.D. Russell, A high dimensional moving mesh strategy, Appl. Numer. Math. 26 (1998) 63–76. [42] W. Huang, R.D. Russell, Adaptive mesh movement-the MMPDE approach and its applications, J. Comput. Appl. Math. 128 (1) (2001) 383–398. [43] W. Huang, R.D. Russell, Adaptive Moving Mesh Methods, Springer Science Business Media, 2010. [44] W. Huang, R.D. Russell, Adaptive Moving Mesh Methods, Springer, New York, 2011.

16

M. Bagherpoorfard and A.R. Soheili / Applied Mathematics and Computation 375 (2020) 125087

[45] Y. Huang, N. Zhang, Y. Pei, Well-balanced finite volume scheme for shallow water flooding and drying over arbitrary topography, Eng. Appl. Comput. Fluid Mech. 7 (1) (2013) 40–54. [46] J.M. Hyman, S. Li, L.R. Petzold, An adaptive moving mesh method with static rezoning for partial differential equations, Comput. Math. Appl. 46 (10–11) (2003) 1511–1524. [47] T. Kärnä, B. de Brye, O. Gourgue, J. Lambrechts, R. Comblen, V. Legat, E. Deleersnijder, A fully implicit wetting-drying method for DG-FEM shallow water models, with an application to the Scheldt estuary, Comput. Methods Appl. Mech. Eng. 200 (5–8) (2011) 509–524. [48] G. Kesserwani, Q. Liang, Well-balanced RKDG2 solutions to the shallow water equations over irregular domains with wetting and drying, Comput. Fluids 39 (10) (2010) 2040–2050. [49] G. Kesserwani, Q. Liang, Conservative high-order discontinuous Galerkin method for the shallow water equations with arbitrary topography, Int. J. Numer. Methods Eng. 86 (2011) 47–69. [50] D.H. Kim, Y.S. Cho, Y.K. Yi, Propagation and run-up of nearshore tsunamis with HLLC approximate Riemann solver, Ocean Eng. 34 (8–9) (2007) 1164–1173. [51] A.G. Kulikovskii, N.V. Pogorelov, A.Y. Semenov, Mathematical Aspects of Numerical Solution of Hyperbolic Systems, Chapman and Hall/CRC, 2001. [52] R.J. LeVeque, High resolution finite volume methods on arbitrary grids via wave propagaton, J. Comput. Phys. 78 (1) (1988) 36–63. [53] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser Base, 1992. [54] R.J. LeVeque, Wave propagation algorithms for multi-dimensional hyperbolic systems, J. Comput. Phys. 131 (1997) 327–353. [55] R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave propagation algorithm, J. Comput. Phys. 146 (1998) 346–365. [56] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, 2002. [57] R.J. LeVeque, Clawpack Software. http://www.amath.washington.edu/∼claw. [58] B.W. Levin, M.A. Nosov, Physics of Tsunamis, Springer, Berlin, 2008. [59] Q. Liang, A.G.L. Borthwick, G. Stelling, Simulation of dam and dyke-break hydrodynamics on dynamically adaptive quadtree grids, Int. J. Numer. Methods Fluids 46 (2) (2004) 127–162. [60] Q. Liang, F. Marche, Numerical resolution of well-balanced shallow water equations with complex source terms, Adv. Water Resour. 32 (6) (2009) 873–884. [61] S.J. Liang, J.H. Tang, M.S. Wu, Solution of shallow-water equations using least-squares finite-element method, Acta Mech. Sinica 24 (5) (2008) 523–532. [62] R. Marlow, M.E. Hubbard, P.K. Jimack, Moving mesh methods for solving parabolic partial differential equations, Comput. Fluids. 46 (1) (2011) 353–361. [63] V. Michel-Dansac, C. Berthon, S. Clain, F. Foucher, A well-balanced scheme for the shallow-water equations with topography or manning friction, J. Comput. Phys. 335 (2017) 115–154. [64] S. Mungkasi, Z. Li, S.G. Roberts, Weak local residuals as smoothness indicators for the shallow water equations, Appl. Math. Lett. 30 (2014) 51–55. [65] S. Mungkasi, S.G. Roberts, Numerical entropy production for shallow water flows, ANZIAM J. 52 (2010) 1–17. [66] S. Mungkasi, S.G. Roberts, A finite volume method for shallow water flows on triangular computational grids, Int. Conf. Adv. Comput. Sci. Inf. Sys. (2011) 79–84. [67] S. Mungkasi, S.G. Roberts, Validation of ANUGA hydraulic model using exact solutions to shallow water wave problems, J. Phys. 423 (2013) 012029. [68] S. Noelle, N. Pankratz, G. Puppo, J.R. Natvig, Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows, J. Comput. Phys. 213 (2006) 474–499. [69] C. Parés, M. Castro, On the well-balance property of Roe’s method for nonconservative hyperbolic systems. applications to shallow water systems, ESAIM Math. Model. Numer. Anal. 38 (5) (2004) 821–852. [70] T. Pongsanguansin, M. Maleewong, K. Mekchay, Consistent weighted average flux of well-balanced TVD-RK discontinuous Galerkin method for shallow water flows, Model. Simul. Eng. 2015 (2015) 591282. [71] S. Popinet, Quadtree-adaptative tsunami modelling, Ocean Dyn. 61 (9) (2011) 1261–1285. [72] P.L. Roe, Approximate riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys. 43 (1981) 357–372. [73] J.A. Rossmanith, A Wave Propagation Method With Constrained Transport for Ideal and Shallow Water Magnetohydrodynamics, PhD thesis, University of Washington, 2002. [74] A.R. Soheili, J.M. Stockie, A moving mesh method with variable mesh relaxation time, Appl. Numer. Math. 58 (3) (2008) 249–263. [75] J.M. Stockie, J.A. Mackenzie, R.D. Russell, A moving mesh method for one dimensional hyperbolic conservation laws, SIAM J. Sci. Comput. 22 (2001) 1791–1813. [76] H. Tang, T. Tang, Adaptive mesh methods for one and two-dimensional hyperbolic conservation laws, SIAM J. Numer. Anal. 41 (2003) 487–515. [77] H. Tang, T. Tang, K.A. Xu, A gas-kinetic scheme for shallow-water equations with source terms, Z. für Angew. Math. 290 Phys. ZAMP 55 (3) (2004) 365–382. [78] E.F. Toro, Shock Capturing Methods for Free Surface Shallow Flows, John Wiley and Sons, 2001. [79] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer-Verlag, 2009. [80] E.T. Toro, P. Garcia-Navarro, Godunov-type methods for free-surface shallow flows: a review, J. Hydraul. Res. 45 (6) (2007) 736–751. [81] A. Valiani, V. Caleffi, A. Zanni, Case study: malpasset dam-break simulation using a two-dimensional finite volume method, J. Hydraul. Eng. 128 (5) (2002) 460–472. [82] S. Vater, N. Beisiegel, J. Behrens, A limiter-based well-balanced discontinuous Galerkin method for shallow-water flows with wetting and drying: one dimensional case, Adv. Water Res. 85 (2015) 1–13. [83] J.J. Wijetunge, Numerical simulation of the 2004 indian ocean tsunami: case study of effect of sand dunes on the spatial distribution of inundation in hambantota, sri lanka, J. Appl. Fluid Mech. 3 (2010) 125–135. [84] K. Wu, H. Tang, A newton multigrid method for steady-state shallow water equations with topography and dry areas, Appl. Math. Mech. 37 (11) (2016) 1441–1466. [85] Y. Xing, C.W. Shu, High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms, J. Comput. Phys. 214 (2) (2006) 567–598. [86] F. Zhou, G. Chen, S. Noelle, H. Guo, A well balanced stable generalized Riemann problem scheme for shallow water equations using adaptive moving unstructured triangular meshes, Int. J. Numer. Methods Fluids 73 (3) (2013) 266–283. [87] Z.G. Zhou, D.M. Causon, C.G. Mingham, D.M. Ingram, The surface gradient method for the treatment of source terms in the shallow-water equations, J. Comput. Phys. 168 (1) (2001) 1–25.