Well-Balanced High-Order Discontinuous Galerkin Methods for Systems of Balance Laws

Abstract

This work introduces a general strategy to develop well-balanced high-order Discontinuous Galerkin (DG) numerical schemes for systems of balance laws. The essence of our approach is a local projection step that guarantees the exactly well-balanced character of the resulting numerical method for smooth stationary solutions. The strategy can be adapted to some well-known different time marching DG discretisations. Particularly, in this article, Runge–Kutta DG and ADER DG methods are studied. Additionally, a limiting procedure based on a modified WENO approach is described to deal with the spurious oscillations generated in the presence of non-smooth solutions, keeping the well-balanced properties of the scheme intact. The resulting numerical method is then exactly well-balanced and high-order in space and time for smooth solutions. Finally, some numerical results are depicted using different systems of balance laws to show the performance of the introduced numerical strategy.

1. Introduction

We consider one-dimensional systems of balance laws of the form

$${\partial }_{t}U+{\partial }_{x}F\left(U\right)=S\left(U\right)H^{\prime }+G\left(U\right).$$

(1)

Here, $U(x,t)$ is the state vector taking values in the space of admissibly states $\Omega \subset {\mathbb{R}}^M$; $F:\Omega \to {\mathbb{R}}^M$ is the flux; $S\left(U\right):\Omega \to {\mathbb{R}}^M,$ is a function that takes into account the geometrical source terms that are multiplied by the known function $H:\mathbb{R}\to \mathbb{R}$; and $G\left(U\right)$ contains algebraic source terms. We suppose that the system is hyperbolic, and therefore the jacobian matrix of the flux has M distinct real eigenvalues.

The system (1) often admits steady-state solutions that, in most cases, are non-trivial. Moreover, these solutions usually convey significant physical meaning. A paradigmatic example of this is the lake-at-rest solution of the well-known shallow-water equations. The correct preservation of these solutions plays an essential role in any simulation. Well-balanced (or exactly well-balanced) numerical schemes aim to preserve all or at least a relevant set of specific steady-state solutions up to machine precision. A precise definition of well-balanced method is given in Section 2.

An essential advantage of well-balanced schemes is that they can accurately resolve small perturbations to steady-state solutions with relatively coarse meshes. Indeed, if a numerical method cannot preserve stationary solutions or recover a stationary solution after a perturbation, the oscillations introduced may destroy the numerical solution for long integration times. Additionally, the numerical solution may be lost entirely if the perturbation is in the same order as the discretisation errors introduced by a non well-balanced method. Then, a finer mesh or a higher-order numerical method may be considered. Still, in this case, the computational effort can quickly become excessive. The ability of the numerical scheme to correctly preserve stationary solutions becomes then a matter of paramount importance.

The first effort to preserve stationary solution was originally introduced by Bermúdez and Vázquez-Cendón in [1] for the sallow water system, where the C-property was defined as the ability of a numerical scheme to solve a set of stationary solutions exactly or, at least, with enhanced accuracy. Since then, multiple works have been published on well-balanced numerical methods [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]. It is of interest to highlight the recent work in [35], where the authors develop high order exactly well-balanced methods for systems of conservation laws for finite volume method with both continuous and discontinuous source terms. A similar strategy is presented in [36] for finite-difference methods. This paper extends the numerical strategy introduced within the framework of finite volume methods in [35] to Discontinuous Galerkin (DG) methods, describing a general approach to preserve stationary solutions for general systems of conservation laws.

The DG methods were introduced by Reed and Hill in [37] and were later extended to nonlinear hyperbolic systems of conservation laws by Cockburn and Shu in a series of publications [38,39,40,41]. Since then, several authors have proposed DG numerical schemes that are exactly well-balanced for a set of stationary solutions but they mainly focus on a particular system, most commonly the shallow water equations (see [42,43,44,45,46,47,48,49,50,51,52,53]).

Many hyperbolic systems of conservation laws are frequently numerically approximated with Discontinuous Galerkin methods as they allow to scale the numerical order of accuracy of the scheme quickly, thus providing an excellent resolution with coarse meshes. That makes these schemes especially well-suited when the computational effort is an important issue. A natural approach when discretising hyperbolic systems would be to use the methods of lines, in which all but the time dimension is approximated by the DG approach. The resulting time dimensional system of ODE can be then discretised using explicit methods such as Runge–Kutta [54,55,56,57,58,59,60]. A more sophisticated approach that allows reaching arbitrary high order in space and time in a fully implicit manner consist on ADER (Arbitrary high order DErivative Riemann problem) approach. The ADER technique was first developed by Toro and Titarev for finite volume methods [61,62,63,64] and later extended to DG methods [55,65,66]. Since then, it has been widely used in order to solve many hyperbolic systems of conservation laws see [67,68,69,70,71,72]. The recent ADER approach for DG methods is based on a predictor-corrector scheme that computes an approximated predictor solution through a fixed point algorithm locally in each element. With the DG method providing arbitrary high order in space, the combination of both techniques allows one-step fully discrete arbitrary high order in space and time numerical schemes.

The DG and ADER-DG methods have been successfully applied to geophysical flows. In the framework of the one-layer shallow water model, we highlight the work in [68,69,70,73] among others, and in the multilayer shallow water case in [74,75]. More recently, advances on hyperbolic reformulations of nonlinear dispersive shallow water systems that make use of the DG or ADER-DG techniques are found on [67,76,77,78].

A recurrent problem when high order methods are used is the appearance of spurious oscillations at discontinuities that appear when the hyperbolic system of conservation laws is considered. To avoid that, several limiting techniques exist. For instance, a multi-dimensional and optimal-order detector (MOOD) can be used (see in [79]). Alternatively, other limiter strategies are the minmod type total variation bounded (TVB) limiter [38,80], the moment-based limiter [81] and the improved moment limiter [82]. Another approach to limiting is the one based on WENO methodology, mainly developed by Shu, Qiu, and Zhu et al. (see in [83,84,85]). The WENO limiter approach uses a convex combination of all available neighbour polynomials, assigning essentially zero weights to stencils containing discontinuities. This work proposes a limiter procedure based on this method but takes advantage of a well-balanced reconstruction operator that preserves the well-balanced properties of the numerical method.

This paper is organised as follows. In Section 2, a brief description of the considered numerical schemes are performed, including Runge–Kutta and ADER discretisations in time. Section 3 describes the followed well-balanced methodology proposed for the numerical schemes studied. Next, in Section 4, a limiter procedure that preserves the well-balanced features of the numerical schemes is proposed. Then, in Section 5 some numerical tests are provided for the Shallow water, Burger, Euler with gravity and the hyperbolic dispersive Gavrilyuk–Favrie equations. Finally, some conclusions are given in Section 6.

2. Discontinuous Galerking Method

We briefly recall the DG method to discretise the general system (1) for one-spatial domains. We consider a set of element cells ${\Omega }_i=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$, with constant length $\Delta x.$ The numerical solution of the PDE (1) at a given time t at each subset ${\Omega }_i$ is represented by $u_h(x,t)$ and is described in terms of piece-wise polynomials of degree N on the spatial direction. We shall denote by ${\mathcal{U}}_h$ the space of piece-wise polynomials up to degree N so that $u_h{(x,t)|}_{{\Omega }_i}\in {\mathcal{U}}_h,\forall i$.

An orthogonal nodal basis is defined by the Lagrange interpolation polynomials over the $(N+1)$ Gauss–Legendre quadrature nodes on the cell element ${\Omega }_i$. We stress that the piece-wise solution $u_h$ may be discontinuous across elements, allowing jump discontinuities at cell interfaces.

Given the element ${\Omega }_i$, the discrete solution $u_h$ projected on the nodal spatial basis polynomials ${\Phi }_l\left(x\right)$ reads

$$u_h(x,t)=\sum _l{\widehat{u}}_{i,l}\left(t\right){\Phi }_l\left(x\right){\widehat{u}}_{i,l}\left(t\right){\Phi }_l\left(x\right),\mathrm{for}x\in {\Omega }_i,$$

(2)

with ${\widehat{u}}_{i,l}\left(t\right)$ being the unknown degrees of freedom, and where the Einstein summation convention over two repeated indices has been considered. In this work, we suppose that $H\left(x\right)$ is a given function and it is known at every point of the domain. In the subsequent notation, we omit the space dependency of ${\Phi }_l\left(x\right)$ for the sake of simplicity.

The DG methods results then from multiplying the PDE system (1) with a polynomial test function ${\Phi }_k\in {\mathcal{U}}_h$ and integrating over the spatial element ${\Omega }_i$. That leads to the following weak formulation where we seek $u_h(x,t)\in {\mathcal{U}}_h$ such that, for all ${\Omega }_i$,

$$\begin{array}{c}\hfill {\int }_{{\Omega }_i}{\Phi }_k{\partial }_t{u}_{h}dx+{\int }_{{\Omega }_i}{\Phi }_k{\partial }_{x}F\left(u_h\right)dx={\int }_{{\Omega }_i}{\Phi }_k\left(S\left(u_h\right)H^{\prime }+G\left(u_h\right)\right)dx,\end{array}$$

(3)

holds for all test functions ${\Phi }_k\left(x\right)\in {\mathcal{U}}_h$. The dependence on x and t has been dropped for the sake of compactness. Using (2), the flux function is integrated by parts and the weak formulation (3) can be rewritten as

$$\begin{array}{c}\left({\int }_{{\Omega }_i}{\Phi }_k^{2}dx\right)\frac{d}{dt}{\widehat{u}}_{i,k}={\Phi }_k\left(x_{i-1/2}\right)F_{i-1/2}-{\Phi }_k\left(x_{i+1/2}\right)F_{i+1/2}\hfill \\ \hfill +{\int }_{{\Omega }_i^{\circ }}F\left(u_h\right){\Phi }_k^{\prime }dx+{\int }_{{\Omega }_i^{\circ }}{\Phi }_k\left(S\left(u_h\right)H^{\prime }+G\left(u_h\right)\right)dx,\end{array}$$

(4)

where the first integral contains the so-called element mass matrix, which is diagonal since the DG polynomial basis is orthogonal. Likewise, ${\Omega }_i^{\circ }$ denotes the set of all interior points of ${\Omega }_i$. The second integral accounts for the smooth part of the flux, while the last integral accounts for the presence of the continuous source terms. Finally, $F_{i-1/2}-F_{i+1/2}$ describes the jump across element interfaces given by a consistent numerical flux of the homogeneous system

$$F_{i+1/2}=\mathbb{F}\left(u_{h,i+1/2}^{-},u_{h,i+1/2}^{+}\right),$$

(5)

where $u_{h,i+1/2}^{-}$ and $u_{h,i+1/2}^{+}$ are the boundary extrapolated states at the left and right of the interface $x_{i+1/2}$, respectively. Note that the numerical fluxes are multiplied by the evaluation of the Lagrange base at the cell interfaces, ${\Phi }_k\left(x_{i\pm 1/2}\right)$. In this work, the standard HLL Riemann solver is used. However, the numerical methods described in this work can be easily extended to more accurate Riemann solvers such as the Polynomial Viscosity Methods or Rational Viscosity Methods (see [86,87,88,89]).

Definition 1.

Given a continuous stationary solution U of system (1), the semi-discrete numerical scheme (4) is said to be well-balanced (or exactly well-balanced) if the set $\{{\widehat{u}}_{i,l}=U\left(x_{i,l}\right)\}$ is an equilibrium of the ODE system, where $\left\{x_{i,l}\right\}$ is the set of nodal spatial points. Note that this definition mimics the one given in [35] or [31] for finite volume methods.

It is clear that the semi-discrete method (4), as it is formulated, is not well-balanced in general. The following section will describe a procedure to define DG well-balanced methods.

The semi-discrete numerical scheme (4) can be then approximated with any time-marching numerical method. This work uses classical Runge–Kutta methods and the family of one-step ADER-DG methods based on a predictor-corrector scheme. In the next Subsection, a brief description of the aforementioned time discretisation methods is given. The interested reader can be referred to [38,39,40,41,80,90] for DG Runge–Kutta methods or [67,68,69,70,71,72] for ADER-DG methods.

2.1. Runge–Kutta Time Discretization

We first remark that due to the choice of an orthogonal DG polynomial basis, the mass-matrix appearing in the semi-discrete system (4) is diagonal, and that leads us to write (4) for each element ${\Omega }_i$ and degree of freedom ${\widehat{u}}_{i,k}$ in the more compact form

$${\partial }_t{\widehat{u}}_{i,k}\left(t\right)={\displaystyle \frac{1}{{\omega }_k}}L\left({\widehat{u}}_{i,k}\left(t\right)\right),{\omega }_k={\int }_{{\Omega }_i}{\Phi }_k^{2}dx,$$

(6)

where $L\left({\widehat{u}}_{i,k}\left(t\right)\right)$ is the spatial discretization operator. To discretise the time variable $t,$ total variation diminishing (TVD) schemes have been considered in this work.

For instance, the third order Runge–Kutta scheme applied to the time interval $[t^n,t^n+\Delta t],$ $\Delta t$ being the time-step, reads

$$\begin{array}{cc}\hfill {\widehat{u}}_{i,k}^{\left(1\right)}& =L\left({\widehat{u}}_{i,k}\left(t^n\right)\right),\hfill \\ \hfill {\widehat{u}}_{i,k}^{\left(2\right)}& =L({\widehat{u}}_{i,k}\left(t^n\right)+\frac{\Delta t}{2}{\widehat{u}}_{i,k}^{\left(1\right)}),\hfill \\ \hfill {\widehat{u}}_{i,k}^{\left(3\right)}& =L({\widehat{u}}_{i,k}\left(t^n\right)-\Delta t{\widehat{u}}_{i,k}^{\left(1\right)}+2\Delta t{\widehat{u}}_{i,k}^{\left(2\right)}),\hfill \\ \hfill {\widehat{u}}_{i,k}^{n+1}& ={\widehat{u}}_{i,k}\left(t^n\right)+\frac{\Delta t}{6}({\widehat{u}}_{i,k}^{\left(1\right)}+4{\widehat{u}}_{i,k}^{\left(2\right)}+{\widehat{u}}_{i,k}^{\left(3\right)}).\hfill \end{array}$$

(7)

2.2. ADER Time Discretization

The predictor-corrector ADER-DG approach leads to a high-order one-step method in both space and time. Here, the high-order accuracy in time is reached with a local space and time predictor guess. The spatial-time predictor $q_h(x,t)$ is computed based on a weak formulation of the PDE. Given the known solution $u_h(x,t^n)$ at time $t^n$, a Cauchy problem in the cell element is solved without considering interactions between the neighbour elements. The predictor solution is then used in (4) to compute the numerical solution $u_h(x,t^{n+1})$ at the subsequent time step.

Here, we briefly recall the numerical procedure. Given a space-time control volume ${\Omega }_i\times [t^n,t^{n+1}]$, the predictor solution $q_h$ is projected onto the local space-time basis,

$$\left\{\begin{array}{c}{q}_h(x,t)=\sum _l{\theta }_l(x,t){\widehat{q}}_l^i{\theta }_l(x,t){\widehat{q}}_l^i,\mathrm{for}x\in {\Omega }_i,t\in (t^n,t^{n+1}],\hfill \\ q_h(x,t^n)=u_h(x,t^n),\hfill \end{array}\right.$$

(8)

where again the Einstein notation has been used. If the system (1) is multiplied by the space-time test function ${\theta }_k\equiv {\theta }_k(x,t)$ and integrated over ${\Omega }_i\times [t^n,t^{n+1}]$, it yields

$$\begin{array}{c}{\int }_{{\Omega }_i}{\theta }_k^{n+1}{q}_h(x,t^{n+1})dx-{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i}{\partial }_t{\theta }_k{q}_h(x,t)dxdt={\int }_{{\Omega }_i}{\theta }_k^n{q}_h(x,t^n)dx\hfill \\ \hfill -{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i}{\theta }_k\left({\partial }_{x}F\left(q_h\right)-S\left(q_h\right)H^{\prime }-G\left(q_h\right)\right)dxdt.\end{array}$$

(9)

Using Equation (8), Equation (9) becomes a nonlinear system for the unknowns ${\widehat{q}}_l^i$ that can be solved via a fixed point algorithm detailed in [65,91]. Notice that (9) is completely local, without any interaction with neighbouring cells. The iterative method converges in $N+1$ steps at most for linear homogeneous systems [92].

Finally, the one-step high-order in both space and time ADER-DG method, after replacing $u_h(x,t)$ by $q_h(x,t)$ in the right-hand side of (4), and integrating in the time interval $[t^n,t^{n+1}],$ reads

$$\begin{array}{c}\left({\int }_{{\Omega }_i}{\Phi }_k^{2}dx\right)\left({\widehat{u}}_{i,k}^{n+1}-{\widehat{u}}_{i,k}^n\right)={\int }_{t^n}^{t^{n+1}}{\Phi }_k\left(x_{i-1/2}\right)F_{i-1/2}-{\Phi }_k\left(x_{i+1/2}\right)F_{i+1/2}dt\hfill \\ \hfill +{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i^{\circ }}F\left(q_h\right){\partial }_x{\Phi }_{k}dxdt+{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i^{\circ }}{\Phi }_k\left(S\left(q_h\right)H^{\prime }+G\left(q_h\right)\right)dxdt,\end{array}$$

(10)

where

$$F_{i+1/2}=\mathbb{F}\left(q_{h,i+1/2}^{-},q_{h,i+1/2}^{+}\right).$$

Remark 1.

DG methods are restricted to the Courant–Friedrichs–Lewy (CFL) condition

$$\Delta t\le \frac{CFL}{2N+1}\frac{\Delta x}{|{\lambda }_{max}|},CFL\in (0,1).$$

(11)

3. Well-Balanced Technique

In this Section, we detail a well-balanced technique for the Runge–Kutta DG methods (4)–(7) and the ADER-DG method (9) and (10).

3.1. Well-Balanced Runge–Kutta DG Methods

We propose now a spatial DG discretization method that exactly preserves stationary solutions. Later, we will discuss the well-balanced technique for the chosen time marching discretization. As in [35], we begin by computing a suited stationary solution in each element ${\Omega }_i$ at each time step $t^n$:

• Seek for the stationary solution $u_i^{*}\left(x\right)$ of the minimization problem in ${\Omega }_i$

  $$\{\begin{array}{}\mathrm{(12a)}& {\partial }_{x}F\left(u_i^{*}\left(x\right)\right)=S\left(u_i^{*}\left(x\right)\right)H^{\prime }\left(x\right)+G\left(u_i^{*}\left(x\right)\right),x\in {\Omega }_i\hfill \mathrm{(12b)}& {\displaystyle {\displaystyle \frac{1}{\Delta x}}{\int }_{{\Omega }_i}{u}_i^{*}\left(x\right)dx={\displaystyle \frac{1}{\Delta x}}{\int }_{{\Omega }_i}{u}_h(x,t^n)dx,}\hfill \mathrm{(12c)}& {\displaystyle \mathrm{that} \mathrm{minimizes}{\int }_{{\Omega }_i}{(u_i^{*}\left(x\right)-u_h(x,t^n))}^{2}dx.}\hfill \end{array}$$

  The previous problem is in general difficult. Nevertheless if $u_i^{*}\left(x\right)$ depends on a set of parameters, problem (12) reduces, in general, to a non-linear system of equations if the integrals in (12b) are approximated by quadrature formulas defined in terms of the Gaussian nodal points of the DG polynomials. Moreover, in some situations, like the ones considered below, (12a) and (12b) are enough to determine uniquely the local stationary solution. For instance, if $u_i^{*}\left(x\right)$ is of the form

  $$u_i^{*}\left(x\right)=c_i+s\left(x\right),x\in {\Omega }_i,c_i\in \mathbb{R},s\left(x\right)\mathrm{known}\mathrm{function},$$

  (13)

  then the integration of (12b) by using quadrature rules leads to

  $$c_i=\sum _l{\omega }_l\left(u_h(x_l,t^n)-s\left(x_l\right)\right),$$

  (14)

  where ${\omega }_l,x_l$ are the Gaussian quadrature rule weights and nodes respectively. This quadrature rule must be of the same order of the DG space approximation at the element ${\Omega }_i.$ Similarly, for stationary solutions $u_i^{*}\left(x\right)$ of the form

  $$u_i^{*}\left(x\right)=c_{i}s\left(x\right),x\in {\Omega }_i,c_i\in \mathbb{R},s\left(x\right)\mathrm{known}\mathrm{function},$$

  (15)

  the constant $c_i$ can be simply computed from (12b) and yields

  $$c_i={\displaystyle \frac{{\sum }_l{\omega }_l{u}_h(x_l,t^n)}{{\sum }_l{\omega }_{l}s\left(x_l\right)}}.$$

  (16)
• Compute the fluctuation DG polynomial:

  $${v_h}_{{|}_{{\Omega }_i}}(x,t)={u_h}_{{|}_{{\Omega }_i}}(x,t)-u_{i,h}^{*}\left(x\right),x\in {\Omega }_i,t\in [t^n,t^{n+1}],$$

  (17)

  where $u_{i,h}^{*}$ is the projection of $u_i^{*}$ into ${\mathcal{U}}_h$. Note that (12c) guarantees that $v_h(x,t^n)$ is minimal in $L^2$-norm.
• To preserve stationary solutions, the numerical method (4) is modified as follows:

  $$\begin{array}{c}\left({\int }_{{\Omega }_i}{\Phi }_k^{2}dx\right)\frac{d}{dt}{\widehat{u}}_{i,k}={\int }_{{\Omega }_i^{\circ }}\left(F\left(u_h\right)-F\left(u_{i,h}^{*}\right)\right){\partial }_x{\Phi }_{k}dx\hfill \\ +{\Phi }_k\left(x_{i-1/2}\right)\left({\tilde{F}}_{i-1/2}-F_{i-1/2}^{*}\right)+{\Phi }_k\left(x_{i+1/2}\right)\left(F_{i+1/2}^{*}-{\tilde{F}}_{i+1/2}\right)\\ \hfill +{\int }_{{\Omega }_i^{\circ }}{\Phi }_k\left(S\left(u_h\right)-S\left(u_{i,h}^{*}\right)\right){\widehat{H}}_{i,l}{\partial }_x{\Phi }_{l}dx+{\int }_{{\Omega }_i^{\circ }}{\Phi }_k\left(G\left(u_h\right)-G\left(u_{i,h}^{*}\right)\right)dx,\end{array}$$

  (18)

  where ${\widehat{H}}_{i,l}$ are the coordinates of the projection of $H\left(x\right)$ into ${\mathcal{U}}_h$, and

  $${\tilde{F}}_{i+1/2}=\mathbb{F}\left(u_i^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{-},u_{i+1}^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{+}\right),$$

  (19)

  is the usual first-order consistent numerical flux evaluated at the states

  $$u_i^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{-},u_{i+1}^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{+},$$

  (20)

  where $u_i^{*}$ and $u_{i+1}^{*}$ are the stationary solutions computed at the step 1 at cells ${\Omega }_i$ and ${\Omega }_{i+1}$, respectively, and $v_{h,i+1/2}^{-}={v_h}_{{|}_{{\Omega }_i}}\left(x_{i+1/2}\right)$ and $v_{h,i+1/2}^{+}={v_h}_{{|}_{{\Omega }_{i+1}}}\left(x_{i+1/2}\right)$ are the extrapolated values of the fluctuation computed at step 2 onto the left and right cell interfaces respectively. Finally, $F_{i+1/2}^{*}$ denotes the evaluation of the physical flux function at the stationary solution at the intercells:

  $$F_{i+1/2}^{*}=F\left(u_i^{*}\left(x_{i+1/2}\right)\right).$$

Remark 2.

The procedure to compute a stationary solution (13)–(16) can be generalized for some other different stationary solutions, or even for a k-parameter family of stationary solutions as described in [35],

$$u_h^{*}(x;c_{i,1},\dots ,c_{i,k}),c_{i,r}\in \mathbb{R}.$$

Remark 3.

We stress that, as $u_i^{*}$ is a solution of (12), then (18) is of the same order of accuracy than (4) (see in [35] for more details).

Theorem 1.

The semi-discrete numerical method (18) is well-balanced for a continuous global stationary solution $u^{*}$, provided that every $u_i^{*}$ defined by (12) verifies $u_{{|}_{{\Omega }_i}}^{*}=u_i^{*}$ for each ${\Omega }_i$, in the sense that ${\widehat{u}}_{i,k}^{*}$, coordinates of the projection of $u^{*}$ into ${\mathcal{U}}_h$, is an equilibrium of the ODE system (18).

Proof. 

We consider as initial condition ${\widehat{u}}_{i,k}\left(0\right)={\widehat{u}}_{i,k}^{*}$, that implies

$${u_h}_{{|}_{{\Omega }_i}}(x,0)=u_{i,h}^{*}\left(x\right),x\in {\Omega }_i,$$

and the three volume integrals are zero. In the same way, the fluctuations ${v_h}_{{|}_{{\Omega }_i}}(x,0)={u_h}_{{|}_{{\Omega }_i}}(x,0)-u_{i,h}^{*}\left(x\right)=0$, therefore the extrapolated values onto the element interfaces, $v_{h,i-1/2}^{-}\left(0\right)$ and $v_{h,i-1/2}^{+}\left(0\right),$ are exactly zero as well. The continuity of $u^{*}$ implies that $u_i^{*}\left(x_{i-1/2}\right)=u_{i-1}^{*}\left(x_{i-1/2}\right)=u^{*}\left(x_{i-1/2}\right)$ and the consistency of the numerical flux implies that

$$\begin{array}{cc}\hfill & {\tilde{F}}_{i-1/2}-F_{i-1/2}^{*}=\hfill \\ \hfill & =\mathbb{F}\left(u_{i-1,h}^{*}\left(x_{i-1/2}\right)+v_{h,i-1/2}^{-}\left(0\right),u_{i,h}^{*}\left(x_{i-1/2}\right)+v_{h,i-1/2}^{+}\left(0\right)\right)-F\left(u_i^{*}\left(x_{i-1/2}\right)\right)\hfill \\ \hfill & =\mathbb{F}\left(u^{*}\left(x_{i-1/2}\right),u^{*}\left(x_{i-1/2}\right)\right)-F\left(u^{*}\left(x_{i-1/2}\right)\right)=0.\hfill \end{array}$$

The same reasoning can be applied to prove that $F_{i+1/2}^{*}-{\tilde{F}}_{i+1/2}=0.$ We conclude that ${\widehat{u}}_{i,k}^{*}$ in an equilibrium of the ODE system given by the semi-discrete numerical method (18). □

Once the semi-discrete well-balanced method (18) has been proposed, the Cauchy problem

$$\left\{\begin{array}{c}{\partial }_t{\widehat{u}}_{i,k}={\displaystyle \frac{1}{{\omega }_k}}\tilde{L}\left({\widehat{u}}_{i,k}\right),t\in [t^n,t^{n+1}]\hfill \\ {\widehat{u}}_{i,k}\left(0\right)={\widehat{u}}_{i,k,0},\hfill \end{array}\right.$$

(21)

for an autonomous ODE system has to be solved. Here $\tilde{L}$ is given by the right-hand side of (18). Clearly, a stationary solution $u_h^{*}$ is solution of (21) if and only if

$$\tilde{L}\left(u_h^{*}\right)=0.$$

(22)

Therefore, given a time discretization of the form

$${\widehat{u}}_{i,k}^{n+1}={\widehat{u}}_{i,k}^n+\Delta t\phi ({\widehat{u}}_{i,k}^n,\Delta t),$$

for (21), the final discrete method is well-balanced for the stationary solution $u_h^{*}$ if and only if

$$\phi (u_h^{*},\Delta t)=0.$$

(23)

Finally, it is clear that (22) implies (23) for any Runge–Kutta method as the ones considered in this work.

3.2. Well-Balanced ADER DG Methods

To design a well-balanced method based on the ADER time discretization, we follow a similar strategy as for the well-balanced DG Runge–Kutta method. Given an element ${\Omega }_i,$ a time-step $t^n,$ and the known numerical solution $u_h(x,t^n):$

• First, we seek for the discrete stationary solution $u_i^{*}\left(x\right)$ of the Cauchy problem in ${\Omega }_i$, as in the first step of the previous Section 3.1.
• Then, we define the fluctuation DG polynomial $v_h$

  $$\left\{\begin{array}{c}{v_h}_{{|}_{{\Omega }_i}}(x,t)={q_h}_{{|}_{{\Omega }_i}}(x,t)-u_{i,h}^{*}\left(x\right),x\in {\Omega }_i,t\in (t^n,t^{n+1}],\hfill \\ {q_h}_{{|}_{{\Omega }_i}}(x,t^n)={u_h}_{{|}_{{\Omega }_i}}(x,t^n),x\in {\Omega }_i,\hfill \\ {v_h}_{{|}_{{\Omega }_i}}(x,t^n)={u_h}_{{|}_{{\Omega }_i}}(x,t^n)-u_{i,h}^{*}\left(x\right),x\in {\Omega }_i,\hfill \end{array}\right.$$

  (24)

  and the local space-time predictor given in (9) is modified as follows:

  $$\begin{array}{cc}\hfill {\int }_{{\Omega }_i}& {\theta }_k^{n+1}{v}_h(x,t^{n+1})dx-{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i}{\partial }_t{\theta }_k{v}_h(x,t)dxdt\hfill \\ \hfill & ={\int }_{{\Omega }_i}{\theta }_k^n{v}_h(x,t^n)dx-{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i}{\theta }_k\left({\partial }_{x}F\left(q_h\right)-{\partial }_{x}F\left(u_{i,h}^{*}\right)\right)\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i}{\theta }_k\left(\left(S\left(q_h\right)-S\left(u_{i,h}^{*}\right)\right){\widehat{H}}_{i,l}{\partial }_x{\Phi }_l+G\left(q_h\right)-G\left(u_{i,h}^{*}\right)\right)\hfill \end{array}$$

  (25)
  Once the nonlinear system for the unknown ${\widehat{v}}_l^i$ has been solved via a fixed point algorithm, the space-time predictor $q_h(x,t)$ is recovered as

  $${q_h}_{{|}_{{\Omega }_i}}(x,t)={v_h}_{{|}_{{\Omega }_i}}(x,t)+u_{i,h}^{*}\left(x\right),x\in {\Omega }_i,t\in [t^n,t^{n+1}].$$
• Finally, the one-step ADER-DG method described in (10) is modified as follows:

  $$\begin{array}{cc}\hfill & \left({\int }_{{\Omega }_i}{\Phi }_k^{2}dx\right)\left({\widehat{u}}_{i,k}^{n+1}-{\widehat{u}}_{i,k}^n\right)={\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i^{\circ }}\left(F\left(q_h\right)-F\left(u_{i,h}^{*}\right)\right){\partial }_x{\Phi }_{k}dxdt\hfill \\ & +{\int }_{t^n}^{t^{n+1}}{\Phi }_k\left(x_{i-1/2}\right)\left({\tilde{F}}_{i-1/2}\left(t\right)-F_{i-1/2}^{*}\right)+{\Phi }_k\left(x_{i+1/2}\right)\left(F_{i+1/2}^{*}-{\tilde{F}}_{i+1/2}\left(t\right)\right)dt\hfill \\ \hfill & +{\int }_{t^n}^{t^{n+1}}{\int }_{{\Omega }_i}{\Phi }_k\left(\left(S\left(q_h\right)-S\left(u_{i,h}^{*}\right)\right){\widehat{H}}_{i,l}{\partial }_x{\Phi }_l+G\left(q_h\right)-G\left(u_{i,h}^{*}\right)\right)dxdt,\hfill \end{array}$$

  (26)

  where

  $${\tilde{F}}_{i+1/2}\left(t\right)=\mathbb{F}\left(u_i^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{-}\left(t\right),u_{i+1}^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{+}\left(t\right)\right),$$

  (27)

  is the usual first order consistent numerical flux evaluated at the reconstructed predictor states

  $$u_i^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{-}\left(t\right),u_{i+1}^{*}\left(x_{i+1/2}\right)+v_{h,i+1/2}^{+}\left(t\right),$$

  where $u_i^{*}$ and $u_{i+1}^{*}$ are the stationary solutions computed at the step 1 at cells ${\Omega }_i$ and ${\Omega }_{i+1}$, respectively, and $v_{h,i+1/2}^{-}\left(t\right)$ and $v_{h,i+1/2}^{+}\left(t\right)$ are the extrapolated values of the fluctuation computed at step 2 onto the left and right cell interfaces, respectively.
  Additionally, $F_{i+1/2}^{*}$ denotes the evaluation of the physical flux at the local stationary solution at the inter-cells:

  $$F_{i+1/2}^{*}=F\left(u_i^{*}\left(x_{i+1/2}\right)\right).$$

Note that the added terms in the DG system (18) and ADER-DG system (26) are zero up to the order of the method. In this way, when the solution is far from a stationary solution, the solution yielded by these systems will be the same as their non-well-balanced counterparts (4) and (10) up to the order of the method. Additionally, in Section 5 a numerical example is given where a comparison between a well-balanced and non well-balanced simulation for a problem far from the equilibrium.

Finally, we stress that, although the techniques here have been described for the particular case of one-dimensional problems, all methods can be easily extended in a dimension by dimension fashion for the case of multiple spatial dimensions. For an implementation of the DG method in the higher spatial dimension case, we refer the reader to the work in [76].

4. A Well-Balanced WENO Limiting Procedure

We now detail the limiting procedure employed in this work, based on Zhong and Shu’s work in [83], where the authors propose a weighted essentially non-oscillatory (WENO) limiter for Runge–Kutta DG methods.

The main difference with our proposal is that we limit the fluctuations with respect to the stationary solutions, denoted by $v_h,$ rather than the solution itself. Once the fluctuation has been limited, the numerical solution can be recovered by adding the stationary solution.

The limiter can be seen as a two-steps procedure as in [83]: in the first stage, the cells that need limiting, or troubled cells, are identified. Then, in the second step, the cells that were deemed as troubled are limited.

For both Runge–Kutta and ADER time-marching discretisation, the limiting procedure is applied a priori and only once at the beginning of each time step.

4.1. Detector of Troubled Cells

We use the TVD minmod limiter detector following ideas detailed in [41,83]. We first define

$${\overline{v}}_i=\frac{1}{\Delta x}{\int }_{{\Omega }_i}{v}_h(x,t^n)dx,$$

(28)

where $v_h(x,t^n)$ denotes the fluctuation with respect to the discrete stationary solution $u_{i,h}^{*}\left(x\right)$ that has been computed as it is detailed in Section 3.1,

$${v_h}_{{|}_{{\Omega }_i}}(x,t^n)={u_h}_{{|}_{{\Omega }_i}}(x,t^n)-u_{i,h}^{*}\left(x\right),x\in {\Omega }_i.$$

(29)

We also define

$${\widehat{v}}_R=v_{i+\frac{1}{2}}^{-}-{\overline{v}}_i,{\widehat{v}}_L={\overline{v}}_i-{\tilde{v}}_{i-\frac{1}{2}}^{+},$$

(30)

where $v_{i-1/2}^{+}$ and $v_{i+1/2}^{-}$ are the evaluation of the polynomial fluctuation $v_h(x,t^n)$ onto the left and right cell interfaces respectively. Then, the minmod is applied as follows:

$${\sigma }_R=\mathrm{minmod}\left({\widehat{v}}_R,{\Delta }^{+}{v}_i,{\Delta }^{-}{v}_i\right),{\sigma }_L=\mathrm{minmod}\left({\widehat{v}}_L,{\Delta }^{+}{v}_i,{\Delta }^{-}{v}_i\right),$$

(31)

where

$${\Delta }^{+}{v}_i={\overline{v}}_{i+1}-{\overline{v}}_i,{\Delta }^{-}{v}_i={\overline{v}}_i-{\overline{v}}_{i-1},$$

(32)

and the following minmod functions are given by

$$\mathrm{minmod}(a,b,c)=\left\{\begin{array}{cc}a\hfill & \mathrm{if}\left|a\right|\le M\Delta x,\hfill \\ \mathrm{m}(a,b,c)\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(33)

$$\mathrm{m}(a,b,c)=\left\{\begin{array}{cc}smin(a,b,c)\hfill & \mathrm{if}s=\mathrm{sign}\left(a\right)=\mathrm{sign}\left(b\right)=\mathrm{sign}\left(c\right),\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(34)

Here, M is a parameter to be properly chosen depending on the PDE system. Finally, a cell is marked as troubled if $\mathrm{minmod}(a,b,c)\ne a.$

4.2. Limiting Procedure

The limiting technique for the cells marked as troubled by the detector described in the previous Subsection is now detailed. The general idea is to reconstruct a new polynomial solution at a troubled cell as a convex combination of its neighbouring polynomials.

If the cell ${\Omega }_i$ is deemed troubled, then the first step consist on projecting the fluctuation at each neighbour element ${\Omega }_j,j\in \{i-1,i,i+1\},$ into a first-order DG polynomial with the following form:

$$v_{p,j}\left(x\right)={\displaystyle \frac{1}{\Delta x}}{\int }_{{\Omega }_j}{v}_{j,h}\left(x\right)dx+{\displaystyle \frac{(x-x_{j,c})}{\Delta x}}{\int }_{{\Omega }_j}\frac{\partial }{{\partial }_x}{v}_{j,h}\left(x\right)dx,x\in {\Omega }_j,$$

(35)

that is a first order polynomial that preserves the cell average of $v_{j,h}\left(x\right)$. Here, the sub-index p denotes the projected polynomial using (35). We remark that we have used again $v_h$ to denote the fluctuation with respect to the stationary solution as in (29). Likewise, $x_{j,c}$ denotes the centre of the element ${\Omega }_j$ and the time dependency has been dropped to simplify the notation.

This projected polynomial is then used to rewrite the new polynomial fluctuation on the element ${\Omega }_i$ as follows:

$${v_h}_{{|}_{{\Omega }_i}}^{\left(\mathrm{new}\right)}\left(x\right)={\omega }_{0,i}{\overline{v}}_{i-1,p}\left(x\right)+{\omega }_{1,i}{v}_{i,p}\left(x\right)+{\omega }_{2,i}{\overline{v}}_{i+1,p}\left(x\right),x\in {\Omega }_i,$$

(36)

where

$${\overline{v}}_{i\pm 1,p}\left(x\right)=v_{i\pm 1,p}\left(x\right)-\frac{1}{\Delta x}{\int }_{{\Omega }_i}{v}_{i\pm 1,p}\left(x\right)dx+\frac{1}{\Delta x}{\int }_{{\Omega }_i}{v}_{i,h}\left(x\right)dx,$$

(37)

which ensures that ${v_h}_{{|}_{{\Omega }_i}}^{\left(\mathrm{new}\right)}\left(x\right)$ has the same cell average as ${v_h}_{{|}_{{\Omega }_i}}\left(x\right)$. As we have made use of the fluctuation to calculate the new limited solution, the final limited solution is computed as

$${u_h}_{{|}_{{\Omega }_i}}^{\left(\mathrm{new}\right)}\left(x\right)=u_{i,h}^{*}\left(x\right)+{v_h}_{{|}_{{\Omega }_i}}^{\left(\mathrm{new}\right)}\left(x\right).$$

(38)

The normalized weights for the WENO convex combination in (36) follows the classical WENO procedure:

$${\omega }_{l,i}={\displaystyle \frac{{\overline{\omega }}_{l,i}}{{\overline{\omega }}_{0,i}+{\overline{\omega }}_{1,i}+{\overline{\omega }}_{2,i}}},{\overline{\omega }}_{l,i}=\frac{{\gamma }_l}{{(ϵ+{\beta }_l)}^r},l=0,1,2.$$

(39)

The parameters ${\gamma }_l$ used in this paper are ${\gamma }_1=0.998,$ and ${\gamma }_{0,2}=0.001;$ $ϵ$ is set to $ϵ={10}^{-6}.$ The employed values for r are in $\{0.25,2\},$ and will be better detailed in the next Subsection. Finally, the smooth indicators ${\beta }_{l,i}$ are defined as

$${\beta }_{l,i}=\sum _{s=1}^2{\int }_{{\Omega }_i}\Delta x^{2s-1}{\left(\frac{{\partial }^s}{\partial x^s}{v}_{p,l}\left(x\right)\right)}^{2}dx.$$

(40)

Again, we stress that the general limiting procedure described here is applied once per time-step for either Runge–Kutta or ADER time discretisation.

4.3. Double-Loop Limiting Procedure

Even though the component-wise WENO technique is essentially non-oscillatory, the resulting scheme may still generate spurious oscillations, especially at very high order. We, therefore, propose the following strategy.

• First, the procedure (35)–(38) with the value $r=2$ is used to produce a candidate solution $u_h^c\left(x\right).$ This candidate solution is evaluated once more by the detector described in Section 4.1.
• Then, if the element is once again considered as troubled, the limiting strategy (35)–(38) will be repeated changing the projector operator. Now, the fluctuation is projected into a constant value given by

  $$v_{j,p}\left(x\right)=\frac{1}{\Delta x}{\int }_{{\Omega }_j}{v}_{j,h}\left(x\right)dx.$$

  (41)
  Additionally, the parameter r will be set to $r=0.25$. This strategy will ensure a smoother discrete solution for highly oscillating regions near strong discontinuities.

That finishes the limiting technique description.

Remark 4.

It is essential to notice that the stationary solutions considered in this work are smooth, and therefore the detector does not mark cells as troubled. However, let us suppose that roundoff errors or local extremes provoke that an element is marked as troubled in the presence of a stationary solution. In that case, the use of the fluctuation (29) will ensure that the stationary solution is preserved, since the procedure (35)–(38) is exact for the null operator.

5. Numerical Test

In this Section, some numerical experiments are conducted to test the proposed numerical first, second, third and fourth-order schemes for Runge–Kutta and ADER time discretisations.

We consider four systems: Burgers, Shallow water, Euler with gravity and the hyperbolic dispersive Gavrilyuk–Favrie equations. For each system, a brief description is given. Then, we test the ability of the numerical scheme to preserve stationary solutions. Afterwards, some numerical tests consisting of a perturbation of a stationary solution are shown alongside the behaviour of the well-balanced limiter technique proposed in this paper.

The CFL number in (11) is always set to $0.9$, meanwhile the limiter value M defined in (33) is set to $1,$ except for the problem (63) that is set to $0.1.$ For the rest of the limiter parameters appearing in Section 4, we maintain the same values that were defined in that Subsection.

Remark that for the two-spatial domain dispersive system, the CFL condition (11) for DG methods in the case of two dimensional domains reads

$$\Delta t\le {\displaystyle \frac{1}{2}}\frac{CFL}{2N+1}\frac{min\{\Delta x,\Delta y\}}{|{\lambda }_{max}|},CFL\in (0,1),$$

and we set $CFL=0.9.$ Besides, for the two dimensional system, the values of the $\gamma$ weights for the limiter are set to ${\gamma }_c=0.996$ for the central element, and $\gamma =0.001$ for the four neighbouring elements.

5.1. Burgers’ Equation

Let us consider the one-dimensional scalar balance law consisting of Burgers’ equation with a source term written in the form (1) with

$$U=u,F\left(U\right)={\displaystyle \frac{u^2}{2}},S\left(U\right)=u^2,G\left(U\right)=0,$$

(42)

whose stationary solutions are given by

$$u\left(x\right)=C{e}^{H\left(x\right)},C\in \mathbb{R},$$

(43)

and with $H\left(x\right)$ a known continuous function.

5.1.1. Preservation of a Stationary Solution

In this numerical experiment we set as initial condition a stationary solution of the form (43) to ensure the well-balancing properties of the proposed methods for the system (42). The initial condition considered is given by

$$u=e^{H\left(x\right)},H=x.$$

(44)

The computational domain is $\Omega =[-1,1].$ Dirichlet and open boundary conditions were imposed at the left and at the right of the domain, respectively. Figure 1 shows the error between the computed numerical solutions and the stationary one given by $e^{H\left(x\right)}$. As it can be seen in Figure 1 and

Table 1. Numerical validation of the well-balanced and non well-balanced methods for the stationary problem (44) (Burgers’ equation). Table contains $L_1$ errors for the variable u at a final time $t=10$ s for both Runge–Kutta and ADER DG schemes of order $N+1=1,2,3,4.$ Uniform Cartesian meshes of $N_x$ elements has been used.

		Non Well-Balanced		Well-Balanced
N	${\mathit{N}}_{\mathit{x}}$	Runge-Kutta	ADER	Runge-Kutta	ADER
0	25	$1.39\times {10}^{-01}$	$1.39\times {10}^{-01}$	0	0
	50	$6.52\times {10}^{-02}$	$6.52\times {10}^{-02}$	$2.66\times {10}^{-17}$	0
	100	$3.17\times {10}^{-02}$	$3.17\times {10}^{-02}$	$1.55\times {10}^{-17}$	0
	200	$1.56\times {10}^{-02}$	$1.56\times {10}^{-02}$	0	0
1	25	$6.47\times {10}^{-04}$	$5.37\times {10}^{-04}$	$8.88\times {10}^{-18}$	0
	50	$1.60\times {10}^{-04}$	$1.34\times {10}^{-04}$	0	0
	100	$3.96\times {10}^{-05}$	$3.33\times {10}^{-05}$	$2.22\times {10}^{-18}$	$1.11\times {10}^{-18}$
	200	$9.85\times {10}^{-06}$	$8.29\times {10}^{-06}$	$2.22\times {10}^{-18}$	$5.55\times {10}^{-19}$
2	25	$4.49\times {10}^{-06}$	$3.80\times {10}^{-06}$	0	$2.47\times {10}^{-18}$
	50	$5.91\times {10}^{-07}$	$5.25\times {10}^{-07}$	$7.40\times {10}^{-18}$	0
	100	$7.57\times {10}^{-08}$	$6.85\times {10}^{-08}$	$2.47\times {10}^{-18}$	0
	200	$9.57\times {10}^{-09}$	$8.75\times {10}^{-09}$	$1.23\times {10}^{-18}$	$9.25\times {10}^{-19}$
3	25	$2.56\times {10}^{-08}$	$1.10\times {10}^{-08}$	$4.63\times {10}^{-18}$	$3.09\times {10}^{-18}$
	50	$1.70\times {10}^{-09}$	$7.85\times {10}^{-10}$	$1.54\times {10}^{-18}$	$7.53\times {10}^{-18}$
	100	$1.09\times {10}^{-10}$	$5.61\times {10}^{-11}$	$2.32\times {10}^{-18}$	$4.15\times {10}^{-18}$
	200	$6.91\times {10}^{-12}$	$4.00\times {10}^{-12}$	$1.54\times {10}^{-18}$	$2.27\times {10}^{-18}$

, the well-balanced methods can preserve the stationary solution at a final time $t=10$ s.

5.1.2. Perturbed Stationary Solutions

We consider now a perturbation of the stationary solution (44) given by

$$u=e^H+0.3e^{-200{(x+0.5)}^2},H=x.$$

(45)

The computational domain $\Omega =[-1,1]$ is covered with a coarse mesh of $N_x=100$ elements. Dirichlet and open boundary conditions are imposed at the left and the right of the domain, respectively. The numerical experiment is run up to time $t=10$ s with the fourth-order DG method and for both Runge–Kutta and ADER time discretisations.

Figure 2 shows the error between the computed numerical solutions and the stationary one given by $e^H$. As it can be seen, once the perturbation has left the domain, the well-balanced numerical scheme is able to recover the stationary solution $e^H$, in contrast with the non well-balanced scheme.

5.1.3. Limited Simulation

In this final example, we aim to check the robustness of the proposed limiter technique described in this paper. To do that, we consider the initial condition

$$u=\frac{1}{2}+sin\left(x\right),H=0.$$

(46)

The computational domain $\Omega =[0,2\pi ]$ is covered with a coarse mesh of $N_x=80$ elements. Periodic boundary conditions were imposed everywhere. Figure 3 shows the numerical approximation of u with different time-marching discretisations at the final time $t=1.5$ s, as well as a reference solution that has been computed by a first order scheme with $N_x=$ 20,000 elements. As it can be seen, the proposed schemes provide an accurate approximation without producing any spurious oscillations, even for the case of coarse meshes.

5.1.4. Well-Balanced and Non-Well-Balanced Schemes Comparison

We perform now a comparison between the well-balanced and non-well-balanced schemes for a solution which is not a perturbation of an equilibrium state. Particularly, the following initial condition is considered:

$$u=0.1e^{-200{(x+0.5)}^2},H=x.$$

(47)

Periodic boundary conditions and 100 discretization points are considered in a computational domain $\Omega =[-1,1]$. Figure 4 depicts the results for a Runge–Kutta second and third order scheme. The reference solution has been computed using a first order method with 20,000 discretization points. As it can be seen, the results shows that the well-balanced and non well-balanced methods are in fact very similar regardless of the order of the method. Similar results are obtained for ADER time marching methods.

5.2. Shallow Water Equations

Let us consider the one-dimensional shallow water equations written in the form (1) with

$$U=\left(\begin{array}{c}h\\ hu\end{array}\right),F\left(U\right)=\left(\begin{array}{c}hu\\ h{u}^2+{\displaystyle \frac{1}{2}}g{h}^2\end{array}\right),S\left(U\right)=\left(\begin{array}{c}0\\ gh\end{array}\right),G\left(U\right)=0,$$

(48)

where $h=h(x,t)$ is the total water depth, u is the depth-averaged velocity on the $x$-direction, and $H=H\left(x\right)$ is the known still water depth that we will suppose to be continuous. Finally, $g=9.81$ is the gravitational acceleration. In this case we adapt the well-balanced strategy described in Section 3 to preserve the stationary solutions corresponding to water at rest:

$$u=0,\eta =cst,$$

(49)

where $\eta =h-H$ is the free-surface elevation.

A well-balance method could also be design following the ideas described in [93].

5.2.1. Order of Accuracy Test

We begin by checking numerically the order of accuracy of the numerical scheme for a Runge–Kutta and ADER time discretisation for a second-, third-, and fourth-order schemes. As pointed before, we consider here only a well-balanced method that preserves the water at rest solution adapting the strategy described in Section 3. The computational domain is $\Omega =[-5,5]$, and periodic boundary conditions are imposed everywhere. We then simulate the propagation of a Gaussian wave, with the initial condition given by

$$u=0,h=H+0.1e^{-x^2},H=1-0.2e^{-x^2},$$

(50)

and until the final time $t=1$ second, with a set of refined meshes.

Table 2. Numerical convergence results for the initial condition (50) (shallow water equations). High-order Runge–Kutta and ADER DG schemes of order $N+1=2,3,4.$ Uniform Cartesian meshes of $N_x$ elements has been used. The $L_1$ errors refer to the variables h and $hu$ at a final time $t=1$ s.

		Runge-Kutta				ADER
		$\mathit{h}$		$\mathit{hu}$		$\mathit{h}$		$\mathit{hu}$
$\mathit{N}$	${\mathit{N}}_{\mathit{x}}$	Error	Order	Error	Order	Error	Order	Error	Order
1	25	$3.93\times {10}^{-03}$	-	$1.27\times {10}^{-02}$		$3.60\times {10}^{-03}$	-	$1.17\times {10}^{-02}$	-
	50	$5.92\times {10}^{-04}$	$2.73$	$1.86\times {10}^{-03}$	$2.77$	$8.91\times {10}^{-04}$	$2.01$	$2.79\times {10}^{-03}$	$2.07$
	100	$1.11\times {10}^{-04}$	$2.41$	$3.49\times {10}^{-04}$	$2.41$	$2.21\times {10}^{-04}$	$2.01$	$6.78\times {10}^{-04}$	$2.04$
	200	$2.16\times {10}^{-05}$	$2.37$	$6.81\times {10}^{-05}$	$2.36$	$5.50\times {10}^{-05}$	$2.00$	$1.67\times {10}^{-04}$	$2.02$
2	25	$7.37\times {10}^{-05}$	-	$2.19\times {10}^{-04}$	-	$2.97\times {10}^{-04}$	-	$9.46\times {10}^{-04}$	-
	50	$3.27\times {10}^{-06}$	$4.49$	$9.72\times {10}^{-06}$	$4.50$	$3.28\times {10}^{-05}$	$3.18$	$1.06\times {10}^{-04}$	$3.15$
	100	$2.47\times {10}^{-07}$	$3.73$	$7.30\times {10}^{-07}$	$3.74$	$3.79\times {10}^{-06}$	$3.11$	$1.26\times {10}^{-05}$	$3.08$
	200	$2.23\times {10}^{-08}$	$3.47$	$6.92\times {10}^{-08}$	$3.40$	$4.60\times {10}^{-07}$	$3.04$	$1.55\times {10}^{-06}$	$3.02$
3	25	$3.52\times {10}^{-06}$	-	$1.02\times {10}^{-05}$	-	$9.48\times {10}^{-06}$	-	$3.43\times {10}^{-05}$	-
	50	$1.41\times {10}^{-07}$	$4.64$	$3.52\times {10}^{-07}$	$4.86$	$6.60\times {10}^{-07}$	$3.84$	$2.24\times {10}^{-06}$	3.94
	100	$8.95\times {10}^{-09}$	$3.98$	$2.21\times {10}^{-08}$	$3.99$	$4.41\times {10}^{-08}$	$3.91$	$1.44\times {10}^{-07}$	3.96
	200	$5.63\times {10}^{-10}$	$3.99$	$1.39\times {10}^{-09}$	$3.99$	$2.75\times {10}^{-09}$	$4.00$	$9.12\times {10}^{-09}$	3.98

shows the $L_1$ errors and the numerical convergence rates. The expected convergence order of the Runge–Kutta and ADER DG methods can be observed.

5.2.2. Preserving Stationary Solutions

In this experiment, we set as initial condition a stationary solution of the form (49) to ensure that the well-balanced methods are able to preserve it up to machine precision, while the non well-balanced methods fail. The initial condition considered is given by

$$u=0,h=H,H=1-0.8e^{-100x^2}.$$

(51)

The computational domain is $\Omega =[-1,1]$ and Dirichlet boundary conditions are imposed everywhere. As it can be seen in Figure 5 and

Table 3. Numerical validation of the well-balanced and non well-balanced methods for the stationary problem (51) (shallow water equations). Table contains $L_1$ errors for the water depth variable h at a final time $t=10$ s for Runge–Kutta and ADER DG schemes of order $N+1=1,2,3,4.$ Uniform Cartesian meshes of $N_x$ elements has been used.

		Non-Well-Balanced		Well-Balanced
N	${\mathit{N}}_{\mathit{x}}$	Runge–Kutta	ADER	Runge-Kutta	ADER
0	25	$1.92\times {10}^{-02}$	$1.92\times {10}^{-02}$	$3.00\times {10}^{-17}$	$3.00\times {10}^{-17}$
	50	$6.08\times {10}^{-03}$	$6.08\times {10}^{-03}$	$1.27\times {10}^{-16}$	$3.34\times {10}^{-16}$
	100	$1.40\times {10}^{-03}$	$1.40\times {10}^{-03}$	$4.86\times {10}^{-17}$	$1.64\times {10}^{-17}$
	200	$5.37\times {10}^{-04}$	$5.32\times {10}^{-04}$	$5.97\times {10}^{-18}$	$2.30\times {10}^{-17}$
1	25	$1.16\times {10}^{-03}$	$2.62\times {10}^{-03}$	$6.66\times {10}^{-17}$	$6.88\times {10}^{-16}$
	50	$2.73\times {10}^{-04}$	$3.47\times {10}^{-04}$	$2.06\times {10}^{-16}$	$8.61\times {10}^{-16}$
	100	$3.76\times {10}^{-05}$	$5.48\times {10}^{-05}$	$2.82\times {10}^{-16}$	$1.42\times {10}^{-15}$
	200	$4.81\times {10}^{-06}$	$8.48\times {10}^{-06}$	$4.37\times {10}^{-16}$	$1.61\times {10}^{-15}$
2	25	$3.22\times {10}^{-04}$	$4.09\times {10}^{-04}$	$1.52\times {10}^{-16}$	$1.36\times {10}^{-16}$
	50	$3.20\times {10}^{-05}$	$4.39\times {10}^{-05}$	$1.85\times {10}^{-16}$	$1.20\times {10}^{-15}$
	100	$2.36\times {10}^{-06}$	$4.47\times {10}^{-06}$	$2.74\times {10}^{-16}$	$1.58\times {10}^{-15}$
	200	$1.51\times {10}^{-07}$	$4.95\times {10}^{-07}$	$3.67\times {10}^{-16}$	$1.96\times {10}^{-15}$
3	25	$1.84\times {10}^{-05}$	$4.18\times {10}^{-05}$	$1.47\times {10}^{-16}$	$1.63\times {10}^{-15}$
	50	$1.64\times {10}^{-06}$	$2.42\times {10}^{-06}$	$3.11\times {10}^{-16}$	$2.50\times {10}^{-15}$
	100	$6.43\times {10}^{-08}$	$6.25\times {10}^{-08}$	$5.78\times {10}^{-16}$	$4.14\times {10}^{-16}$
	200	$2.08\times {10}^{-09}$	$3.57\times {10}^{-09}$	$6.75\times {10}^{-16}$	$6.98\times {10}^{-16}$

and

Table 4. Numerical validation of the well-balanced and non well-balanced methods for the stationary problem (51) (shallow water equations). Table contains $L_1$ errors for the conserved variable $hu$ at a final time $t=10$ s for Runge–Kutta and ADER DG schemes of order $N+1=1,2,3,4.$ Uniform Cartesian meshes of $N_x$ elements has been used.

		Non-Well-Balanced		Well-Balanced
N	${\mathit{N}}_{\mathit{x}}$	Runge–Kutta	ADER	Runge–Kutta	ADER
0	25	$6.12\times {10}^{-02}$	$6.12\times {10}^{-02}$	$4.40\times {10}^{-15}$	$4.39\times {10}^{-15}$
	50	$3.96\times {10}^{-02}$	$3.96\times {10}^{-02}$	$1.06\times {10}^{-15}$	$1.11\times {10}^{-15}$
	100	$2.20\times {10}^{-02}$	$2.21\times {10}^{-02}$	$4.09\times {10}^{-16}$	$3.52\times {10}^{-16}$
	200	$1.20\times {10}^{-02}$	$1.20\times {10}^{-02}$	$1.69\times {10}^{-16}$	$1.79\times {10}^{-16}$
1	25	$3.56\times {10}^{-03}$	$6.49\times {10}^{-03}$	$5.68\times {10}^{-16}$	$1.07\times {10}^{-14}$
	50	$3.35\times {10}^{-04}$	$7.43\times {10}^{-04}$	$5.33\times {10}^{-16}$	$2.19\times {10}^{-15}$
	100	$2.28\times {10}^{-05}$	$1.38\times {10}^{-04}$	$8.28\times {10}^{-16}$	$9.95\times {10}^{-15}$
	200	$1.44\times {10}^{-06}$	$3.38\times {10}^{-05}$	$1.30\times {10}^{-15}$	$3.89\times {10}^{-15}$
2	25	$5.27\times {10}^{-04}$	$7.48\times {10}^{-04}$	$3.97\times {10}^{-16}$	$4.50\times {10}^{-16}$
	50	$1.40\times {10}^{-04}$	$1.82\times {10}^{-04}$	$6.83\times {10}^{-16}$	$4.28\times {10}^{-15}$
	100	$1.96\times {10}^{-05}$	$2.69\times {10}^{-05}$	$7.47\times {10}^{-16}$	$4.23\times {10}^{-15}$
	200	$2.52\times {10}^{-06}$	$3.45\times {10}^{-06}$	$9.91\times {10}^{-16}$	$2.79\times {10}^{-15}$
3	25	$1.23\times {10}^{-04}$	$3.02\times {10}^{-04}$	$3.62\times {10}^{-16}$	$1.51\times {10}^{-14}$
	50	$2.59\times {10}^{-06}$	$1.02\times {10}^{-05}$	$7.01\times {10}^{-16}$	$4.29\times {10}^{-15}$
	100	$4.85\times {10}^{-08}$	$3.25\times {10}^{-07}$	$9.93\times {10}^{-16}$	$1.08\times {10}^{-15}$
	200	$7.75\times {10}^{-10}$	$2.20\times {10}^{-08}$	$1.87\times {10}^{-15}$	$1.55\times {10}^{-15}$

, the well-balanced methods are able to preserve the stationary solution at a final time $t=10$ s for both time-marching discretisations and relative coarse meshes. Again, the numerical solutions are compared with the stationary solution (51) as reference.

5.2.3. Perturbed Stationary Solutions

We consider now a perturbation of the stationary solution (51) given by

$$u=0,h=0.1e^{-100x^2}+H,H=1-0.8e^{-100x^2}.$$

(52)

The computational domain $\Omega =[-1,1]$ is covered with $N_x=100$ elements. Dirichlet and free-outflow boundary conditions are imposed at the right and the left of the domain, respectively.

As it can be seen in Figure 6, the well-balanced scheme can recover the stationary solution at a final time $t=10$ s with the fourth-order DG method and for both Runge–Kutta and ADER time discretisations. In this case, we plot the difference between the numerically computed free surfaces and the exact stationary (51), which is zero.

Moreover, in order to show the better performance of a well-balanced method against a non well-balanced method, we propose a similar numerical test on a very coarse mesh composed by 20 elements. Here, we consider the initial condition given by

$$u=0,h={10}^{-4}{e}^{-100x^2}+H,H=1-0.8e^{-100x^2}.$$

(53)

The numerical computation is then run up to the short time $t=0.5$ s; the perturbation does not reach boundaries in this case, and therefore periodic boundary conditions were imposed everywhere.

As the perturbation is of the order ${(\Delta x)}^{N+1}$, $N+1=4$ being the order of the scheme, we can expect that only the well-balanced method yields satisfactory results. In contrast, the non well-balanced scheme introduces oscillations of that order that destroy the solution. Here we only show the numerical results in Figure 7 for the Runge–Kutta time discretisation. Similar results can be observed for ADER time discretisations.

5.2.4. Limited Simulation

In order to check the robustness of the limiting procedure described in Section 4.3, we consider now a situation far from the equilibrium, where the use of some limiting technique is mandatory. The computational domain $\Omega =[-5,5]$ is covered now with a very refined mesh of $N_x=400$ elements. The initial condition is set to

$$u=0,h=0.1e^{-5x^2}+H,H=1-0.8e^{-x^2}.$$

(54)

Periodic boundary conditions are imposed to obtain more interactions between waves. Some snapshots of the computed free surface are depicted in Figure 8. As it can be observed, there are no oscillations at the neighbourhood of the shock waves, which indicates that the limiting approach provides a high overall quality of the computed numerical solution. The solution is depicted alongside a reference solution that has been computed by a first order scheme with $N_x=$ 20,000 elements.

Note that the small oscillations appearing with the ADER method could be prevented with an ad hoc tuning of the parameters associated with the limiter described in Section 4. Here, we have maintained the same parameters for both approaches to compare both Runge–Kutta and ADER methods.

5.3. Compressible Euler Equations with Gravitational Force

We consider the gas dynamic Euler equations written as in (1):

$$U=\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right),F\left(U\right)=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ u\left(E+p\right)\end{array}\right),S\left(U\right)=\left(\begin{array}{c}0\\ -\rho \\ -\rho u\end{array}\right),G\left(U\right)=0,$$

(55)

$\rho \ge 0$ being the density, u the velocity, E the total energy, and $H\left(x\right)$ a given continuous gravitational potential. The internal energy e is given by $\rho e=E-\frac{1}{2}\rho u^2$. The pressure $p\ge 0$ is given by e through the equation of state (ideal gas),

$$p=(\gamma -1)\rho e,$$

(56)

where $\gamma >1$ is the adiabatic constant that will take the value $\frac{5}{3}$ unless stated otherwise.

In this work, we will focus on the hydrostatic steady states solutions given by

$$u=0,{\partial }_{x}p=-\rho H^{\prime }\left(x\right).$$

(57)

According to [35], a family of stationary solutions of (55) with $u=0$ are

$$u=0,\rho =\alpha (H,C_1)\ge 0,p\left(x\right)=\beta (H,C_1,C_2)\ge 0,E=\frac{p\left(x\right)}{\gamma -1},$$

(58)

where $\alpha$ is a given continuous function and $\beta$ is given by

$${\partial }_{\sigma }\beta (H,C_1,C_2)=-\alpha (H,C1).$$

The numerical methods are then chosen to preserve one of those families given by the choice of $\alpha :$

$$\alpha (H,C_1)=C_1{e}^{-H},$$

(59)

that results in the following set of stationary solutions with $u=0$,

$$u=0,\rho =C_1{e}^{-H}\ge 0,p=C_1{e}^{-H}+C_2\ge 0,E=\frac{p}{\gamma -1}.$$

(60)

5.3.1. Preserving Stationary Solutions

In this experiment, we set as initial condition a stationary solution of the form (60) to check the well-balancing properties of the proposed methods for the system (55). The initial condition is given by

$$u=0,\rho =e^{-H},p(x,0)=e^{-H},E={\displaystyle \frac{p}{\gamma -1}},H=x.$$

(61)

The computational domain is $\Omega =[-1,1]$ and Dirichlet boundary conditions were imposed everywhere. Figure 9 shows the error between the compute numerical solutions and the stationary one given by (61). As it can be seen in Figure 9 and

Table 5. Numerical validation of the well-balanced and non well-balanced methods for the stationary problem (61) (Euler equations). Table contains $L_1$ errors for the density variable $\rho$ at a final time $t=10$ s for both Runge–Kutta and ADER DG schemes of order $N+1=1,2,3,4.$ Uniform Castesian meshes of $N_x$ elements has been used.

		Non-Well-Balanced		Well-Balanced
N	${\mathit{N}}_{\mathit{x}}$	Runge–Kutta	ADER	Runge–Kutta	ADER
0	25	$5.13\times {10}^{-02}$	$1.03\times {10}^{-02}$	$7.99\times {10}^{-17}$	$1.31\times {10}^{-16}$
	50	$3.18\times {10}^{-02}$	$5.46\times {10}^{-03}$	$2.31\times {10}^{-16}$	$4.42\times {10}^{-16}$
	100	$1.77\times {10}^{-02}$	$3.21\times {10}^{-03}$	$1.94\times {10}^{-16}$	$6.21\times {10}^{-16}$
	200	$9.41\times {10}^{-03}$	$1.85\times {10}^{-03}$	$3.26\times {10}^{-16}$	$6.61\times {10}^{-16}$
1	25	$2.35\times {10}^{-04}$	$3.08\times {10}^{-04}$	$7.46\times {10}^{-16}$	$1.15\times {10}^{-15}$
	50	$7.10\times {10}^{-05}$	$8.81\times {10}^{-05}$	$1.20\times {10}^{-15}$	$1.47\times {10}^{-15}$
	100	$1.94\times {10}^{-05}$	$2.34\times {10}^{-05}$	$2.15\times {10}^{-15}$	$2.65\times {10}^{-15}$
	200	$9.26\times {10}^{-06}$	$6.03\times {10}^{-06}$	$3.17\times {10}^{-15}$	$4.21\times {10}^{-15}$
2	25	$3.52\times {10}^{-06}$	$2.98\times {10}^{-06}$	$1.06\times {10}^{-15}$	$2.84\times {10}^{-15}$
	50	$5.40\times {10}^{-07}$	$4.67\times {10}^{-07}$	$1.78\times {10}^{-15}$	$5.01\times {10}^{-15}$
	100	$7.44\times {10}^{-08}$	$6.58\times {10}^{-08}$	$4.40\times {10}^{-15}$	$9.39\times {10}^{-15}$
	200	$9.75\times {10}^{-09}$	$8.73\times {10}^{-09}$	$4.93\times {10}^{-15}$	$1.36\times {10}^{-14}$
3	25	$8.20\times {10}^{-09}$	$1.81\times {10}^{-08}$	$1.26\times {10}^{-15}$	$1.35\times {10}^{-14}$
	50	$7.14\times {10}^{-10}$	$1.53\times {10}^{-09}$	$2.74\times {10}^{-15}$	$2.83\times {10}^{-14}$
	100	$5.18\times {10}^{-11}$	$1.09\times {10}^{-10}$	$7.71\times {10}^{-15}$	$4.53\times {10}^{-14}$
	200	$3.48\times {10}^{-12}$	$7.00\times {10}^{-12}$	$9.54\times {10}^{-15}$	$7.16\times {10}^{-14}$

, the well-balanced methods can preserve the stationary solution at a final time $t=10$ s.

5.3.2. Limited Simulation

We check again the well behaviour of the proposed mixed-order technique described in this paper. To do that, we set $H=0,\gamma =1.4,$ and consider the following discontinuous initial condition given by

$$\left\{\begin{array}{c}\left(\begin{array}{ccc}{\rho }_L,& u_L,& p_L\end{array}\right)\mathrm{if}x<0,\hfill \\ \\ \left(\begin{array}{ccc}{\rho }_L,& u_R,& p_L\end{array}\right)\mathrm{if}x\ge 0.\hfill \end{array}\right.$$

(62)

Then, we set the values for the two well-known numerical tests:

• The sod problem [94]:

  $$\left(\begin{array}{ccc}{\rho }_L,& u_L,& p_L\end{array}\right) = \left(\begin{array}{ccc}1,& 0,& 1\end{array}\right),\left(\begin{array}{ccc}{\rho }_R,& u_R,& p_R\end{array}\right) = \left(\begin{array}{ccc}0.125& 0& 0.1\end{array}\right).$$

  (63)
• The Lax problem [95]:

  $$\left(\begin{array}{ccc}{\rho }_L& u_L& p_L\end{array}\right) = \left(\begin{array}{ccc}0.445,& 0.698,& 3.528\end{array}\right),\left(\begin{array}{ccc}{\rho }_R,& u_R,& p_R\end{array}\right) = \left(\begin{array}{ccc}0.5,& 0,& 0.571\end{array}\right).$$

  (64)

In both cases, the computational domain $\Omega =[-5,5]$ is covered with $N_x=300$ elements. The limiter value M defined in (33) is set to $0.01,$ and the CFL number to $0.5.$ Free outflow boundary conditions were imposed everywhere.

Figure 10 and Figure 11 show the results with the fourth order Runge–Kutta and ADER DG methods for the Sod and Lax problems respectively. The numerical results are compared against a reference solution that has been computed using a first order method with $N_x=$ 20,000 elements.

5.4. Dispersive Gavrilyuk–Favrie Equations

Finally, we will consider a hyperbolic dispersive shallow water system that was firstly introduced in [96] for the case of a flat bottom bathymetry, and later extended for the non-flat bottom case in [78]. The system can be seen as a hyperbolic relaxation approach of the elliptic Serre–Green–Naghdi equations, and is formally deduced from a master Lagrangian formalism. The system for the case of two-spatial domains reads

$${\partial }_{t}U+{\partial }_x{F}_1\left(U\right)+{\partial }_y{F}_2\left(U\right)=S_1\left(U\right){\partial }_{x}H+S_2\left(U\right){\partial }_{y}H+G\left(U\right),$$

(65)

where

$$\begin{array}{cc}\hfill & U=\left(\begin{array}{c}h\\ hu\\ hv\\ hw\\ h\xi \end{array}\right),F_1=\left(\begin{array}{c}hu\\ h{u}^2+{\displaystyle \frac{1}{2}}g{h}^2+hp\\ uhv\\ uhw\\ uh\xi \end{array}\right),F_2=\left(\begin{array}{c}hv\\ vhu\\ h{v}^2+{\displaystyle \frac{1}{2}}g{h}^2+hp\\ vhw\\ vh\xi \end{array}\right)\hfill \\ \hfill & S_1=\left(\begin{array}{c}0\\ gh+{\displaystyle \frac{3h}{2\xi }}p\\ 0\\ 0\\ {\displaystyle \frac{3}{2}}hu\end{array}\right),S_2=\left(\begin{array}{c}0\\ 0\\ gh+{\displaystyle \frac{3h}{2\xi }}p\\ 0\\ {\displaystyle \frac{3}{2}}hv\end{array}\right),G=\left(\begin{array}{c}0\\ 0\\ 0\\ \lambda \left(1-{\displaystyle \frac{\xi }{h}}\right)\\ hw\end{array}\right).\hfill \end{array}$$

(66)

Here, $h=h(x,y,t)$ is the water depth, $H=H(x,y)$ is the known still water depth that we will assume to be continuous, and $u,v,w$ are the depth averaged velocities at the $x,y$ and z direction respectively. Furthermore, $p(h,\xi )$ is the non-hydrostatic pressure

$$p={\displaystyle \frac{3\xi }{h}}\left(1-{\displaystyle \frac{\xi }{h}}\right),$$

$g=9.81$ is the gravitational acceleration, and $\lambda$ is a large pressure constant parameter that makes the system (66) formally converges to the elliptic Serre–Green–Naghdi equations [97] when $\lambda \to \infty .$ Finally, $\xi (x,y,t)$ is an auxiliary equilibrium variable introduced in the relaxation procedure. In the following, we will set $\lambda =300$ for the numerical test considered.

The system satisfied an extra energy conservation law, and the eigenstructure is fully well known. For more details, the reader is refereed to [78]. As for Shallow-water equations, we will consider the family of stationary lake at rest solutions given by

$$u=v=w=0,\eta =cst,h\xi =h^2,$$

(67)

where $\eta =h-H$ is the free-surface elevation.

There are many crucial reasons for considering the numerical validation of the proposed methods for that system here: we will validate the proposed numerical methodology for dispersive shallow flows numerical solutions, and the case of two spatial dimensions; the system introduces a non-homogeneous geometrical source term $G\left(U\right)$; the usual steady-states lake at rest solutions for this system contains nonlinear relations between the conserved variables, and that makes standard well-balancing techniques no longer valid to preserve them. Finally, up to our knowledge, this is the first time that a well-balanced high-order method has been proposed for this particular dispersive hyperbolic system.

5.4.1. Preserving Stationary Solutions

In this first numerical experiment, we set as initial condition a stationary solution of the form (67) to check the well-balancing abilities of the proposed methods for the system (66). The initial condition considered (see Figure 12) is given by

$$u=v=w=0,h=H,h\xi =h^2,H=1-0.2e^{-2(x^2+y^2)}.$$

(68)

The computational domain is $\Omega =[-5,5]\times [-5,5]$ and periodic boundary conditions are imposed. As it can be seen in Figure 13 and Figure 14 and

Table 6. Numerical validation of the well-balanced and non well-balanced methods for the stationary problem (68) (dispersive Gavrilyuk–Favrie equations). Table contains $L_1$ errors for the water depth variable h at a final time $t=10$ s for Runge–Kutta and ADER DG schemes of order $N+1=1,2,3,4.$ Uniform Cartesian meshes of $N_x\times N_y$ elements has been used.

		Non-Well-Balanced		Well-Balanced
N	${\mathit{N}}_{\mathit{x}}\times {\mathit{N}}_{\mathit{y}}$	Runge–Kutta	ADER	Runge–Kutta	ADER
0	20 × 20	$8.06\times {10}^{-01}$	$1.37\times {10}^{-01}$	0	0
	30 × 30	$5.80\times {10}^{-01}$	$3.72\times {10}^{-01}$	0	0
	40 × 40	$4.25\times {10}^{-01}$	$5.36\times {10}^{-01}$	0	0
	50 × 50	$2.28\times {10}^{-01}$	$2.21\times {10}^{-01}$	0	0
1	20 × 20	$1.48\times {10}^{-02}$	$1.71\times {10}^{-02}$	0	0
	30 × 30	$4.20\times {10}^{-03}$	$4.82\times {10}^{-03}$	$5.30\times {10}^{-15}$	$3.31\times {10}^{-14}$
	40 × 40	$1.75\times {10}^{-03}$	$2.01\times {10}^{-03}$	0	0
	50 × 50	$8.68\times {10}^{-04}$	$1.02\times {10}^{-03}$	$1.11\times {10}^{-14}$	$7.56\times {10}^{-14}$
2	20 × 20	$2.38\times {10}^{-03}$	$2.62\times {10}^{-03}$	0	0
	30 × 30	$8.46\times {10}^{-04}$	$8.85\times {10}^{-04}$	$2.77\times {10}^{-14}$	$1.24\times {10}^{-13}$
	40 × 40	$3.70\times {10}^{-04}$	$3.80\times {10}^{-04}$	0	0
	50 × 50	$1.87\times {10}^{-04}$	$1.92\times {10}^{-04}$	$4.45\times {10}^{-14}$	$1.97\times {10}^{-13}$
3	20 × 20	$3.52\times {10}^{-04}$	$4.09\times {10}^{-04}$	0	0
	30 × 30	$4.30\times {10}^{-05}$	$4.60\times {10}^{-05}$	$2.05\times {10}^{-13}$	$2.76\times {10}^{-13}$
	40 × 40	$9.03\times {10}^{-06}$	$9.75\times {10}^{-06}$	$1.97\times {10}^{-13}$	$2.97\times {10}^{-13}$
	50 × 50	$2.82\times {10}^{-06}$	$3.23\times {10}^{-06}$	$3.57\times {10}^{-13}$	$3.62\times {10}^{-13}$

, the well-balanced methods are able to preserve the stationary solution at a final time $t=10$ s for both time-marching discretization and relative coarse meshes. Note that the solutions are compared against the stationary solution (68) as reference.

5.4.2. Limited Simulation

Next, we check the behaviour of the proposed mixed-order technique described in this paper. To do that, we set the following circular dam-break initial condition:

$$u=v=w=0,h\xi =h^2,H=1-0.2e^{-2x^2},$$

$$h=\left\{\begin{array}{cc}H,\hfill & \hfill \mathrm{for}{x}^2+y^2\ge 1\\ 0.25+H,\hfill & \hfill \mathrm{for}{x}^2+y^2<1.\end{array}\right.$$

(69)

The computational domain $\Omega =[-5,5]\times [-5,5]$ is covered with a mesh of $N_x\times N_y=200\times 200$ elements. We set periodic boundary conditions everywhere. Some snapshots of a section at $y=0$ of the computed free-surface are depicted in Figure 15 with a fourth-order Runge–Kutta time marching discretization. As it can be observed, the comparison between a limited and a non-limited simulation shows that there is no oscillations at the neighbourhood of the dispersive shock waves. Note that the high-oscillatory and over-amplified waves appearing in Figure 15 from the unlimited numerical solution, are not of the nature of this hyperbolic dispersive system.

6. Conclusions

We have presented a novel and general technique for preserving stationary solutions for Discontinuous Galerkin methods. This general approach is able to transform standard DG schemes in exactly well-balanced arbitrary high-order Discontinuous Galerkin methods for both Runge–Kutta and ADER time-marching discretisations.

For continuous source terms, the methodology only uses a well-suited definition of a well-balanced reconstruction operator and not on the Riemann-solver. Therefore, the well-balanced technique does not rely on the particular details of the system. That makes the method very general and can be applied to any system of hyperbolic PDEs, as long as the desired stationary solutions to be maintained are known.

Furthermore, the standard WENO limiting procedure is carefully considered. Here, we extend classical ideas to propose a novel WENO limiter that preserve stationary solutions. That makes the final numerical methods proposed here arbitrary high-order, essentially non-oscillatory and well-balanced.

We have tested the numerical methodology derived here with some standard systems, proving the well-balancing properties of the approach, the high-order accuracy, and the general well behaviour of the technique in the presence of strong discontinuities.