An efficient ADER discontinuous Galerkin scheme for directly solving Hamilton-Jacobi equation1) footnote 1) 1) footnote 1) It is an invited paper for the special issue on “The Sixth Chinese-German Workshop on Computational and Applied Mathematics”.

# An efficient ADER discontinuous Galerkin scheme for directly solving Hamilton-Jacobi equation1) 1) 1) It is an invited paper for the special issue on “The Sixth Chinese-German Workshop on Computational and Applied Mathematics”.

Junming Duan LMAM, School of Mathematical Sciences, Peking University, Beijing 100871, China    Huazhong Tang HEDPS, CAPT & LMAM, School of Mathematical Sciences, Peking University, Beijing 100871; School of Mathematics and Computational Science, Xiangtan University, Hunan Province, Xiangtan 411105, P.R. China
###### Abstract

This paper proposes an efficient ADER (Arbitrary DERivatives in space and time) discontinuous Galerkin (DG) scheme to directly solve the Hamilton-Jacobi equation. Unlike multi-stage Runge-Kutta methods used in the Runge-Kutta DG (RKDG) schemes, the ADER scheme is one-stage in time discretization, which is desirable in many applications. The ADER scheme used here relies on a local continuous spacetime Galerkin predictor instead of the usual Cauchy-Kovalewski procedure to achieve high order accuracy both in space and time. In such predictor step, a local Cauchy problem in each cell is solved based on a weak formulation of the original equations in spacetime. The resulting spacetime representation of the numerical solution provides the temporal accuracy that matches the spatial accuracy of the underlying DG solution. The scheme is formulated in the modal space and the volume integral and the numerical fluxes at the cell interfaces can be explicitly written. The explicit formulae of the scheme at third order is provided on two-dimensional structured meshes. The computational complexity of the ADER-DG scheme is compared to that of the RKDG scheme. Numerical experiments are also provided to demonstrate the accuracy and efficiency of our scheme.

An efficient ADER discontinuous Galerkin scheme for directly solving Hamilton-Jacobi equationthanks: Received xxx / Revised version received xxx / Accepted xxx /thanks: It is an invited paper for the special issue on “The Sixth Chinese-German Workshop on Computational and Applied Mathematics”.

Junming Duan
LMAM, School of Mathematical Sciences, Peking University, Beijing 100871, China

and Huazhong Tang
HEDPS, CAPT & LMAM, School of Mathematical Sciences, Peking University, Beijing 100871; School of Mathematics and Computational Science, Xiangtan University, Hunan Province, Xiangtan 411105, P.R. China

Mathematics subject classification: 65M06, 35F21, 70H20.

Key words: Hamilton-Jacobi equation, ADER, discontinuous Galerkin methods, local continuous spacetime Galerkin predictor, high order accuracy.

## 1 Introduction

Consider the Hamilton-Jacobi (HJ) equation

 φt+H(∇xφ,x)=0,φ(x,0)=φ0(x),x∈Ω∈Rd, (1.1)

with suitable boundary conditions, where denotes the Hamiltonian. The HJ equations are used in many application areas, such as optimal control theory, geometrical optics, crystal growth, image processing and computer vision. The solutions of such equations are continuous but their derivatives could be discontinuous even if the initial condition is smooth. Viscosity solutions were firstly introduced and studied in [6, 7], which are the unique physically relevant solutions.

It is well known that the HJ equations are closely related to hyperbolic conservation laws, thus many successful numerical methods for the conservation laws can be adapted for solving the HJ equations. In [7], a monotone finite difference scheme was introduced and proved to be convergent to the viscosity solution. A second order finite difference essentially non-oscillatory (ENO) scheme was developed in [16], and then a higher-order weighted ENO (WENO) scheme is proposed in [14]. Tang and his collaborators [18] developed an adaptive mesh redistribution method for the HJ equations in two- and three-dimensions. Qiu et al. [17] developed the Hermite WENO (HWENO) schemes of the HJ equations. The high order finite difference WENO scheme on unstructured meshes was developed in [21], but its implementation is a bit complicated.

Alternatively, a DG method was designed in [13] to solve the HJ equations, and its reinterpretation and simplified implementation was given in [15]. Those DG methods were based on the fact that the derivatives of the solution satisfied the conservation laws. It was correct in the one-dimensional case but at risk in the multi-dimensional case because corresponding multi-dimensional conservation laws is only weakly hyperbolic in general. Later, a DG method for directly solving the HJ equations with convex Hamiltonians was proposed in [3]. It was further improved and a new DG method was derived for directly solving the general HJ equations with nonconvex Hamiltonians in [4]. This paper will construct the scheme based on the RKDG scheme in [4]. The RKDG method [5] was originally designed to solve conservation laws, which has the advantages of flexibility on complicated geometries and a compact stencil, and is easy to obtain high order accuracy.

Most of the above methods use the multi-stage Runge-Kutta time discretization, thus have the advantage of simplicity but are time-consuming because at each stage, the volume integration and the numerical fluxes at cell interfaces have to be calculated and the nonlinear limiters should be performed to suppress the numerical oscillations. Thus, in order to save the computational cost, it is desirable to use an alternative to the multi-stage Runge-Kutta method. One choice is the Lax-Wendroff type time discretization, which converts all (or partial, when approximations with certain accuracy are expected) time derivatives in a temporal Taylor expansion of the solution into spatial derivatives by repeatedly using the underlying differential equation and its differentiated forms [12]. In [12], a local-structure-preserving DG method with Lax-Wendroff type time discretization was proposed for solving the HJ equations. It is shown that such method is relatively more efficient than the RKDG method in [15]. But the Cauchy-Kowalewski procedure may become a little cumbersome when we want to construct a high order scheme. This paper will use the time discretization (named ADER) proposed in [8, 9]. The ADER scheme has been successfully applied to the (magneto) hydrodynamics and relativistic (magneto) hydrodynamics with stiff or non-stiff source terms [1, 2, 8, 9, 11]. It is based on a local spacetime Galerkin predictor step, at which a local Cauchy problem is solved in each cell, based on a weak formulation of the original partial differential equations in spacetime. Through the above procedure, the resulting spacetime representation of the numerical solution provides the temporal accuracy that matches the spatial accuracy of the underlying DG solution. The ADER scheme is a one-step one-stage time discretization, which means that the volume integration and the numerical fluxes terms at cell interfaces are only calculated once at each time step. Our ADER-DG scheme is formulated in modal space. Thanks to the spacetime representation of the numerical solution, we can write down explicit formulae of the scheme using the strategy presented in [1], and we will provide the implementation details of the scheme at third order on two-dimensional structured meshes. Our ADER-DG scheme can capture the viscosity solution accurately and efficiently, and will be validated by the analysis of the computational complexity and the numerical experiments.

The paper is organized as follows. Section 2 presents the general formulation of our one- and two-dimensional ADER-DG schemes. Section 3 introduces the local spacetime continuous Galerkin predictor, and gives a detailed description of the two-dimensional predictor step at third order. Section 4 describes the calculation of the volume integration and the numerical fluxes terms at the cell interfaces, and the computational complexity of our ADER-DG scheme will be compared to the RKDG scheme in [4]. Section 5 presents numerical experiments and the concluding remarks are given in Section 6.

## 2 General formulation of the ADER-DG schemes

This section will present the general formulation of the ADER-DG schemes, in which the numerical fluxes and the penalty terms adding to the numerical fluxes are firstly developed in [4].

Let us consider the one-dimensional HJ equation at first. In this case, (1.1) becomes

 φt+H(φx,x)=0,φ(x,0)=φ0(x). (2.1)

Assume the computational domain is divided into cells, , , where

 a=x12

Denote the center of and the mesh size as and respectively. Moreover, use to denote the partial derivative of the Hamiltonian with respect to .

The spatial DG approximation space is

 Vkh={v:v|Ii∈Pk(Ii),i=1,2,…,N}, (2.2)

where denotes all polynomials of degree at most on . Assume the current time interval is , and the time stepsize is . Following [4], if multiplying (2.1) with the test function , introducing the numerical fluxes, adding penalty terms for the numerical fluxes at the interfaces of computational cells, and integrating it over the spacetime control volume , then one has

 ∫tn+1tn∫Iiϕm(x)(∂tφ(x,t)+H(∂xφ(x,t),x))dxdt +∫tn+1tnmin(~H1,φ(xi+12,t),0)[φ]i+12(ϕm)−i+12dt+∫tn+1tnmax(~H1,φ(xi−12,t),0)[φ]i−12(ϕm)+i−12dt −∫tn+1tnCΔxi(S1,φ(xi+12,t)−|~H1,φ(xi+12,t)|)[(φ)x]i+12(ϕm)−i+12dt −∫tn+1tnCΔxi(S1,φ(xi−12,t)−|~H1,φ(xi−12,t)|)[(φ)x]i−12(ϕm)+i−12dt =0,∀i=1,2,…,N, (2.3)

where denotes the jump of function at the cell interface, the superscripts denote the right, and left limits of a function, is a positive penalty parameter chosen as in this paper, and and are the Roe speed and the parameter to identify the entropy violating cells [4]. For all , assume is a point located at the cell interface, then and are defined by

 ~H1,φ(x∗,t)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩H(φx(x+∗,t),x+∗)−H(φx(x−∗,t),x−∗)φx(x+∗,t)−φx(x−∗,t),ifφx(x−∗,t)≠φx(x+∗,t),12(H1(φx(x+∗,t),x+∗)+H1(φx(x−∗,t),x−∗)),ifφx(x−∗,t)=φx(x+∗,t), δ1,φ(x∗,t)=max(0,~H1,φ(x∗,t)−H1((φ)x(x−∗,t),x−∗),H1((φ)x(x+∗,t),x+∗)−~H1,φ(x∗,t)), S1,φ(x∗,t)=max(δ1,φ(x∗,t),|~H1,φ(x∗,t)|). (2.4)

It is worth noting that the above definitions only make sense for . That is why we choose the DG space as .

Calculate the time derivative parts in (2.3), restrict the solutions to , thus (2.3) becomes

 ∫Iiϕm(x)(φh(x,tn+1)−φh(x,tn))dx+∫tn+1tn∫Iiϕm(x)H(∂xqh(x,t),x)dxdt +∫tn+1tnmin(~H1,qh(xi+12,t),0)[qh]i+12(ϕm)−i+12dt+∫tn+1tnmax(~H1,qh(xi−12,t),0)[qh]i−12(ϕm)+i−12dt −∫tn+1tnCΔxi(S1,qh(xi+12,t)−|~H1,qh(xi+12,t)|)[(qh)x]i+12(ϕm)−i+12dt −∫tn+1tnCΔxi(S1,qh(xi−12,t)−|~H1,qh(xi−12,t)|)[(qh)x]i−12(ϕm)+i−12dt =0,∀i=1,2,⋯,N, (2.5)

where an element-local spacetime predictor solution is introduced to replace the solution in the integral of the Hamiltonian and the numerical fluxes, which is a high order approximation polynomial obtained by using the local spacetime Galerkin predictor, and will be presented in detail in Section 3. At time , the DG solution is known, if we get , then the DG solution can be evolved to time as from (2.5).

### 2.2 Two-dimensional ADER-DG scheme on structured meshes

Consider the two-dimensional HJ equation

 φt+H(φx,φy,x,y)=0,φ(x,y,0)=φ0(x,y), (2.6)

and the computational domain is divided into cells. with , , , and , , . Use and to denote the partial derivatives of the Hamiltonian with respect to and respectively.

The spatial DG approximation space is

 Vkh={v:v|Ii,j∈Pk(Ii,j),i=1,2,…,Nx,j=1,2,…,Ny}, (2.7)

where denotes all polynomials of degree at most on . Assume the current time interval is , and the time step is . Multiplying (2.6) with test functions , introducing the numerical fluxes and adding penalty terms for the numerical fluxes at the cell interfaces, and then integrating it over the spacetime control volume can give

 ∫tn+1tn∫Ii,jϕm(x,y)(∂tφ(x,y,t)+H(∂xφ(x,y,t),∂yφ(x,y,t),x,y))dxdydt +∫tn+1tn∫Kjmin(~H1,φ(xi+12,y,t),0)[φ](xi+12,y,t)ϕm(x−i+12,y)dydt +∫tn+1tn∫Kjmax(~H1,φ(xi−12,y,t),0)[φ](xi−12,y,t)ϕm(x+i−12,y)dydt +∫tn+1tn∫Jimin(~H2,φ(x,yj+12,t),0)[φ](x,yj+12,t)ϕm(x,y−j+12)dxdt +∫tn+1tn∫Jimax(~H2,φ(x,yj−12,t),0)[φ](x,yj−12,t)ϕm(x,y+j−12)dxdt −CΔxi∫tn+1tn∫Kj(S1,φ(xi+12,y,t)−|~H1,φ(xi+12,y,t)|)[φx](xi+12,y,t)ϕm(x−i+12,y)dydt −CΔxi∫tn+1tn∫Kj(S1,φ(xi−12,y,t)−|~H1,φ(xi−12,y,t)|)[φx](xi−12,y,t)ϕm(x+i−12,y)dydt −CΔyj∫tn+1tn∫Ji(S2,φ(x,yj+12,t)−|~H2,φ(x,yj+12,t)|)[φy](x,yj+12,t)ϕm(x,y−j+12)dxdt −CΔyj∫tn+1tn∫Ji(S2,φ(x,yj−12,t)−|~H2,φ(x,yj−12,t)|)[φy](x,yj−12,t)ϕm(x,y+j−12)dxdt =0,∀i=1,2,…,Nx,j=1,2,…,Ny. (2.8)

For all , if assuming is a point located at the cell interface in the -direction, then the Roe speed and the parameters to identify the entropy violating cells in the scheme are given by

 ~H1,φ(x∗,y,t) =⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩H(φx(x+∗,y,t),¯¯¯¯¯¯φy,x+∗,y)−H(φx(x−∗,y,t),¯¯¯¯¯¯φy,x−∗,y)φx(x+∗,y,t)−φx(x−∗,y,t),φx(x−∗,y,t)≠φx(x+∗,y,t),12(H1(φx(x+∗,y,t),¯¯¯¯¯¯φy,x+∗,y)+H1(φx(x−∗,y,t),¯¯¯¯¯¯φy,x−∗,y)),φx(x−∗,y,t)=φx(x+∗,y,t), δ1,φ(x∗,y,t) =max(0,~H1,φ(x∗,y,t)−H1(φx(x−∗,y,t),¯¯¯¯¯¯φy,x−∗,y),H1(φx(x+∗,y,t),¯¯¯¯¯¯φy,x+∗,y)−~H1,φ(x∗,y,t)), S1,φ(x∗,y,t)=max(δ1,φ(x∗,y,t),|~H1,φ(x∗,y,t)|),

where is the average of the tangential derivative.

Similarly, for all , if denoting as a point located at the cell interface in the -direction, then the Roe speed and the parameters are given by

 ~H2,φ(x,y∗,t) =⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩H(¯¯¯¯¯¯φx,φy(x,y+∗,t),x,y+∗)−H(¯¯¯¯¯¯φx,φy(x,y−∗,t),x,y−∗)φy(x,y+∗,t)−φy(x,y−∗,t),φy(x,y−∗,t)≠φy(x,y+∗,t),12(H2(¯¯¯¯¯¯φx,φy(x,y+∗,t),x,y+∗)+H2(¯¯¯¯¯¯φx,φy(x,y−∗,t),x,y−∗)),φy(x,y−∗,t)=φy(x,y+∗,t), δ2,φ(x,y∗,t) =max(0,~H2,φ(x,y∗,t)−H2(¯¯¯¯¯¯φx,φy(x,y−∗,t),x,y−∗),H2(¯¯¯¯¯¯φx,φy(x,y+∗,t),x,y+∗)−~H2,φ(x,y∗,t)), S2,φ(x,y∗,t)=max(δ2,φ(x,y∗,t),|~H2,φ(x,y∗,t)|),

where .

After calculating the time derivative parts in (2.8) and restricting the solutions to , then (2.8) becomes

 ∫Ii,jϕm(x,y)(φh(x,y,tn+1)−φh(x,y,tn))dxdy +∫tn+1tn∫Ii,jϕm(x,y)H(∂xqh(x,y,t),∂yqh(x,y,t),x,y)dxdydt +∫tn+1tn∫Kjmin(~H1,qh(xi+12,y,t),0)[qh](xi+12,y,t)ϕm(x−i+12,y)dydt +∫tn+1tn∫Kjmax(~H1,qh(xi−12,y,t),0)[qh](xi−12,y,t)ϕm(x+i−12,y)dydt +∫tn+1tn∫Jimin(~H2,qh(x,yj+12,t),0)[qh](x,yj+12,t)ϕm(x,y−j+12)dxdt +∫tn+1tn∫Jimax(~H2,qh(x,yj−12,t),0)[qh](x,yj−12,t)ϕm(x,y+j−12)dxdt −CΔxi∫tn+1tn∫Kj(S1,qh(xi+12,y,t)−|~H1,qh(xi+12,y,t)|)[(qh)x](xi+12,y,t)ϕm(x−i+12,y)dydt −CΔxi∫tn+1tn∫Kj(S1,qh(xi−12,y,t)−|~H1,qh(xi−12,y,t)|)[(qh)x](xi−12,y,t)ϕm(x+i−12,y)dydt −CΔyj∫tn+1tn∫Ji(S2,qh(x,yj+12,t)−|~H2,qh(x,yj+12,t)|)[(qh)y](x,yj+12,t)ϕm(x,y−j+12)dxdt −CΔyj∫tn+1tn∫Ji(S2,qh(x,yj−12,t)−|~H2,qh(x,yj−12,t)|)[(qh)y](x,yj−12,t)ϕm(x,y+j−12)dxdt =0,∀i=1,2,…,Nx,j=1,2,…,Ny, (2.9)

where the element-local spacetime predictor solution will be introduced in Section 3. The remaining task is to give .

## 3 Local spacetime continuous Galerkin predictor

Unlike the classical ADER schemes in [19, 20] using Cauchy-Kovalewski procedure, which may become cumbersome for high order schemes, the new formulation of ADER schemes proposed in [8] is based on a local weak formulation of the governing PDE in spacetime. The new ADER schemes rely on an iterative predictor step to obtain the spacetime representation of the solution within each cell, i.e., the previous mentioned local spacetime predictor solution . This part will construct the predictor step, and give the implementation details of the predictor step in the two-dimensional case at third order.

### 3.1 General formulation of continuous Galerkin predictor

For the sake of convenience, we will only consider the two-dimensional case. Assume the spatial coordinates in the reference element is , and the temporal coordinates in the reference element is . In the reference element, Eq. (2.6) can be written as

 ∂φ∂τ+h(1Δx∂φ∂ξ,1Δy∂φ∂η,ξ,η)=0, (3.1)

where , and are the mesh sizes of the cell. The ADER scheme used here is a modal variant of the ADER scheme with a continuous Galerkin representation in time described in [8]. Assume that there are spacetime basis functions in the reference element, . The continuous Galerkin approach requires that the first elements in the set of basis functions only depend on the space but not on time , that is to say, only depend on the space, . Now the numerical solution can be represented in the basis space as

 qh(ξ,η,τ)=L−1∑l=0^qlθl(ξ,η,τ), (3.2)

where is a vector of modes. Similarly, the Hamiltonian can also be represented in the form of (3.2), . The transcription from to will be given in the next subsection. Another simplification of the continuous Galerkin approach is that the solution is continuous with the initial condition at , which means we only have to calculate once at . If the initial condition can be represented in the modal space as

 φnh(ξ,η)=Ls−1∑l=0^wlθl(ξ,η,τ=0), (3.3)

then at , .

Applying the Galerkin approach to (3.1) gives

 ⟨θk,∂θl∂τ⟩^ql+⟨θk,θl⟩^hl=0, (3.4)

where the angled brackets denote the spacetime integration over the reference element, and the Einstein summation convection is used. Eq. (3.4) can be rewritten in the matrix-vector form

 Kτ^q+M^h=0, (3.5)

where and are the time-stiffness matrix and the mass matrix respectively, and the -th elements of them are

 Kτ;k,l=⟨θk,∂θl∂τ⟩,Mk,l=⟨θk,θl⟩. (3.6)

From the previous assumption and simplification, we know that only the last elements of are needed to be determined in the continuous Galerkin predictor step. So we can split into two parts , where is the first components and is the last components. A similar split can be done for , then the mass matrix and the time-stiffness matrix can be written as

 M=[M00M01M10M11],Kτ=[K00τK01τK10τK11τ], (3.7)

where the dimensions of sub-matrices are respectively, and it is similar for the sub-matrices of . In (3.5), only the last components are useful, and they can be written as

 ^q1=−^M^h1−^M0^h0, (3.8)

where . We can obtain from through one iteration using the above equation. In the continuous Galerkin predictor step, times iterations of Eq. (3.8) are adequate for a -th order scheme [1], thus the cost of the iterative part in our scheme is not high. Once the basis functions are determined, the matrices in (3.8) are known, and the whole iterative scheme can be explicitly written down. We are going to describe the implementation details in the next subsection.

### 3.2 Implementation details of 2D third order continuous Galerkin predictor

This subsection will give the implementation details at third order. Other cases can be completed similarly and we will provide some difference of the implementation in other cases in Remark 3.1. Assume that the basis functions in the reference element are orthogonal Legendre polynomials

 P0(ξ)=1,P1(ξ)=ξ,P2(ξ)