# An Energy-Based Discontinuous Galerkin Method with Tame CFL Numbers for the Wave Equation

We extend and analyze the energy-based discontinuous Galerkin method for second order wave equations on staggered and structured meshes. By combining spatial staggering with local time-stepping near boundaries, the method overcomes the typical numerical stiffness associated with high order piecewise polynomial approximations. In one space dimension with periodic boundary conditions and suitably chosen numerical fluxes, we prove bounds on the spatial operators that establish stability for CFL numbers c Δ t/h < C independent of order when stability-enhanced explicit time-stepping schemes of matching order are used. For problems on bounded domains and in higher dimensions we demonstrate numerically that one can march explicitly with large time steps at high order temporal and spatial accuracy.

## Authors

• 11 publications
• 78 publications
• 4 publications
• 5 publications
• ### An energy-based discontinuous Galerkin method for semilinear wave equations

We generalize the energy-based discontinuous Galerkin method proposed in...
01/27/2020 ∙ by Daniel Appelo, et al. ∙ 0

• ### Hermite-Discontinuous Galerkin Overset Grid Methods for the Scalar Wave Equation

We present high order accurate numerical methods for the wave equation t...
01/14/2020 ∙ by Oleksii Beznosov, et al. ∙ 0

• ### Energy stable Runge-Kutta discontinuous Galerkin schemes for fourth order gradient flows

We present unconditionally energy stable Runge-Kutta (RK) discontinuous ...
01/01/2021 ∙ by Hailiang Liu, et al. ∙ 0

• ### Energy-based discontinuous Galerkin difference methods for second-order wave equations

We combine the newly-constructed Galerkin difference basis with the ener...
05/04/2021 ∙ by Lu Zhang, et al. ∙ 0

• ### A short note on the accuracy of the discontinuous Galerkin method with reentrant faces

We study the convergence of the discontinuous Galerkin (DG) method appli...
04/16/2021 ∙ by Will Pazner, et al. ∙ 0

• ### Hybridizable Discontinuous Galerkin Methods for Helmholtz Equation with High Wave Number. Part I: Linear case

This paper addresses several aspects of the linear Hybridizable Disconti...
04/30/2020 ∙ by Bingxin Zhu, et al. ∙ 0

• ### Efficient Explicit Time Stepping of High Order Discontinuous Galerkin Schemes for Waves

This work presents algorithms for the efficient implementation of discon...
05/09/2018 ∙ by Svenja Schoeder, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Discontinuous Galerkin methods [cockburn1989tvb, HesthavenWarburton02] have emerged as one of the most popular discretization techniques for simulating many physical and engineering phenomena including various linear and nonlinear wave models. The methods have excellent dispersive properties, are geometrically flexible, do not have a global mass matrix, can be implemented at any order of accuracy, and, being Galerkin methods, have robust stability properties.

Although discontinuous Galerkin methods are spectrally convergent with the order of the approximation, very high order methods, say , are seldom used in practice. The primary reason for this is that the spectral radius of the discrete spatial derivative operator grows as , where

is an element length scale. This rapidly growing numerical stiffness forces the use of excessively small time steps, effectively prohibiting the use of very high order methods. The source of this numerical stiffness is the approximation by polynomials which must be sampled throughout an element. Heuristically this can be understood by comparing a wave

and its times larger derivative with a typical orthogonal polynomial, say a Chebyshev polynomial, and its derivative . Clearly for and for all , as for the wave, but the derivative , is times larger at the edges.

This numerical stiffness is particularly troublesome for linear wave propagation problems where solutions are typically smooth and do not form discontinuities. Fortunately this numerical stiffness can be circumvented in several ways, for example by allowing the polynomial approximations to the solution to spread out over many elements as in the traditional finite difference methods or as in the more recent Galerkin difference methods [BANKS2016310], or by only sampling the derivatives at the cell center as in Hermite methods [secondHermite]. It is also possible to remove this numerical stiffness within the discontinuous Galerkin framework either by co-volume filtering as proposed in [TAMECFL] or by carrying two approximate solutions on staggered grids as in central discontinuous Galerkin methods [liu20082]. Central discontinuous Galerkin methods combine features of discontinuous Galerkin methods and central schemes [nessyahu1990non], and they are shown in [liu20082], via Fourier analysis, to allow a larger time step size than standard upwind discontinuous Galerkin methods (except for the lowest first order case), when applied to the linear advection equation. In [fli]

these results were established quantitatively for (upwind) discontinuous Galerkin methods and central discontinuous Galerkin methods by estimating the dependence of the operator bound of the respective spatial discretizations on the approximation order.

In this paper we review and extend the energy-based discontinuous Galerkin method [Upwind2] to a staggered grid. For the resulting staggered formulation of the energy-based discontinuous Galerkin method we prove, for problems in one space dimension and with periodic boundary conditions, that for upwind fluxes the norm of the semi-discrete-in-space operator generally grows quadratically in the approximation orders, and , of the displacement and velocity, just as in the non-staggered formulation. When the numerical fluxes are either purely central or upwinded with order-dependent upwinding parameters, however, the result greatly improves for the staggered formulation, in that the norm of the spatial operator grows linearly in and . This, in combination with the Kreiss-Wu theory [kreiss1993stability], indicates that the Courant-Friedrich-Levy (CFL) number is constant independent of order of accuracy as long as high order locally stable time stepping methods with large stability domains are applied. Such methods can be constructed at arbitrary order by adding additional stages to enhance the stability of standard methods; see, for example, [JolyRodriguezLeapFrog] where stability-enhanced leap-frog schemes are proven to exist at all orders and optimized at orders up to sixteen.

We note that, in contrast with [TAMECFL, liu20082] which deal with first order systems, our method focuses on second order wave equations. In addition, again in contrast with those, no additional fields are required. Thus, computing the time derivatives for the method we propose here has the same cost as a standard upwind DG method of the same order implemented on a non-staggered mesh. As a result we attain the full efficiency afforded by the larger time steps.

At physical boundaries it is no longer possible to stagger the mesh and the CFL constraint becomes order dependent again. Fortunately we present numerical experiments illustrating that the local timestepping methods of Diaz and Grote [diaz2009energy] can be used in a thin layer of elements near the boundary and that this approach allows us to retain the large time steps in the interior. The resulting method, while having some additional computational overhead near boundaries, has the same computational complexity as for the periodic case.

The rest of the paper is organized as follows. In section 2 we review the formulation of energy-DG methods for the scalar wave equation and extend them to staggered, structured meshes. In section 3 we establish bounds on the norm of the spatial operator in one space dimension and with periodic boundary conditions. In section 4 we briefly discuss our time-stepping schemes and the corrections needed to maintain large time steps in the presence of boundaries. Lastly, in section 5 we demonstrate the accuracy and stability of the method in two space dimensions by means of some simple numerical experiments.

## 2 Energy-Based Discontinuous Galerkin Method for the Wave Equation

We consider the scalar wave equation written as a first order system in time

 ∂u(x,t)∂t=v(x,t), (1) ∂v(x,t)∂t=∇⋅(c2(x)∇u(x,t))+f(x,t), (2)

on the domain

 x=(x1,…,xd)∈Ω⊂Rd,t>0,

with initial conditions

 u(x,0)=g(x),  ∂u∂t(x,0)=v(x,0)=h(x), (3)

and boundary conditions

 γ∂u∂t+κc→n⋅∇u=0,  x∈∂Ω. (4)

Here is the speed of sound and is the outward pointing unit normal. For the boundary conditions we assume the normalization and that . Then the choice corresponds to a homogeneous Dirichlet boundary condition on and corresponds to a homogeneous Neumann boundary condition. Any choice with being positive will dissipate the energy of the system and can be thought of as a low order non-reflecting boundary condition.

The energy associated with the scalar wave equation is

 E(t)=12∫Ω(∂u∂t)2+c2(x)|∇u|2dx, (5)

and it is a discrete version of this energy that our energy-based discontinuous Galerkin method is built from.

We now present the non-staggered and staggered formulations of the method. A more thorough analysis of the non-staggered method can be found in [Upwind2], but we include it here to illustrate the differences between the two formulations.

### 2.1 Non-staggered Formulation

Let the finite element mesh, , with

 Ω=⋃jΩj,

be a discretization of consisting of geometry-conforming and non-overlapping (possibly curved) elements with piecewise smooth boundaries.

On each element , the approximation to the displacement, , and the approximation to the velocity, , are elements of some finite dimensional spaces and respectively. Then, the non-staggered energy-based discontinuous Galerkin method can be stated as follows. On each element , require that for all test functions

 ϕ∈Uh,  ψ∈Vh,

the following variational formulation holds:

 ∫Ωjc2∇ϕ⋅(∂∇uh∂t−∇vh)dx = ∫∂Ωj(c2∇ϕ⋅→n)(v∗−vh)ds, (6a) ∫Ωjψ∂vh∂t+∇ψ⋅(c2∇uh)dx = ∫∂Ωjψ((c2∇u)⋅→n)∗ds. (6b)

As described in [Upwind2] the energy is invariant to constants and this necessitates an additional equation complementing (6a)

 ∫Ωj(∂uh∂t−vh)dx=0,∀j. (7)

Here and are numerical fluxes computed from the averages and jumps of function values and derivatives. Arbitrarily labeling values from adjacent elements and we recall the standard notation:

 {{vh}}α = 12(αvh,1+(1−α)vh,2) [[vh]] = vh,1→n1+vh,2→n2, (8) {{c2∇uh}}α = 12(αc2∇uh,1+(1−α)c2∇uh,2), [[c2∇uh]] = c2∇uh,1⋅→n1+c2∇uh,2⋅→n2. (9)

We then set

 v∗={{vh}}α−β[[c2∇uh]], (10)
 (c2∇u⋅→n)∗={{c2∇uh}}1−α⋅→n−τ[[vh]]⋅→n. (11)

Here is an upwinding parameter with units of and is an upwinding parameter with units of . When , one can recover the commonly used central fluxes by choosing , and alternating fluxes with or .

### 2.2 Staggered Formulation

We now consider two structured finite element meshes, and

 Ω=⋃jΩj=⋃kΩ⋄k.

We assume each mesh consists of geometry-conforming and non-overlapping (possibly curved) quadrilaterals (or hexahedra) with piecewise smooth boundaries. We assume that the meshes are staggered. More precisely, away from non-periodic boundaries we assume that all quadrilaterals (hexahedra) are straight sided and convex and that all vertices have valence 4 (6). By staggering we mean that, away from boundaries, the vertices of the mesh coincide with the centers (defined as the vertex, side or area/volume centroid) of the elements in .

For consistency with the theoretical and computational results to follow, we take the approximation to the velocity, , restricted to an element in

, to be a tensor product polynomial in

while the approximation to the displacement, , restricted to an element in , is taken to be a tensor product polynomial in . Here , .

The staggered energy-based discontinuous Galerkin method then can be stated as follows. On each element and , require that for all test functions

 ϕ∈(Qqu(Ωj))d,  ψ∈(Qqv(Ω⋄k))d,

the following variational formulation holds:

 ∫Ωj(c2∇ϕ⋅∂∇uh∂t+∇⋅(c2∇ϕ)vh)dx = ∫∂Ωj(c2∇ϕ⋅→n)v∗ds, (12a) ∫Ω⋄kψ∂vh∂t+∇ψ⋅(c2∇uh)dx = ∫∂Ω⋄kψ((c2∇u)⋅n)∗ds. (12b)

As with the non-staggered formulation we must complement (12a) with the equation

 ∫Ωj(∂uh∂t−vh)dx=0,∀j. (13)

Again, here and are numerical fluxes as in (10)-(11). However, taking account of the staggering, we note that is single valued at and is single valued at so the choice of is not relevant. Lastly we note that the integrals of gradients in the variational form as well as in the calculations below are understood to be piecewise-defined in subdomains where the functions are smooth. For example, the integral in (12b) includes boundaries of elements in across which is discontinuous. We do not, here, interpret in a distributional sense and so no additional boundary terms are implied.

We note the difference between (6a) and (12a). If the term is integrated by parts in (12a), terms involving the jump in across boundaries of dual mesh elements will appear. These play a role in the energy estimate we now derive.

Define the discrete energy to be

 Eh(t)=12∑k∫Ω⋄k(vh)2 dx+12∑j∫Ωjc2⏐∇uh⏐2 dx. (14)

To start, we assume periodic boundary conditions. Choosing in (12a), integrating by parts, and using the fact that is single valued on we find

 12ddt∑j∫Ωjc2⏐∇uh⏐2 dx = ∑j∫Ωjc2∇uh⋅∇vh dx−∑k∫∂Ω⋄kc2∇uh⋅[[vh]] ds −β∑j∫∂Ωj[[c2∇uh]]2 ds,

Similarly we find

 12ddt∑k∫Ω⋄k(vh)2 dx = −∑k∫Ω⋄kc2∇uh⋅∇vh dx+∑k∫∂Ω⋄kc2∇uh⋅[[vh]] ds −τ∑k∫∂Ω⋄k⏐[[vh]]⏐2 ds.

Summing these equations, we see that the left-hand side is the time derivative of the discrete energy. Since and recalling the piecewise definition of the integrals we conclude that the terms involving cancel. Thus we conclude:

 dEhdt=−β∑j∫∂Ωj[[c2∇uh]]2 ds−τ∑k∫∂Ω⋄k⏐[[vh]]⏐2 ds. (15)

At nonperiodic boundaries we alter the staggered mesh so that elements from both and conform to . Now the imposition of the boundary conditions is the same as for the non-staggered formulation. For example, recalling (4) we may set

 v∗ = κγ+κ(vh−c∇uh⋅→n) dx, (16) (c2∇uh⋅→n)∗ = γγ+κ(c2∇uh⋅→n−cvh). (17)

Then the contribution of the nonperiodic boundaries to the energy derivative can be shown to be nonpositive. The mesh modification at these boundaries will preclude taking global large time steps. To maintain the efficiency of the staggered scheme we will then use local time stepping in the vicinity of the boundaries; see section 4 for details.

## 3 Operator Bounds on Periodic Domains in One Space Dimension

In this section we use the techniques from [TAMECFL, fli] to establish bounds for the energy-based DG and the staggered energy-based DG spatial operator for the second-order wave equation (18) and (19) in one space dimension. This allows us to invoke the Kreiss-Wu theory [kreiss1993stability] combined with the energy estimates to establish the stability of fully discrete locally stable explicit time-stepping schemes.

We restrict the analysis to uniform grids, periodic boundary conditions, and constant coefficients. As the key ingredient to taming the CFL condition is to evaluate certain terms with derivatives only near the element centers, we expect that the analysis can be extended to smoothly varying grids and to variable coefficients. The numerical experiments demonstrate the efficiency of the method for a variable coefficient problem. It may also be possible to extend the analysis to problems with Dirichlet or Neumann boundary conditions by using the image principle, however we don’t pursue this here.

### 3.1 Operator Bounds for the Non-staggered Formulation

Now, consider the one dimensional wave equation in a uniform medium

 ut = v, (18) vt = c2uxx, (19)

on the domain . Let the domain be discretized by a grid , and . Associated with the grid, we define two broken finite element spaces

 Uquh={w:w|Ij∈Qqu(Ij)∀j},    Vqvh={w:w|Ij∈Qqv(Ij)∀j}.

Here and below is the space of polynomials of degree up to in , , and . In addition, we denote .

The energy-based DG scheme then consists of finding and such that for any and and for all

 ∫xj+1xjc2ϕx(∂uhx∂t−vhx)dx = c2ϕ−x(v∗−vh,−)∣∣xj+1−c2ϕ+x(v∗−vh,+)∣∣xj, (20a) ∫xj+1xjψ∂vh∂t+c2ψxuhxdx = c2ψ−u∗x∣∣xj+1−c2ψ+u∗x∣∣xj, (20b) ∫xj+1xj∂uh∂t−vhdx = 0. (20c)

Assuming periodic boundary conditions we may add up the equations (20a)-(20b) in to find

 ∫Ωc2ϕx∂uhx∂t+ψ∂vh∂tdx=∫Ωc2ϕxvhx−c2ψxuhxdx (21) +∑jc2(ϕ−x(v∗−vh,−)∣∣xj+1−ϕ+x(v∗−vh,+)∣∣xj+ψ−u∗x∣∣xj+1−ψ+u∗x∣∣xj).

Throughout, the spatial derivative of functions in any broken finite element space shall be understood as being defined element by element. To connect the element solutions in a stable fashion we use the numerical fluxes defined in (10)-(11) and introduce the notation

 v∗ = H(vh,uhx)={{vh}}α−βc2(uh,−x−uh,+x), (22a) u∗x = G(uhx,vh)={{uhx}}1−α−τc2(vh,−−vh,+). (22b)

Then the energy estimate (15) holds. We note that it can also be used to establish error estimates for different choices of , and ; see [Upwind2] for details.

We now establish bounds on the spatial operators which constrain the allowable time step sizes for explicit marching schemes. In particular we are interested in the dependence of these bounds on the approximation orders, and and will follow a similar analysis as in [TAMECFL, fli]. With the choice of the numerical fluxes in (10)-(11), an important observation is that the first two equations in (20) are coupled with (20c) in a one-way manner. That is, (20a)-(20b) will uniquely determine and . Once , are available, one can further recover the missing constant in on (i.e. in the form of the cell average of ) through (20c) for all . As this last step is simply an integration in time it can not affect the numerical stability; see also the discussion in section 4.

These considerations motivate us to define the operator ,

 ∫ΩL(cw,v)(cφ,ψ)dx= c2∫Ωφvx−ψxwdx+c2∑j(ψ−G(w,v)∣∣xj+1−ψ+G(w,v)∣∣xj) (23) +c2∑j(φ−(H(v,w)−v−)∣∣xj+1−φ+(H(v,w)−v+)∣∣xj)

for any and with the operator norm as

 ∥L∥≡supw,φ∈Uqu−1h,v,ψ∈Vqvh(w,v)≠(0,0),(φ,ψ)≠(0,0)∫ΩL(cw,v)(cφ,ψ)dx(∥cw∥2L2(Ω)+∥v∥2L2(Ω))1/2(∥cφ∥2L2(Ω)+∥ψ∥2L2(Ω))1/2. (24)

Once the bound is established for , time step condition, , for method of lines discretization combined with locally stable one-step temporal methods will follow from Kreiss-Wu theory [kreiss1993stability]. Here is defined as the radius of the largest semidisk in the closed left half complex plane contained in the stability domain of the method. It is well known, [ketcheson2015absolute], that one-step methods based on Taylor expansion with terms are locally stable. For of moderate size they have stability domains which grow with order. Thus, if we can establish a bound on that grows linearly in and , we should expect that the fully discrete method can time-march at a CFL condition of when the spatial and temporal orders are matched. As the order increases this does not hold, but, as discussed in section 4, with the introduction of additional stages the size of the domain can be greatly increased. In what follows we will see that such a bound can be established for the staggered method with suitably chosen numerical fluxes but not for the non-staggered method. For the non-staggered method, the bound on is quadratic in and , and this will be proved next.

###### Theorem 1.

Let the energy-based DG spatial operator be defined as in (23) with periodic boundary conditions and with numerical fluxes defined by (10)-(11). Then the following estimate holds:

 ∥L∥ ≤ C1chmax{(qu−1)2,q2v} +C2ch(2max{cβq2u,τc(qv+1)2}+(|α|+|1−α|)(qv+1)qu).

Here and are two universal positive constants, independent of , , , and .

###### Proof.

Consider any and . Applying element-wise integration by parts and the triangle inequality, we have

 ∫ΩL(cw,v)(cφ,ψ)dx= (26) +c2∑j(ψ−G(w,v)∣∣xj+1−ψ+G(w,v)∣∣xj)≤Λ1+Λ2+Λ3,

with

 Λ1 = c2∑j∫xj+1xj|φxv|dx+∫xj+1xj|ψxw|dx, Λ2 = c2∑j∣∣φ−H(v,w)∣∣∣∣xj+1+∣∣φ+H(v,w)∣∣∣∣xj, Λ3 = c2∑j∣∣ψ−G(w,v)∣∣∣∣xj+1+∣∣ψ+G(w,v)∣∣∣∣xj.

We now bound each of the terms, starting with the volume term . By applying the Cauchy-Schwarz inequality, we have

 Λ1≤c∑j∥cφx∥L2(Ij)∥v∥L2(Ij)+∥ψx∥L2(Ij)∥cw∥L2(Ij).

For and , we use the definitions of and as well as the triangle inequality, and arrive at

 Λ2≤c∑j ∣∣cφ−∣∣(|α|∣∣v−∣∣+|1−α|∣∣v+∣∣+cβ∣∣cw−∣∣+cβ∣∣cw+∣∣)∣∣xj+1 + ∣∣cφ+∣∣(|α|∣∣v−∣∣+|1−α|∣∣v+∣∣+cβ∣∣cw−∣∣+cβ∣∣cw+∣∣)∣∣xj, Λ3≤c∑j ∣∣ψ−∣∣(|1−α|∣∣cw−∣∣+|α|∣∣cw+∣∣+τc∣∣v−∣∣+τc∣∣v+∣∣)∣∣xj+1 + ∣∣ψ+∣∣(|1−α|∣∣cw−∣∣+|α|∣∣cw+∣∣+τc∣∣v−∣∣+τc∣∣v+∣∣)∣∣xj.

Now we recall some standard inverse inequalities for polynomials spaces [fli]; there exist positive constants , , such that ,

 ∥px∥L2([−1,1])≤^C1k2∥p∥L2([−1,1]),p(x)≤^C2(k+1)∥p∥L2([−1,1])∀x∈[−1,1]. (27)

By applying these inverse inequalities, with a linear scaling from to , and Cauchy-Schwarz inequality, we find that

 Λ1 ≤2^C1ch∑j((qu−1)2∥cφ∥L2(Ij)∥v∥L2(Ij)+q2v∥ψ∥L2(Ij)∥cw∥L2(Ij)) ≤C1chmax{(qu−1)2,q2v}(∥cw∥2L2(Ω)+∥v∥2L2(Ω))1/2(∥cφ∥2L2(Ω)+∥ψ∥2L2(Ω))1/2,

and similarly, using a linear scaling from to (with ), we have

 Λ2≤ 2^C22ch∑j(qu(qv+1)∥cφ∥L2(Ij)(|α|∥v∥L2(Ij−1∪Ij)+|1−α|∥v∥L2(Ij∪Ij+1)) +cβq2u∥cφ∥L2(Ij)(2∥cw∥L2(Ij)+∥cw∥L2(Ij−1∪Ij+1))),
 Λ3≤ 2^C22ch∑j(qu(qv+1)∥ψ∥L2(Ij)(|1−α|∥cw∥L2(Ij−1∪Ij)+|α|∥cw∥L2(Ij∪Ij+1)) +τ(qv+1)2c∥ψ∥L2(Ij)(2∥v∥L2(Ij)+∥v∥L2(Ij−1∪Ij+1))),

and hence

 Λ2+Λ3≤ C2ch(2max{cβq2u,τc(qv+1)2}+(|α|+|1−α|)(qv+1)qu) (∥cw∥2L2(Ω)+∥v∥2L2(Ω))1/2(∥cφ∥2L2(Ω)+∥ψ∥2L2(Ω))1/2.

Finally by adding up , and based on the operator norm in (24), we reach the bound in (1). ∎

### 3.2 Operator Bounds for the Staggered Formulation

For the staggered version of the method we introduce element centers as well as the staggered grid composed of the elements . Associated with both grids, we define two broken finite element spaces

 Uquh={w:w|Ij∈Qqu(Ij)∀j},    ˜Vqvh={w:w|Ij+12∈Qqv(Ij+12)∀j}.

The staggered energy-based DG scheme then consists of finding and such that for any and and for all

 ∫xj+1xjc2ϕx∂uhx∂tdx+∫ρjxjc2ϕxx~vhdx+∫xj+1ρjc2ϕxx~vhdx∫xj+1xjc2ϕxx~vhdx =c2ϕ−xv∗∣∣xj+1−c2ϕ+xv∗∣∣xj, (28a) ∫ρj+1ρjψ∂~vh∂tdx+∫xj+1ρjc2ψxuhxdx+∫ρj+1xj+1c2ψxuhxdx∫ρj+1ρjc2ψxuhxdx =c2ψ−u∗x∣∣ρj+1−c2ψ+u∗x∣∣ρj, (28b) ∫xj+1xj(∂uh∂t−~vh)dx=0. (28c)

Note that the second and third integrals in (28a) (resp. in (28b)) are against (resp. ) from two elements.

Explicitly we write the flux terms

 v∗ = H(~vh,uhx)=~vh