Abstract

High order methods are becoming increasingly popular in shallow water flow modeling motivated by their high computational efficiency (i.e. the ratio between accuracy and computational cost). In particular, Discontinuous Galerkin (DG) schemes are very well suited for the resolution of the Shallow Water Equations and related models, being a competitive alternative to the traditional finite volume schemes. In this work, a novel framework for the construction of DG schemes using augmented Riemann solvers is proposed. Such solvers incorporate the source term at cell interfaces in the definition of the Riemann problem, allowing definition of two different inner states in the so-called star region. The benefits of this family of solvers lie in the exact preservation of the Rankine–Hugoniot condition at cell interfaces at the discrete level, ensuring the preservation of equilibrium solutions (i.e. the well-balanced property) without requiring extra corrections of the numerical fluxes. The semi-discrete DG operator becomes nil automatically under equilibrium conditions, provided the use of suitable quadrature rules. The proposed scheme is applied to the Burgers' equation with geometric source term and to the SWE. The numerical results evidence that the proposed scheme achieves the prescribed convergence rates and preserves the equilibrium states of relevance with machine precision.

INTRODUCTION

Many engineering and environmental problems involve the study of steady and transient free surface water flows where the vertical scale is much smaller than the horizontal dimensions. Such phenomena are usually described by the Shallow Water Equations (SWE) (Cunge et al. 1980), a depth-averaged model that considers a hydrostatic pressure distribution in the vertical direction. In realistic scenarios, the flow is often dominated by the so-called source terms (e.g. bed topography, friction, Coriolis forces, etc.), which may appear as a source of mass and/or momentum on the right-hand side of the equations. The methods herein presented focus on the numerical resolution of the SWE with bed slope source term, which is a non-linear, non-homogeneous hyperbolic system of conservation laws.

Suitable discretization techniques for the source term are thus required to provide physically feasible solutions. For a numerical scheme, the preservation of equilibrium solutions is a feature of paramount importance, as many realistic problems of physical interest can be regarded as perturbations over an equilibrium state (Chertock et al. 2017). For the SWE system, the most relevant equilibrium condition is known as the lake-at-rest equilibrium (Bermudez & Vázquez-Cendón 1994). This work aims at the design of a novel scheme for the resolution of the SWE and related models in the framework of the Discontinuous Galerkin (DG) method. To this end, the equilibrium condition is generalized for an arbitrary number of spatial dimensions and equations, so that the methods can be presented in a more general form. As a result, the proposed scheme can also be applied to other problems such as the Burgers' equation with source term.

The DG method is a common choice for the numerical resolution of hyperbolic partial differential equations. The origin of this method is usually attributed to Reed and Hill in a paper published in 1973 on the numerical approximation of the neutron transport equation (Reed & Hill 1973). Later on, Cockburn and Shu presented the fundamentals of the scheme in a series of papers published in the 1980s and the 1990s (Cockburn & Shu 1989; Cockburn & Shu 1990). After some years, the DG method became more popular and was successfully applied to the resolution of evolution problems in continuum mechanics, among others. The DG scheme can be seen as a compromise between finite elements (compactness of the stencils, easy -refinement and parallel scalability) and finite volumes (intrinsic conservativeness, the presence of discontinuities between cells and the use of numerical fluxes). Some of the advantages of the DG method, namely its high accuracy, inherent conservativeness, compactness, high parallel efficiency, flexibility for -adaptivity and possibility of using arbitrary geometry and meshes, make it very well suited for the numerical resolution of the SWE (Xing et al. 2010; Caviedes-Voullième et al. 2020), which is the focus of this work.

The most relevant characteristic that finite volumes and DG methods share may presumably be the use of numerical fluxes at cell interfaces, which allow to exchange information between cells. A suitable calculation of such fluxes is essential for an accurate approximation of the numerical solution ensuring conservation and equilibrium properties. In a DG and/or finite volume framework, numerical fluxes are usually computed by means of the resolution of the Riemann problem (RP) at cell interfaces, which is defined by the original set of equations together with a piecewise constant initial condition given by the DG reconstructed data at cell interfaces. Under equilibrium, the Rankine–Hugoniot (RH) condition must hold across discontinuities in the solution (e.g. steady shock waves or contact waves). In the presence of geometric source terms, the piecewise DG representation of the geometric variable (i.e. polynomial reconstruction with discontinuities at cell interfaces) involves the presence of contact waves at cell interfaces. In this case, the RH condition states that the difference of fluxes across the interface must balance the contribution of the source term, introduced by such a contact wave (Murillo & García-Navarro 2010; Rosatti & Begnudelli 2010; LeVeque 2011).

Most traditional approaches found in the literature for the exact preservation of the RH condition (and ultimately of the equilibrium solution) in the presence of geometric source terms are based on extra modifications of the numerical fluxes (e.g. using the hydrostatic reconstruction, enforcing continuity of the bed elevation at cell interfaces, etc.), which proved to be very robust and useful (Audusse et al. 2004; Xing & Shu 2006). However, the solvers therein used do not explicitly account for the full eigenstructure of the system, i.e. they do not account for the contact wave due to the source term in the numerical solution of the RP. To address this issue, augmented Riemann solvers were introduced in an attempt to represent the effect of the source term in the solution of the RP and automatically fulfill RH conditions. When dealing with geometric source terms, it was shown that the consideration of the source term in the definition of the RP is more convenient in order to ensure certain properties of the numerical solution (George 2008; Murillo & García-Navarro 2010). Augmented solvers have been successfully applied to high order schemes (e.g. WENO-ADER and WENO-RK type schemes) (Navas-Montilla & Murillo 2016; Murillo et al. 2019; Navas-Montilla et al. 2019), allowing well-balanced and energy-balanced numerical solutions that hold the RH condition at cell interfaces with machine precision under equilibrium without extra corrections. Moreover, the augmented solver approach has also been useful to circumvent numerical shockwave anomalies in the finite volume framework (Navas-Montilla & Murillo 2019) and to model realistic scenarios with complex physics (Franzini & Soares-Frazao 2018; Gordillo et al. 2019; Martínez-Aranda et al. 2019).

As mentioned before, the preservation of certain equilibrium states at the discrete level is a fundamental requirement for a numerical scheme. When considering the SWE with bottom topography, the relevant equilibrium state is the so-called quiescent equilibrium, also known as ‘lake-at-rest’, where the water surface elevation is constant. Numerical schemes satisfying such an equilibrium state are called well-balanced methods (Bermudez & Vázquez-Cendón 1994; Greenberg & Leroux 1996). There is a large variety of well-balanced methods based on Riemann solvers that ensure the preservation of the still water steady state (García-Navarro & Vázquez-Cendón 2000; Alcrudo & Benkhaldoun 2001; Burguete & García-Navarro 2001, 2004; Sanders et al. 2003; Chinnayya et al. 2004; Yulong & Shu 2005; LeFloch & Thanh 2007, 2011; Bernetti et al. 2008; Catella et al. 2008; Kesserwani et al. 2008; Lee & Wright 2010; Rosatti & Begnudelli 2010; Haleem et al. 2015). In 2015, Caleffi et al. present an exhaustive study on the treatment of the bed step source term for well-balanced DG schemes, both in the framework of the classical finite volume approach and path-conservative methods. They examine the technique by Kesserwani & Liang (2011) that modifies the discrete bed profile data imposing continuity at cell interfaces, the widespread hydrostatic reconstruction method (Audusse et al. 2004) and a path-conservative scheme based on the Dumbser–Osher–Toro (DOT) Riemann solver (Dumbser & Toro 2011). The latter is based on rewriting the SWE as a non-conservative system which includes the trivial equation for the temporal variation of the bed elevation (i.e. representing the conservative product). This approach has been successfully applied to the SWE and enhanced models (Escalante et al. 2019) and represents an alternative to augmented solvers, explicitly including the contact wave of the source term in the eigenstructure of the Jacobian. It must be noted that for a hyperbolic system which includes non-conservative products, the choice of the correct path, necessary to define the correct weak solution, is generally non-trivial when complex equilibrium states are sought (Caleffi et al. 2015).

The application of augmented Riemann solvers for the computation of the numerical fluxes at cell interfaces represents an alternative to the aforementioned methods and, to the knowledge of the authors, it has never been used in a DG framework. The work herein described focuses on the construction of a DG scheme using augmented Riemann solvers for the resolution of hyperbolic conservation laws with geometric source terms, motivated by the benefits of such family solvers. It is herein evidenced that DG schemes are automatically well balanced when using augmented solvers (i.e. when the numerical fluxes satisfy the RH condition at cell interfaces), provided (i) the use of an exact quadrature rule for the surface and volume integrals and (ii) the reconstructed data satisfy the discrete equilibrium. The time stepping is carried out using an explicit strong stability preserving (SSP) Runge–Kutta (RK) integrator. The proposed methods are applied to the one-dimensional (1D) Burgers' equation with source term and to the SWE with bed topography. For the latter, the HLLS solver is applied, which approaches the eigenstructure of the system by two non-linear waves associated with the convective flux and an extra contact wave associated with the bed discontinuity at cell interfaces. The performance of the HLLS solver in the context of DG-RK schemes is assessed when considering an additional transport equation in the SWE, thus an additional contact wave.

The analytical derivation of the proposed scheme, as well as the proof of the well-balanced property, is carried out in 1D and extended in the Supplementary material, Appendix B, for a general hyperbolic system of equations without any restriction on the spatial dimension. Furthermore, the concept of isotropic equilibrium is introduced as a generalization of the hydrostatic equilibrium. It appears when the matrix of physical fluxes, evaluated at an equilibrium state, is spherical (i.e. proportional to the identity matrix).

The performance of the proposed schemes is herein assessed using different test cases, including steady and transient solutions, that allow to evaluate the well-balanced property and the empirical rates of convergence. The numerical results evidence that the proposed schemes preserve the quiescent equilibrium with machine precision and achieve the expected convergence rates up to third order of accuracy, limited here by the use of the SSPRK3 method.

The paper is organized as follows. In the ‘Mathematical model’ section, the mathematical formulation of the problems of interest is presented. In the ‘General formulation of the well-balanced DG-RK method’ section, DG schemes are introduced and the proposed DG-SSPRK3 augmented scheme is presented. The ‘Application to the Burgers’ equation’ and ‘Application to the SWE’ sections include the numerical results of the application of the proposed scheme to Burgers' and the SWE, respectively, and finally, in the last section, the concluding remarks are presented.

MATHEMATICAL MODEL

The methods herein considered focus on the resolution of non-homogeneous hyperbolic systems of conservation laws, which are generally expressed as follows: 
formula
(1)
where is the vector of conserved quantities that takes values on , with is the number of equations. is the matrix of physical fluxes. The source term is defined as a geometric source term: 
formula
(2)
with a geometric function that depends upon the position (e.g. z is the bed elevation in the SWE). In this work, is defined as follows: 
formula
(3)
where is the identity matrix, u is a component of related to the equilibrium as it will be shown later (e.g. the water depth, h, in the SWE) and is a non-linear function of u. We will generally consider that , where a is constant.
In this work, we are interested in isotropic equilibrium (i.e. steady) states where only one variable of the system, u, is related to the equilibrium. Let be a steady solution of (1). The equilibrium state will be called isotropic if the physical flux matrix at equilibrium is of the form: 
formula
(4)
The left-hand side of (1) under equilibrium reads 
formula
(5)
which can be combined with (2) to yield 
formula
(6)
with the equilibrium value of one of the components of vector that participates in the equilibrium condition. Note that the hydrostatic equilibrium is a particular case of the isotropic equilibrium.
We can define the so-called equilibrium variable as 
formula
(7)
which is constant under equilibrium. Note that u is just the variable of the problem (a component of the vector ) which participates in the equilibrium condition.
In what follows, the 1D version of (1) will be considered as follows: 
formula
(8)
noticing that in this case, the flux matrix is reduced to vector , which stands for the physical flux in the x-direction according to the standard notation and will include a single derivative of the geometric variable, .

GENERAL FORMULATION OF THE WELL-BALANCED DG-RK METHOD

Let us consider the 1D problem in Equation (8) (see Supplementary material, Appendix B, for the multi-dimensional extension). The spatial domain, , is discretized in N volume cells, denoted by , such that . The DG method can be obtained by applying a traditional Galerkin projection method to each element. Inside each element, the solution is approximated by a linear combination of basis functions as follows: 
formula
(9)
where is the l-th degree of freedom and is the number of degrees of freedom on each element. It is worth noticing that the basis functions only depend upon space and the degrees of freedom are a function of time only.
Projecting the governing equation onto each element of the basis set, the weak form of the problem can be written as follows: 
formula
(10)
for , and integrating by parts the second term of (10) using Green's first identity: 
formula
(11)
As the solution at the interfaces and is bi-valued, the concept of numerical flux has to be used and the previous identity is rewritten as follows: 
formula
(12)
with and are the cell entering and leaving numerical fluxes, respectively. Note that we will consider different left and right numerical fluxes computed from the solution of the RP at the interface, denoted by and , respectively. Using the equality in Equation (12), we can rewrite Equation (10) as follows: 
formula
(13)
for .
Substitution of in Equation (13) yields 
formula
(14)
When choosing an orthogonal set of basis functions , the following property is satisfied 
formula
(15)
for , where is the Kronecker delta and is the value of the integral over . Using the result in (15), (14) becomes a decoupled system of equations for the unknown degrees of freedom that reads 
formula
(16)
for . Equation (16) can be rewritten as follows: 
formula
(17)
where 
formula
(18)
To compute the integrals in (18), Gaussian quadrature is employed. For each cell , the volume integrals in (18) are computed as follows: 
formula
(19)
and 
formula
(20)
where is the number of the Gaussian quadrature points; , and are the evaluation of the flux, source term and the basis function at the Gaussian quadrature points q; is the Gaussian weights and is the determinant of the Jacobian matrix associated with the change of coordinates from the reference element, , to the real element, . In the 1D cases presented in this work, we consider , provided the reference interval is .
Equation (17) can be integrated over the time step using a TVD RK method. The third order SSP RK method will be herein used: 
formula
(21)
which allow to update the degrees of freedom at , denoted by .

Augmented Riemann solvers

When dealing with hyperbolic conservation laws with geometric source terms (e.g. Equation (1)), non-homogeneous RPs naturally arise at quadrature points across cell interfaces motivated by the discontinuous representation of the solution provided by the DG approach. In such cases, the numerical fluxes at cell interfaces, and , must be computed using augmented Riemann solvers.

Non-homogeneous RPs with geometric source terms are 1D initial value problems generally expressed as follows: 
formula
(22)
where and are the left and right reconstruction of the variables at cell interfaces, provided by the DG polynomials in cells i and , respectively; only includes the -component and is a singular source term defined at the interface such that , with the property .
Under equilibrium (i.e. steady state), the RH condition 
formula
(23)
holds for Equation (22). For a numerical scheme based on augmented solvers, preserving the RH condition at cell interfaces is a necessary condition for well-balancing.
Augmented Riemann solvers include the contribution of the source term as an extra wave of zero velocity at the interface (i.e. ), which originates a jump of the fluxes and conserved variables across the interface. The approximate solution, hereafter referred to as , is therefore composed by at least (depending on the number of waves) two different internal states in the so-called star region, separated by a contact discontinuity at produced by the source term. In the vicinity of , left and right states of the approximate solution will be hereafter denoted by and , respectively, expressed as follows: 
formula
(24)
Analogously, an approximate flux function can be also constructed from flux values and with a similar structure than . In this case, also intercell values for the fluxes can be defined at both sides of the -axis as follows: 
formula
(25)
Augmented Riemann solvers are designed to satisfy the following relation: 
formula
(26)
which can also be expressed as follows: 
formula
(27)
noticing that, under steady state conditions (i.e. and ), Equation (23) is recovered. Note that the numerical fluxes will be generally expressed using the subscript notation as and .

Resolution of the augmented RP for a scalar equation

When addressing the resolution of the augmented RP for a scalar equation, the approach proposed in Murillo & García-Navarro (2010) can be applied. The scalar version of (22) reads 
formula
(28)
For this problem, the left and right numerical fluxes in (25) can be computed as follows (Murillo & García-Navarro 2010): 
formula
(29)
with and 
formula
(30)
where is the approximate wave celerity at the interface.

Resolution of the augmented RP for a system of conservation laws: the HLLS solver

The HLLS solver can be used for the approximation of the solution of (22). It is a non-linear solver that approaches the wave structure of the system by a 2-wave structure. When solving systems with waves (e.g. 1D SWE), the HLLS solver accounts for the full eigenstructure of the system. The numerical fluxes are computed as follows: 
formula
(31)
 
formula
(32)
where the HLLS fluxes are given by the following equation: 
formula
(33)
 
formula
(34)
where and are the left- and right-hand limits to the cell edge of the conserved quantities () and their -th time derivatives and is a matrix related to the source term. Further details can be found in Navas-Montilla & Murillo (2016).

Well-balancing the scheme

Well-balanced schemes are able to preserve the isotropic equilibrium in Equation (6) with machine accuracy. This means that the spatial operator in Equation (18) must vanish in such an equilibrium condition. The particular requirements a DG method must fulfill in order to be well-balanced are stated in Xing & Shu (2006) and Caleffi & Valiani (2012) as follows:

  • Under equilibrium, the reconstructed data (using the DG projection) must satisfy the equilibrium condition: this is achieved by applying the DG reconstruction over the equilibrium variable and the geometric variable, and then computing the conserved variable from the difference between these.

  • Under equilibrium, the numerical flux at the interface must be equal to the physical flux: this is automatically achieved by augmented solvers.

  • Under equilibrium, the quadrature rule for the surface and volume integrals must be exact: this is achieved by choosing the adequate number of quadrature points, provided the polynomial order.

Note that other different approaches for the construction of well-balanced DG schemes that do not require some of the conditions above can be found in the literature, in particular Wintermeyer et al. (2017) proposed a well-balanced DG method for the SWE with inexact integration. These approaches have recently been extended to solve complex multiphase flow problems (Manzanero et al. 2020).

In this section, the aforementioned requirements are proved to be fulfilled. The well-balanced equilibrium condition is first derived using the analytical form of the spatial operator in Equation (18). Then, the conditions for the numerical approximation of the integrals in order to preserve the discrete equilibrium with machine accuracy are stated.

Under steady state (e.g. the isotropic equilibrium in Equation (6)), the -th cell-leaving numerical fluxes in Equation (18) are equal to the physical fluxes in the cell , according to Equation (27). This allows Equation (18) to be rewritten as follows: 
formula
(35)
where is a steady-state equilibrium solution.
The integrals in (35) must be computed numerically using the Gaussian quadrature following Equations (19) and (20). For the numerical scheme to satisfy the well-balanced property at the discrete level, the quadrature rules must be of sufficient order of accuracy to exactly integrate the product of polynomial functions. In that case, Equation (35) would yield to 
formula
(36)
which guarantees the preservation of the discrete equilibrium with machine precision. Table 1 summarizes the polynomial degree, truncation order and maximum integrand degree in the spatial operator, up to fourth order of accuracy.
Table 1

Polynomial degree, truncation order and maximum integrand degree (MID) in the spatial operator for 1D and 2D/3D problems

p (degree)Truncation errorMID (1D)MID (2D,3D)
 
 
 
p (degree)Truncation errorMID (1D)MID (2D,3D)
 
 
 

The choice of a suitable orthogonal basis: Legendre polynomials

Legendre polynomials are herein used to construct a polynomial basis for the DG formulation. They are suitable as they are orthogonal with respect to the norm, which satisfies Equation (15). A recursive formula for Legendre polynomials, hereafter denoted by , is given by the following equation: 
formula
(37)
Analogously, a recursive formula for the derivatives of the polynomials, , can be obtained as follows: 
formula
(38)

APPLICATION TO THE BURGERS’ EQUATION

The application of the novel scheme to the scalar inviscid Burgers' equation with the geometric source term is considered in this section. This equation represents a scalar version of the problem in Equations (1) and (2) and reads 
formula
(39)
where is the problem variable and is a geometric variable. The equilibrium condition for Equation (39) is obtained when vanishes, yielding 
formula
(40)
which allows definition of an equilibrium variable as . Steady solutions of (39) will satisfy that e is constant.
For the Burgers' equation, the wave celerity is given by and the approximate wave celerity yields . The approximation of the integral of the source term at cell interfaces is selected as follows: 
formula
(41)
with . This allows the numerical fluxes to be written as follows: 
formula
(42)
where it is straightforward to notice that the fluctuations cancel out under the equilibrium condition of Equation (40). The numerical fluxes thus hold the condition in Equation (23), allowing the well-balancing of the DG spatial operator. It must be borne in mind that the DG projection must be directly applied to the equilibrium variable, e, and the geometric variable, z. Then, the projection of the advected variable (see Equation (9)) is computed as follows: 
formula
(43)

Case 1. Preservation of an initial equilibrium

The preservation of an initial equilibrium condition at the discrete level is tested here. The following initial condition 
formula
(44)
is considered in the domain and is depicted in Figure 1. The solution is computed at using the DG2, DG3 and DG4 schemes in combination with an SSPRK3 integrator, setting and and . The initial condition is in equilibrium and must thus be preserved by the well-balanced version of the schemes. To assess the performance of the schemes in the preservation of the initial equilibrium, the infinity error norm for the equilibrium variable e, denoted as , is used. It is defined as an integral error norm as follows: 
formula
(45)
with q looping over the quadrature points.
Figure 1

Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 1

Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Results for are presented for the well-balanced and non-well-balanced version of the scheme shown in Table 2. Such a non-well-balanced version of the scheme is constructed by means of a homogeneous flux, that is, without including the source term in the definition of the RP, thus preventing the exact balance between fluxes and sources. It is observed that the well-balanced version of the schemes successfully ensures the well-balanced property with machine accuracy. Differently, the non-well-balanced schemes exhibit substantial errors, which are reduced as the order is increased.

Table 2

error at for the well-balanced and non-well-balanced DG2-SSPRK3, DG3-SSPRK3 and DG4-SSPRK3 schemes in two grids with and

Non-well-balanced
Well-balanced
MeshDG2DG3DG4DG2DG3DG4
 9.22 × 10−03 2.78 × 10−03 7.89 × 10−06 8.88 × 10−16 4.44 × 10−16 4.88 × 10−15 
 2.28 × 10−03 6.96 × 10−04 4.90 × 10−07 1.77 × 10−15 4.44 × 10−16 5.77 × 10−15 
Non-well-balanced
Well-balanced
MeshDG2DG3DG4DG2DG3DG4
 9.22 × 10−03 2.78 × 10−03 7.89 × 10−06 8.88 × 10−16 4.44 × 10−16 4.88 × 10−15 
 2.28 × 10−03 6.96 × 10−04 4.90 × 10−07 1.77 × 10−15 4.44 × 10−16 5.77 × 10−15 

Case 2. Perturbation of the equilibrium state

Most events of physical interest are often considered perturbations of an equilibrium state. In this case, the propagation of a small perturbation over the equilibrium state presented in Equation (44) is computed. The initial condition for e, depicted in Figure 2, is now given by 
formula
(46)
Figure 2

Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 2

Initial condition for u (blue) and e (green) as well as the geometric variable, z (black). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

is given by 
formula
(47)
and . The computational domain is given by . The solution is computed at () using the DG2, DG3 and DG4 schemes in combination with an SSPRK3 integrator, in two different meshes with and .

Figure 3 shows the exact solution and numerical results for e computed by the different schemes mentioned above. It is observed that the well-balancedness of the scheme is crucial for an adequate resolution of the propagation of the perturbation. The results provided by the non-well-balanced schemes show the presence of spurious waves (which are not noticeable for the fourth order scheme) that are only reduced with mesh refinement and by increasing the order of accuracy. Note that the exact solution displayed in Figure 3 is computed using the procedure in Algorithm 1 (Supplementary material, Appendix A), which is based on a fine evolution of the solution using the conservation property along the characteristic curves.

Figure 3

Numerical solution for e computed by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (purple), DG3-SSPRK3 (yellow) and DG4-SSPRK3 schemes (green) at in a (left) and grid (right). The initial condition is depicted using the black dashed line and the exact solution using the black solid line. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 3

Numerical solution for e computed by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (purple), DG3-SSPRK3 (yellow) and DG4-SSPRK3 schemes (green) at in a (left) and grid (right). The initial condition is depicted using the black dashed line and the exact solution using the black solid line. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 4

Evolution of the computed e provided by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (left), DG3-SSPRK3 schemes (right) at in a grid.

Figure 4

Evolution of the computed e provided by the non-well-balanced (top) and well-balanced version (bottom) of the DG2-SSPRK3 (left), DG3-SSPRK3 schemes (right) at in a grid.

A representation of the evolution of e in time is depicted in Figure 4. The numerical results evidence that two different spurious perturbations appear in the non-well-balanced solutions: (i) a traveling perturbation that moves at velocity u and (ii) a steady perturbation over the geometrical hump around . Such spurious waves do not appear in the well-balanced solutions.

Case 3. Convergence rate test

The order of the truncation error of the solution can be estimated using the double-mesh principle (Murillo et al. 2019). A mesh refinement ratio of 2 is considered, leading to the following set of mesh's number of cells: . For two arbitrary meshes and , with the number of cells, the norm of the difference of the computed u in such meshes is calculated as follows: 
formula
(48)
with , 
formula
(49)
where stands for the grid in which the solution is computed and l is the cell number. Additionally, the error norm can also be defined as follows: 
formula
(50)
The test case is configured as follows. The initial condition is given by 
formula
(51)
which is defined in the domain . The solution is computed at () using the DG2, DG3 schemes in combination with an SSPRK3 integrator, in the meshes indicated above. Table 3 shows the relative numerical errors, computed using the and error norms, and empirical convergence rates for u. The numerical results evidence that the prescribed order of accuracy is achieved.
Table 3

Convergence rates for the and error norms, computed by the well-balanced DG2-SSPRK3 and DG3-SSPRK3 schemes

WB DG2-SSPRK3
WB DG3-SSPRK3
Meshes ()
 2.12 × 10−05 – 1.27 × 10−04 – 2.14 × 10−07 – 1.08 × 10−06 – 
 2.95 × 10−06 2.85 1.85 × 10−05 2.78 8.80 × 10−09 4.61 5.71 × 10−08 4.25 
 4.18 × 10−07 2.82 2.53 × 10−06 2.86 6.01 × 10−10 3.87 4.05 × 10−09 3.82 
 6.48 × 10−08 2.69 3.78 × 10−07 2.74 6.06 × 10−11 3.31 4.09 × 10−10 3.31 
 1.15 × 10−08 2.49 7.03 × 10−08 2.42 7.16 × 10−12 3.08 4.80 × 10−11 3.09 
 2.375 × 10−09 2.28 1.56 × 10−08 2.16 8.84 × 10−13 3.02 5.80 × 10−12 3.05 
 5.39375 × 10−10 2.13 3.73 × 10−09 2.06 1.10 × 10−13 3.01 7.12 × 10−13 3.03 
WB DG2-SSPRK3
WB DG3-SSPRK3
Meshes ()
 2.12 × 10−05 – 1.27 × 10−04 – 2.14 × 10−07 – 1.08 × 10−06 – 
 2.95 × 10−06 2.85 1.85 × 10−05 2.78 8.80 × 10−09 4.61 5.71 × 10−08 4.25 
 4.18 × 10−07 2.82 2.53 × 10−06 2.86 6.01 × 10−10 3.87 4.05 × 10−09 3.82 
 6.48 × 10−08 2.69 3.78 × 10−07 2.74 6.06 × 10−11 3.31 4.09 × 10−10 3.31 
 1.15 × 10−08 2.49 7.03 × 10−08 2.42 7.16 × 10−12 3.08 4.80 × 10−11 3.09 
 2.375 × 10−09 2.28 1.56 × 10−08 2.16 8.84 × 10−13 3.02 5.80 × 10−12 3.05 
 5.39375 × 10−10 2.13 3.73 × 10−09 2.06 1.10 × 10−13 3.01 7.12 × 10−13 3.03 

APPLICATION TO THE SWE

The 1D, -split SWE represent the projection of the 2D SWE system onto the -axis, as a result of the rotational invariance property of the SWE (Murillo & Navas-Montilla 2016). They are usually solved at cell interfaces when addressing the computation of 2D shallow water flows. The 1D, -split SWE over the irregular bed are expressed as follows: 
formula
(52)
where is the vector of conserved variables, namely the water depth, h, the unit discharge in x, , and the unit discharge in the transverse direction, ; is the physical flux; u is the depth-averaged velocity and m/s2 is the acceleration of gravity. The unit discharge, , is also denoted by q. The source term accounts for the variations in bed geometry: 
formula
(53)
where z is the geometric variable, representing the bed elevation. Note that the third equation in the system in (52) involves a contact wave in the eigenstructure of the Jacobian matrix of and is associated with the passive transport of with celerity u.
The application of the novel solver to the SWE with the bed slope (geometric) source term in Equation (52) is considered in this section. The HLLS solver is used to compute the numerical fluxes at cell interfaces. The following approximation of the integral of the source term at cell interfaces is chosen: 
formula
(54)
which was proven to hold the equality in Equation (23) under quiescent equilibrium conditions.
As in the Burgers' equation, the DG reconstruction must be applied to the equilibrium variable, , the geometric variable, z, and discharge, . The degrees of freedom of the water depth will be computed as follows: 
formula
(55)

Recall that the conservation of the specific discharge in the momentum equation is affected by the use of different left and right numerical fluxes at cell interfaces, allowing to ensure the exact equilibrium between momentum fluxes and the force exerted by the change of bathymetry. On the other hand, the conservation of mass, i.e. variable h, is not affected as the source term only acts in the momentum equation.

Case 1. Preservation of an initial equilibrium

The well-balanced property of the proposed scheme is herein verified. A quiescent initial condition, , is set, where z is defined as follows: 
formula
(56)
is considered in the domain and is depicted in Figure 5. The solution is computed at using the DG2 and DG3 schemes in combination with an SSPRK3 integrator, setting and and . For the latter mesh, more than 2,000 time steps are taken. The initial condition is in equilibrium and must thus be preserved by the well-balanced version of the schemes. To assess the performance of the schemes, the infinity error norm for the equilibrium variable, , is used.
Figure 5

Initial (green) and bed elevation z (black).

Figure 5

Initial (green) and bed elevation z (black).

Figure 6

Numerical solution for (left) and (right) computed by the well-balanced first order Godunov's scheme (green), DG2-SSPRK3 (purple) and DG3-SSPRK3 (orange) at in a grid. The reference solution, computed by a third order WENO-ADER scheme (Navas-Montilla & Murillo 2016) in a very fine grid (), is depicted in black. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 6

Numerical solution for (left) and (right) computed by the well-balanced first order Godunov's scheme (green), DG2-SSPRK3 (purple) and DG3-SSPRK3 (orange) at in a grid. The reference solution, computed by a third order WENO-ADER scheme (Navas-Montilla & Murillo 2016) in a very fine grid (), is depicted in black. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Results for are presented for the well-balanced version of the schemes shown in Table 4. It is observed that the proposed schemes successfully ensure the well-balanced property with machine accuracy.

Table 4

error at for the well-balanced DG2-SSPRK3 and DG3-SSPRK3 schemes in two grids with and

MeshDG2-SSPRK3DG3-SSPRK3
 2.22 × 10−16 1.00 × 10−16 
 8.88 × 10−16 3.11 × 10−15 
MeshDG2-SSPRK3DG3-SSPRK3
 2.22 × 10−16 1.00 × 10−16 
 8.88 × 10−16 3.11 × 10−15 

Case 2. Convergence rate test

In this test case, the convergence rates of the scheme are evaluated for the variables h and , involved in the non-linear dynamics of Equation (52). To empirically assess the order of convergence of the proposed scheme, the following initial condition is chosen: 
formula
(57)
which is defined in the domain . The topography is given by the smooth function . The solution is computed at () using the HLLS, DG2 and DG3 schemes in combination with an SSPRK3 integrator, in eight different meshes with the number of cells 20, 40, 80, 160, 320, 640, 1,280, 2,560. The numerical solution using 20 cells is depicted in Figure 6. Tables 5 and 6 show the relative numerical errors, computed using the and error norms, and empirical convergence rates for h and , respectively. The numerical results evidence that the prescribed order of accuracy is achieved.
Table 5

Convergence rates for the water depth and discharge, h and , associated with the and error norms, computed by the well-balanced DG3-SSPRK3 schemes



Meshes ()
 1.47 × 10−05 – 7.17 × 10−05 – 6.30 × 10−05 – 4.22 × 10−04 – 
 7.98 × 10−07 4.20 5.34 × 10−06 3.74 3.54 × 10−06 4.15 3.18 × 10−05 3.73 
 4.67 × 10−08 4.09 3.47 × 10−07 3.94 2.09 × 10−07 4.08 1.66 × 10−06 4.25 
 4.17 × 10−09 3.48 2.97 × 10−08 3.54 1.86 × 10−08 3.48 1.40 × 10−07 3.571 
 4.72 × 10−10 3.14 3.30 × 10−09 3.17 2.11 × 10−09 3.14 1.55 × 10−08 3.17 
 5.74 × 10−11 3.03 4.00 × 10−10 3.04 2.57 × 10−10 3.03 1.88 × 10−09 3.04 
 7.12 × 10−12 3.01 4.93 × 10−11 3.02 3.19 × 10−11 3.01 2.33 × 10−10 3.01 


Meshes ()
 1.47 × 10−05 – 7.17 × 10−05 – 6.30 × 10−05 – 4.22 × 10−04 – 
 7.98 × 10−07 4.20 5.34 × 10−06 3.74 3.54 × 10−06 4.15 3.18 × 10−05 3.73 
 4.67 × 10−08 4.09 3.47 × 10−07 3.94 2.09 × 10−07 4.08 1.66 × 10−06 4.25 
 4.17 × 10−09 3.48 2.97 × 10−08 3.54 1.86 × 10−08 3.48 1.40 × 10−07 3.571 
 4.72 × 10−10 3.14 3.30 × 10−09 3.17 2.11 × 10−09 3.14 1.55 × 10−08 3.17 
 5.74 × 10−11 3.03 4.00 × 10−10 3.04 2.57 × 10−10 3.03 1.88 × 10−09 3.04 
 7.12 × 10−12 3.01 4.93 × 10−11 3.02 3.19 × 10−11 3.01 2.33 × 10−10 3.01 
Table 6

Convergence rates for the water depth and discharge, h and , associated with the and error norms, computed by the well-balanced DG2-SSPRK3 schemes



Meshes ()
 4.43 × 10−04 – 2.30 × 10−03 – 2.20 × 10−03 – 1.16 × 10−02 – 
 7.80 × 10−05 2.51 4.01 × 10−04 2.49 3.63 × 10−04 2.57 2.20 × 10−03 2.39 
 1.08 × 10−05 2.84 5.48 × 10−05 2.87 5.04 × 10−05 2.85 3.31 × 10−04 2.74 
 1.57 × 10−06 2.79 7.31 × 10−06 2.90 7.31 × 10−06 2.79 5.16 × 10−05 2.68 
 2.68 × 10−07 2.56 1.13 × 10−06 2.69 1.22 × 10−06 2.57 8.73 × 10−06 2.56 
 5.51 × 10−08 2.23 2.24 × 10−07 2.34 2.47 × 10−07 2.31 1.74 × 10−06 2.32 
 1.27 × 10−08 2.12 5.09 × 10−08 2.14 5.65 × 10−08 2.13 3.98 × 10−07 2.13 


Meshes ()
 4.43 × 10−04 – 2.30 × 10−03 – 2.20 × 10−03 – 1.16 × 10−02 – 
 7.80 × 10−05 2.51 4.01 × 10−04 2.49 3.63 × 10−04 2.57 2.20 × 10−03 2.39 
 1.08 × 10−05 2.84 5.48 × 10−05 2.87 5.04 × 10−05 2.85 3.31 × 10−04 2.74 
 1.57 × 10−06 2.79 7.31 × 10−06 2.90 7.31 × 10−06 2.79 5.16 × 10−05 2.68 
 2.68 × 10−07 2.56 1.13 × 10−06 2.69 1.22 × 10−06 2.57 8.73 × 10−06 2.56 
 5.51 × 10−08 2.23 2.24 × 10−07 2.34 2.47 × 10−07 2.31 1.74 × 10−06 2.32 
 1.27 × 10−08 2.12 5.09 × 10−08 2.14 5.65 × 10−08 2.13 3.98 × 10−07 2.13 

Case 3. Riemann problems

The numerical resolution of different RPs, provided in Table 7, is considered in this section. The spatial domain is given by , the discontinuity is positioned at and the acceleration of gravity is set to = 9.81 m/s2. Numerical solutions are plotted for , and and compared with exact solutions. For all the problems, we set .

Table 7

Initial condition for the RPs

LabelT
RP1 1.0 0.1614067989 0.0 0.0 1.0 0.1614067989 0.0 0.05 0.02 
RP2 1.0 0.6245627691 1.0 1.0 0.0 0.0 0.0 0.3 0.01 
RP3 1.0 1.0 0.0 0.0 1.0 2.0 0.0 0.0  
RP4 0.4 1.6 2.3 1.686811595600704 1.0 0.0 0.0 0.0 0.1835 
LabelT
RP1 1.0 0.1614067989 0.0 0.0 1.0 0.1614067989 0.0 0.05 0.02 
RP2 1.0 0.6245627691 1.0 1.0 0.0 0.0 0.0 0.3 0.01 
RP3 1.0 1.0 0.0 0.0 1.0 2.0 0.0 0.0  
RP4 0.4 1.6 2.3 1.686811595600704 1.0 0.0 0.0 0.0 0.1835 

Units in m, m2/s and s.

The first RP (RP1) is computed using two different grids, composed of 160 and 320 cells, respectively, to show the convergence to the exact solution with mesh refinement. The solution for the water surface elevation, normal discharge and shear discharge are plotted in Figure 7. It is observed that the intermediate states are accurately predicted, as well as the location of the shock wave and the rarefaction wave. As no slope limiter is used, the second order scheme exhibits a spurious spike in the cell containing the right-moving shock. However, the presence of this spike does not affect the rest of the solution. When moving from the first order scheme to the second order scheme, a significant gain in accuracy is evidenced. A detail of the solution around the contact discontinuity is also presented in Figure 7, showing that the second order scheme accurately reproduces such a discontinuous wave without oscillations, even when the slope limiter is not used.

Figure 7

RP1. Numerical solution for (top), (middle) and (bottom) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 160 (left) and 320 cells (right), compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 7

RP1. Numerical solution for (top), (middle) and (bottom) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 160 (left) and 320 cells (right), compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

The second RP (RP2) consists of an equilibrium RP, as the initial data is a discontinuous solution which satisfies the RH condition. The numerical solution, depicted in Figure 8, is computed using 160 cells. As expected, both the first order and second order schemes preserve the steady moving water equilibrium with machine precision.

Figure 8

RP2. Numerical solution for (left) and (right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 160 cells, compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 8

RP2. Numerical solution for (left) and (right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 160 cells, compared with the exact solution (black solid line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

The third RP (RP3) consists of another equilibrium RP representing an ideal shear layer from Navas-Montilla & Juez (2020). The left and right regions of the flow have a zero velocity in the -direction, but there is velocity in the transverse (shear) direction with a discontinuity at the interface. Ideally, the left and right regions of the flow should be decoupled as there is no viscosity in the equations and thus no explicit shear stresses. If a complete Riemann solver such as the (A)Roe or HLLC(S) solver were used, the exact solution would be obtained and the left and right regions would be decoupled in the numerical solution. When using the HLLS solver, which does not explicitly account for the contact wave, the shear discontinuity in is smeared out due to the numerical diffusion intrinsic to this particular solver (Navas-Montilla & Juez 2020). The numerical solution for this problem is computed using 80 and 160 cells at times and s, and is depicted in Figure 9. Only is shown in the figure, as the other variables are constant and have no interest. The exact solution and the analytical solution of the equivalent equation of the first order HLLS solver is also depicted. The latter is given by 
formula
(58)
where is the wave celerity (Navas-Montilla & Juez 2020). Figure 9 shows that the first order scheme adds a large amount of numerical diffusion, which would involve a strong mixing between the left and right regions in the context of the resolution of a shear layer. This diffusion increases as the cell size or the simulation time is increased. On the other hand, the DG2-RK3 scheme reproduces a sharper transition across the contact discontinuity (roughly two-cell wide) and shows no dependence on the simulation time (this is noticed as the DG polynomials show no discontinuities across cell interfaces). This is a clear advantage when compared to other high-order numerical methods, such as WENO-ADER schemes, which show a more smeared profile (Navas-Montilla & Juez 2020). A small overshoot is shown in the DG2-RK3 solution as no slope limiter is used. However, the performance of such a scheme is by far better than the first order version of the HLLS solver and even better than other higher order schemes.
Figure 9

RP3. Numerical solution for computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at (left) and s (right) using 80 (top) and 160 cells (bottom), compared with the exact solution (black solid line) and the analytical solution of the equivalent equation for the first order scheme (dashed line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 9

RP3. Numerical solution for computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at (left) and s (right) using 80 (top) and 160 cells (bottom), compared with the exact solution (black solid line) and the analytical solution of the equivalent equation for the first order scheme (dashed line). Units are given in m and m/s. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

The fourth RP (RP4) involves the computation of a slowly moving hydraulic jump (strong shock), usually linked to the slowly moving shock anomaly in the framework of FV schemes (Navas-Montilla & Murillo 2019). The aim of this case is to assess the robustness of the augmented DG-RK scheme without the use of a limiter and to determine whether or not this scheme suffers from the slowly moving shock anomaly. Figure 10 shows the numerical solution for h and computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 80 cells. It is observed that both schemes produce a spurious oscillation downstream the hydraulic jump, however, of a different nature. In Figure 11, the evolution in time of the average value of provided by the schemes at a cell on the left side of the initial discontinuity () is depicted. It is observed that the first order scheme produces a spurious spike of the unit discharge (up to m2/s) when the shock crosses the aforementioned cell, which corresponds to the slowly moving shock anomaly (Navas-Montilla & Murillo 2019). Contrarily, the DG2-SSPRK3 does not show this anomaly, being the oscillations only produced by the lack of a limiter. This evidences that the proposed DG2-SSPRK3 scheme is free from the slowly moving shock anomaly.

Figure 10

RP4. Numerical solution for h (left) and (right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 80 cells, compared with the exact solution (black solid line). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 10

RP4. Numerical solution for h (left) and (right) computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 (purple) at s using 80 cells, compared with the exact solution (black solid line). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 11

RP4. Evolution in time of the average value of at a cell on the left side of the initial discontinuity () computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 scheme (purple). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

Figure 11

RP4. Evolution in time of the average value of at a cell on the left side of the initial discontinuity () computed by the well-balanced first order Godunov's scheme (green) and DG2-SSPRK3 scheme (purple). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2020.206.

The reason for this is that FV schemes inherently enforce the width of shock waves due to the cell-averaged representation of data, whereas DG schemes do not. When cell-averaging the data around the shock, non-physical values of the conserved variables may appear in the case of a non-monotone (non-linear) Hugoniot locus (Navas-Montilla & Murillo 2019), which eventually produce the observed spike in the SWE. Contrarily, DG methods are based on a modal/nodal reconstruction of pointwise data, without averaging, and therefore do not exhibit this problem. Nevertheless, this behavior should be further investigated in future works.

CONCLUDING REMARKS

In this work, the application of augmented Riemann solvers in the framework of the DG method is addressed. An explicit DG-SSPRK3 scheme based on augmented Riemann solvers is constructed. The proposed scheme is applied to the resolution of hyperbolic conservation laws with geometric source terms, in particular, to the Burgers' equation with geometric source term and to the SWE with bed topography.

The use of augmented Riemann solvers (e.g. the HLLS solver) proves to be very useful and effective for the preservation of equilibrium states. Thanks to the definition of the RP considering the contribution of the source term at cell interfaces, the exact RH condition is satisfied at cell interfaces by the Riemann solver without requiring extra flux corrections. It is proved that the DG spatial operator becomes nil under equilibrium, provided (i) the use of an exact quadrature rule for the surface and volume integrals and (ii) the reconstructed data satisfy the discrete equilibrium. Such a discrete equilibrium property of the scheme is proved with the independence of the number of spatial dimensions.

For a proper analysis of the fundamental solutions obtained by the combination of augmented solvers with the DG method, the scheme herein presented includes neither extra flux corrections nor slope limiting algorithms. Nevertheless, the incorporation of additional numerical fixes (i.e. slope limiters, wet-dry treatment, etc.) to address realistic problems of engineering and environmental interest is straightforward. Furthermore, the application of this scheme to higher spatial dimensions is also straightforward thanks to the rotational invariance property of the systems (i.e. 1D RPs are always solved, with the independence of the number of spatial dimensions).

The numerical results evidence that the proposed scheme achieves the prescribed convergence rates and preserves the equilibrium states of relevance with machine precision, for the problems considered in this work. Convergence rates have been assessed up to third order of accuracy, as this is the limiting order due to the SSPRK3 time integrator. Arbitrary order numerical schemes can be constructed using the DG-ADER approach (Dumbser & Munz 2005; Fambri et al. 2017), where augmented solvers can be applied as presented here for the resolution of the Derivative RP using the LFS solver in Navas-Montilla & Murillo (2016). The results also evidence that the proposed DG2-SSPRK3 scheme is free from the slowly moving shock anomaly.

ACKNOWLEDGEMENTS

This work was funded by the Spanish Ministry of Science and Innovation under the research project PGC2018-094341-B-I00. This work has also been partially funded by Gobierno de Aragón through Fondo Social Europeo (Feder 2014-2020 ‘Construyendo Europa desde Aragón’).

SUPPLEMENTARY MATERIAL

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/hydro.2020.206.

REFERENCES

REFERENCES
Audusse
E.
,
Bouchut
F.
,
Bristeau
M. O.
,
Klein
R.
&
Perthame
B. T.
2004
A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows
.
SIAM J. Sci. Comput.
25
,
2050
2065
.
Bermudez
A.
&
Vázquez-Cendón
M. E.
1994
Upwind methods for hyperbolic conservation laws with source terms
.
Comput. Fluids
23
,
1049
1071
.
Caleffi
V.
,
Valiani
A.
&
Li
G.
2015
A comparison between bottom-discontinuity numerical treatments in the DG framework
.
Appl. Math. Model.
40
(
17–18
),
7516
7531
.
Catella
M.
,
Paris
E.
&
Solari
L.
2008
Conservative scheme for numerical modeling of flow in natural geometry
.
J. Hydraul. Eng. – ASCE
134
(
6
),
736
748
.
Caviedes-Voullième
D.
,
Gerhard
N.
,
Sikstel
A.
&
Müller
S.
2020
Multiwavelet-based mesh adaptivity with discontinuous Galerkin schemes: exploring 2D shallow water problems
.
Adv. Water Resour.
138
,
103559
.
Chertock
A.
,
Dudzinski
M.
,
Kurganov
A.
&
Lukácová-Medvidóvá
M.
2017
Well-balanced schemes for the shallow water equations with Coriolis forces
.
Numer. Math.
138
(
4
),
939
973
.
Chinnayya
A.
,
LeRoux
A. Y.
&
Seguin
N.
2004
A well-balanced numerical scheme for the approximation of the shallow-water equations with topography: the resonance phenomenon
.
Int. J. Finite
1
,
1
33
.
Cockburn
B.
&
Shu
C.-W.
1989
TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework
.
Math. Comput.
52
,
411
435
.
Cockburn
B.
&
Shu
C.-W.
1990
TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case
.
Math. Comput.
54
,
545
581
.
Cunge
J. A.
,
Holly
F. M.
&
Verwey
A.
1980
Practical Aspects of Computational River Hydraulics
.
Pitman
,
London
,
UK
.
Dumbser
M.
&
Munz
C. D.
2005
ADER discontinuous Galerkin schemes for aeroacoustics
.
C. R. Mécanique
333
,
683
687
.
García-Navarro
P.
&
Vázquez-Cendón
M. E.
2000
On numerical treatment of the source terms in the shallow water equations
.
Comput. Fluids
29
,
951
979
.
Gordillo
G.
,
Morales-Hernández
M.
&
García-Navarro
P.
2019
Finite volume model for the simulation of 1D unsteady river flow and water quality based on the WASP
.
J. Hydroinform
.
22
(
2
),
327
345
.
Haleem
D. A.
,
Kesserwani
G.
&
Caviedes-Voullième
D.
2015
Haar wavelet-based adaptive finite volume shallow water solver
.
J. Hydroinform.
17
,
857
873
.
Kesserwani
G.
,
Ghostine
R.
,
Vazquez
J.
,
Ghenaim
A.
&
Mosé
R.
2008
Application of a second-order Runge-Kutta discontinuous Galerkin scheme for the shallow water equations with source terms
.
Int. J. Numer. Meth. Fluids
56
,
805
821
.
Lee
S. H.
&
Wright
N. G.
2010
Simple and efficient solution of the shallow water equations with source terms
.
Int. J. Numer. Meth. Fluids
63
,
313
340
.
Manzanero
J.
,
Rubio
G.
,
Kopriva
D. A.
,
Ferrer
E.
&
Valero
E.
2020
A free-energy stable nodal discontinuous Galerkin approximation with summation-by-parts property for the Cahn-Hilliard equation
.
J. Comput. Phys.
340
,
200
242
.
Martínez-Aranda
S.
,
Murillo
J.
&
García-Navarro
P.
2019
A comparative analysis of capacity and non-capacity formulations for the simulation of unsteady flows over finite-depth erodible beds
.
Adv. Water Resour.
130
,
91
112
.
Reed
W. H.
&
Hill
T. R.
1973
Triangular mesh methods for the neutron transport equation
. In:
National Topical Meeting on Mathematical Models and Computational Techniques for Analysis of Nuclear Systems
.
Los Alamos Scientific Laboratory, NM. Available from https://www.osti.gov/servlets/purl/4491151
.
Sanders
B. F.
,
Jaffe
D. A.
&
Chu
A. K.
2003
Discretization of integral equations describing flow in nonprismatic channels with uneven beds
.
J. Hydraul. Eng. – ASCE
129
(
3
),
235
244
.
Xing
Y.
&
Shu
C. W.
2006
A new approach of high order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms
.
Comput. Phys.
1
,
100
134
.

Supplementary data