Sparse grid central discontinuous Galerkin method for linear hyperbolic systems in high dimensions

# Sparse grid central discontinuous Galerkin method for linear hyperbolic systems in high dimensions

Zhanjing Tao Department of Mathematics, Michigan State University, East Lansing, MI 48824 U.S.A. taozhanj@msu.edu    Anqi Chen Department of Mathematics, Michigan State University, East Lansing, MI 48824 U.S.A. chenanq3@msu.edu.    Mengping Zhang School of Mathematical Sciences, University of Science and Technology of China, Hefei, Anhui, 230026 People’s Republic of China. mpzhang@ustc.edu.cn. Research supported by NSFC grant 11471305.    Yingda Cheng Department of Mathematics, Department of Computational Mathematics, Science and Engineering, Michigan State University, East Lansing, MI 48824 U.S.A. ycheng@msu.edu. Research is supported by NSF grants DMS-1453661, DMS-1720023 and the Simons Foundation.
July 1, 2019
###### Abstract

In this paper, we develop sparse grid central discontinuous Galerkin (CDG) scheme for linear hyperbolic systems with variable coefficients in high dimensions. The scheme combines the CDG framework with the sparse grid approach, with the aim of breaking the curse of dimensionality. A new hierarchical representation of piecewise polynomials on the dual mesh is introduced and analyzed, resulting in a sparse finite element space that can be used for non-periodic problems. Theoretical results, such as stability and error estimates are obtained for scalar problems. CFL conditions are studied numerically comparing discontinuous Galerkin (DG), CDG, sparse grid DG and sparse grid CDG methods. Numerical results including scalar linear equations, acoustic and elastic waves are provided.

Key words. central discontinuous Galerkin method; sparse grid; linear hyperbolic system; stability; error estimate

## 1 Introduction

In this paper, we develop sparse grid central discontinuous Galerkin (CDG) method for the following time-dependent linear hyperbolic system with variable coefficients

 ∂u∂t+d∑i=1∂(Ai(t,x)u)∂xi=0,x∈Ω, \hb@xt@.01(1.1)

subject to appropriate initial and boundary conditions. In the expression above, is the spatial dimension of the problem, is the unknown function, are the given smooth variable coefficients. We assume in the paper, but the discussion can be easily generalized to arbitrary box-shaped domains. The model (LABEL:eq:model0) arises in many contexts , such as simulations of acoustic, elastic waves, and Maxwell’s equations in free space. The scheme we develop in this paper can also apply to the case when is defined through another set of equations that can be nonlinearly coupled with , such as the models in kinetic plasma waves and incompressible flows.

Many numerical methods, including finite difference, finite volume, finite element, spectral methods etc., have been developed in the literature for (LABEL:eq:model0) addressing different challenges in various applications. The focus of this paper, is to design a class of conservative numerical schemes, with high computational efficiency, for system (LABEL:eq:model0) when is large. It is well known that any grid based method suffers from the curse of dimensionality . This term refers to the fact that the computational cost and storage requirements scale as for a -dimensional problem, where denotes the mesh size in one coordinate direction, while the approximation accuracy is independent of To overcome this bottleneck, sparse grid methods [36, 3, 7] were introduced to reduce the degrees of freedom for high-dimensional numerical simulations. Sparse grid techniques have been incorporated in collocation methods for high-dimensional stochastic differential equations [35, 34, 25, 23], finite element methods [36, 3, 28], finite difference methods [9, 11], finite volume methods , and spectral methods [10, 8, 29, 30] for high-dimensional PDEs.

In recent years, we initiate a line of research on the development of sparse grid discontinuous Galerkin (DG) methods [33, 12, 13]. The DG method  is a class of finite element methods using discontinuous approximation space for the numerical solution and the test functions. The Runge-Kutta DG scheme  developed in a series of papers for hyperbolic equations became very popular due to its provable convergence, excellent conservation properties and accommodation for adaptivity and parallel implementations. The sparse grid DG method designed in  is well suited for time-dependent transport problems in high dimensions, reducing degrees of freedom of from to , maintaining conservation, with provable convergence rate of in norm when the solution is smooth. Similar to , in this paper, we restrict our attention to smooth solutions of (LABEL:eq:model0). It is known that for non-smooth solutions, adaptivity should be invoked to capture discontinuity like structures. This can be achieved using the idea in  and is left for our future work.

Based on the scheme constructed in , the goal of the present paper is to design and analyze the sparse grid CDG method. The CDG schemes [18, 20, 21] are a class of DG schemes on overlapping cells that combine the idea of the central schemes [24, 16, 19] with the DG weak formulation. Such methods are intrinsically Riemann solver free, therefore no costly flux evaluations are needed in the computation. It is well known that the CDG schemes allow larger CFL numbers than the standard DG methods except for piecewise constant approximations [20, 26]. This compensates the increased cost caused by duplicate representation of the solution on the dual mesh. Motivated by this, we develop sparse grid CDG method that avoids the evaluation of numerical fluxes. We investigate stability, convergence rate and CFL condition of the resulting scheme. A novelty of this work is the design of the scheme for non-periodic problems, where a new hierarchical representation of the solution is presented, which results in a sparse finite element space that can be defined on the dual mesh. projection results are studied for this space, which helps the convergence proof of the schemes for initial-boundary value problems.

The rest of this paper is organized as follows: in Section LABEL:sec:method, we construct the sparse grid CDG formulations for periodic and non-periodic problems, and perform numerical study of the CFL conditions. In Section LABEL:sec:analysis, we prove stability and error estimates for scalar equations. The numerical performance is validated in Section LABEL:sec:numerical by several benchmark tests, including scalar transport equations, acoustic and elastic waves. Conclusions and future work are discussed in Section LABEL:sec:conclu.

## 2 Numerical methods

In this section, we define and discuss the properties of the proposed sparse grid CDG methods. For convenience of notations, we rewrite (LABEL:eq:model0) in a component-wise form as

 ∂ul∂t+∇⋅(Al(t,x)u)=0,l=1,⋯,m,x∈Ω, \hb@xt@.01(2.1)

where denotes a collection of the -th row of each matrix . The problem is solved with given initial value and periodic or Dirichlet type boundary conditions.

We proceed as follows. First, we introduce the scheme for periodic problems. In this setting, the finite element space on the primal and dual mesh can be defined in similar ways. Then, we discuss the implementation details and perform numerical study of the CFL conditions. Finally, we consider the more complicated non-periodic problems, for which a new sparse finite element space will be introduced on the dual mesh.

### 2.1 Periodic problems

To define the sparse finite element space, we first review the hierarchical decomposition of piecewise polynomial space in one dimension . Consider a general interval , we define the -th level mesh to be a uniform partition of cells with length and , for any Let

 Vkn([a,b]):={v:v∈Pk(Ijn),∀j=0,…,2n−1}

be the usual piecewise polynomials of degree at most on . Then, we have the nested structure

 Vk0([a,b])⊂Vk1([a,b])⊂Vk2([a,b])⊂Vk3([a,b])⊂⋯

Similar to , we can now define the multiwavelet subspace , as the orthogonal complement of in with respect to the inner product on , i.e.,

 Vkn−1([a,b])⊕Wkn([a,b])=Vkn([a,b]),Wkn([a,b])⊥Vkn−1([a,b]).

For notational convenience, we let , which is the standard piecewise polynomial space of degree on . This gives the hierarchical decomposition on as .

For a dimensional domain we recall some basic notations about multi-indices. For a multi-index , where denotes the set of nonnegative integers, the and norms are defined as

 |α|1:=∑di=1αi,|α|∞:=max1≤i≤dαi.

The component-wise arithmetic operations and relational operations are defined as

 α⋅β:=(α1β1,…,αdβd),c⋅α:=(cα1,…,cαd),2α:=(2α1,…,2αd),
 α≤β⇔αi≤βi,∀i,α<β⇔α≤β and α≠β.

By making use of the multi-index notation, we denote by the mesh level in a multivariate sense. We define the tensor-product mesh grid and the corresponding mesh size Based on the grid , we denote by as an elementary cell, and

 Vkl([a,b]d):={v:v(x)∈Qk(Ijl),0≤j≤2l−1}=Vkl1,x1([a,b])×⋯×Vkld,xd([a,b])

as the standard tensor-product piecewise polynomial space on this mesh, where denotes the collection of polynomials of degree up to in each dimension on cell . If , the grid and space will be further denoted by and , respectively.

Based on a tensor-product construction, the multi-dimensional increment space can be defined as

 Wkl([a,b]d)=Wkl1,x1([a,b])×⋯×Wkld,xd([a,b]).

Therefore, we have The sparse finite element approximation space we consider, is defined by

 ^VkN([a,b]d):=⨁|l|1≤Nl∈Nd0Wkl([a,b]d).

This is a subset of , and its number of degrees of freedom scales as , which is significantly less than that of with exponential dependence on . This is the key to computational savings in high dimensions.

The standard CDG schemes [18, 20] is characterized by numerical approximations on two sets of overlapping grids: primal and dual meshes. Now, we are ready to incorporate the sparse finite element space defined above into the CDG framework. For the domain under consideration we let be the primal mesh and , which is the periodic extension of restricted to , be the dual mesh. Similarly, we let and to be the periodic extension of restricted to . Here and below, the subscripts and represent the quantities defined on the primal and dual mesh, respectively.

The approximation properties for the sparse finite element space have been established in previous work [33, 12]. By using a lemma in , we can have estimates for projection operator onto the spaces

To facilitate the discussion, below we introduce some notations about norms and semi-norms. Let , on primal or dual mesh , we use to denote the standard broken Sobolev norm, i.e. where is the standard Sobolev norm on (and is used to denote the norm). Similarly, we use to denote the broken Sobolev semi-norm, and to denote the broken Sobolev norm and semi-norm that are supported on a general grid . For any set , we define to be the complement set of in For a non-negative integer and set , we define the semi-norm on any domain denoted by

 |v|Hα,L(Ω′):=∥∥ ∥∥(∂α∂xαi1⋯∂α∂xαir)v∥∥ ∥∥L2(Ω′),

and

 |v|Hq+1(Ω′):=max1≤r≤d⎛⎜⎝maxL⊂{1,2,⋯,d}|L|=r|v|Hq+1,L(Ω′)⎞⎟⎠,

which is the norm for the mixed derivative of of at most degree in each direction. In this paper, we use the notation to represent , where the constant is independent of and the mesh level considered. The following results are obtained from Lemma 3.2 in .

###### Lemma 2.1 (L2 projection estimate)

Let be projections onto the spaces , respectively, then for , , and which is periodic on , , , we have for

 |PGv−v|Hs(ΩN,G)≲{Nd2−N(q+1)|v|Hq+1(Ω)s=0,2−Nq|v|Hq+1(Ω)s=1. \hb@xt@.01(2.2)

This lemma shows that the norm and semi-norm of the projection error scale like and with respect to when the function has bounded mixed derivatives up to enough degrees. This lemma will be used in Theorem LABEL:thm:lterr to establish convergence of the scheme.

Now, we are ready to formulate the sparse grid CDG scheme. Below we review some standard notations about jumps and averages of piecewise functions. With or , let be the collection of all elementary cell , be the union of the interfaces for all the elements in (here we have taken into account the periodic boundary condition when defining ) and be the set of functions defined on . For any and , we define their averages and jumps on the interior edges as follows. Suppose is an interior edge shared by elements and , either on primal or dual mesh, we define the unit normal vectors and on pointing exterior of and , respectively, then

 [q] = q−n−+q+n+, {q}=12(q−+q+), = q−⋅n−+q+⋅n+, {q}=12(q−+q+).

The semi-discrete sparse grid CDG scheme for (LABEL:eq:model), based on the weak formulation introduced in [18, 20], is defined as follows: we find and , such that

 ∫Ω(ulh)tφhdx= 1τmax∫Ω(vlh−ulh)φhdx+∫ΩAl(t,x)vh⋅∇φhdx−∑e∈ΓN,P∫eAl(t,x)vh⋅[φh]ds, \hb@xt@.01(2.3) ∫Ω(vlh)tψhdx= 1τmax∫Ω(ulh−vlh)ψhdx+∫ΩAl(t,x)uh⋅∇ψhdx−∑e∈ΓN,D∫eAl(t,x)uh⋅[ψh]ds, \hb@xt@.01(2.4)

for any and , where and is an upper bound for the time step due to the CFL restriction (see Section LABEL:sec:cfl for detailed discussions).

### 2.2 Discussions on implementations

Here, we briefly discuss some details about the implementation of the scheme. We perform the computation by using orthonormal multiwavelet bases constructed by Alpert . In 1D, the bases of are denoted by

 vjp,l(x),p=1,⋯,k+1,j=0,⋯,2l−1−1

and they satisfy . Figures LABEL:fig:nonp_basis_k01 and LABEL:fig:nonp_basis_k11 provide illustrations of the basis functions for and The bases in in multi-dimensions are defined by tensor products

 vs=vjp,l:=d∏i=1vjipi,li(xi),pi=1,⋯,k+1,ji=0,⋯,max(0,2li−1−1),

where we have used the notation and to denote the multi-index for the bases.

As for temporal schemes, we can use the total variation diminishing Runge-Kutta (TVD-RK) methods  to solve the ordinary differential equations for the coefficients resulting from the discretization. To calculate the right-hand-side of (LABEL:eq:DGformulation_p)-(LABEL:eq:DGformulation_d), the fast matrix-vector product by LU split or LU decomposition algorithms [30, 31, 27] can be applied, by which one can decompose all calculations into one dimensional operations. Below, we briefly describe the LU decomposition algorithm for the calculation of the following matrix-vector product which appears at the right-hand-side of (LABEL:eq:DGformulation_p)-(LABEL:eq:DGformulation_d)

 bj=∑s:|l|1≤Nfst1s1,j1⋯tdsd,jd,

where can be the coefficient of the basis in sparse grid space and are the corresponding one-dimensional transform of coefficients from basis to basis in the -th dimension in our scheme. Note that we have one-dimensional bases in each dimension, and we use to denote the -th basis. The bases are ordered according to grid increment. Using Algorithm 1 in , we should calculate all the one-dimensional transform along each direction associated with a block lower triangular matrix, and then calculate all the one-dimensional transforms having a block upper triangular structure. The fast matrix-vector product on sparse grid with LU decomposition can be proceeded as follows.

1. Calculate (block) LU decomposition for where are the permutation matrices, are lower and upper triangular matrices.

2. Compute the transform with a (block) lower triangular matrix for

3. Compute the transform with a (block) upper triangular matrix for

Note that in step 1, the LU decomposition pivots only from rows or columns in the same mesh level to maintain the hierarchical structure. This pivoting can be successfully done in the sparse grid CDG scheme, but not in the sparse grid DG scheme, for which additional splitting of the flux terms are deemed necessary for variable coefficient case.

For the integrals involving variable-coefficient, we use Gaussian quadrature to compute these terms. Since these integrals are multi-dimensional integrations, we use the so-called unidirectional principle to separate the integration into multiplication of one-dimensional integrals. For example, if is separable,

 ∫Ωϕ(x)=∫[a,b]ϕ1(x1)⋯∫[a,b]ϕd(xd).

When the variable coefficient is separable, we can use unidirectional principle directly. If it is not separable, we can find as the projection of onto the sparse grid finite element space, and then use to compute the integrals.

### 2.3 Discussions on CFL conditions

It is well known that the CDG schemes allow larger CFL numbers than the standard DG methods except for piecewise constant approximations [20, 26]. Here, we perform a numerical study of the CFL conditions of DG , CDG , sparse grid DG , and the sparse grid CDG schemes. We only consider the two-dimensional case solving constant coefficient equation for now. The results are listed in Table LABEL:table:linear_cfl. The CFL number of DG method is obtained from Table 2.2 in . The rest of the table is computed by eigenvalue analysis of the discretization matrix, and by requiring the amplification of the eigenvalues to be bounded by 1 in magnitude. We observe that the sparse grid DG method has CFL number that is about two times the CFL number of the standard DG method. The sparse grid CDG method offers the largest CFL conditions among all four methods. Here, as a side note, we find that the CFL number for two-dimensional CDG method is larger than the CFL number for one-dimensional CDG method in . This table shows that one advantage of the sparse grid CDG method is the ability to take large time steps for time evolution problems. In general, further numerical results suggest that for equation the CFL number for sparse grid DG and sparse grid CDG method will change with the value of the coefficients Results in higher dimensions are yet to be studied. A preliminary calculation shows that for equation the CFL conditions for CDG, sparse grid DG and sparse grid CDG methods in 3D are all higher than those for the 2D case in Table LABEL:table:linear_cfl. The sparse grid CDG method still possesses the largest CFL number among all four methods. Those interesting issues will be investigated in our future work.

### 2.4 Non-periodic problems

Here, we consider non-periodic problems, where equation (LABEL:eq:model0) or (LABEL:eq:model) is supplemented by Dirichlet boundary condition on the inflow edges. In this case, we can no longer use periodicity to define the finite element space on the dual mesh, and a new grid hierarchy needs to be introduced.

Recall that for standard CDG methods with non-periodic boundary condition on the domain the finite element space on dual mesh with cell size is represented by

 Vkn,D={v:v∈Pk(Ijn,D),∀j=0,…,2n}, \hb@xt@.01(2.5)

where the mesh is partitioned as

 I0n,D=[0,hn2],Ijn,D=[(j−12)hn,(j+12)hn],j=1,…,2n−1,I2nn,D=[1−hn2,1],

which consists of cells of size and two cells at the left and right ends of size It is easy to see that this space does not have nested structures, i.e. Therefore, we need a new hierarchy to define the increment polynomial spaces.

For a fixed refined mesh level we define the following grid on level by a collection of cells as

 I0l,N,D=[0,hl−hN2],Ijl,N,D=[jhl−hN2,(j+1)hl−hN2],j=1,…,2l−1,I2ll,N,D=[1−hN2,1],

which consists of cells of size and a cell at the left end of size and a cell at the right end of size This grid structure is naturally nested, and therefore which consists of piecewise polynomials of degree defined on are also nested, and as defined in (LABEL:eq:sd).

Then the definitions of sparse finite element space in Section LABEL:sec:periodic can be naturally extended here. We let , be a complement set of in i.e.

 Vkl−1,N,D⊕Wkl,N,D=Vkl,N,D.

However, we no longer require to be orthogonal to , because such definition will be difficult to implement in practice. Instead, we define to be a span of basis functions that are shifted basis functions of space defined in Section LABEL:sec:periodic, namely,

 Wkl,N,D=Wkl([−hN2,1−hN2])∣∣[0,1],l≥1.

By denoting we have decomposed Illustration of basis functions by such definitions for and can be found in Figures LABEL:fig:nonp_basis_k02 and LABEL:fig:nonp_basis_k12. The dimension of is while the dimensions of are

Finally, the sparse finite element space on the dual mesh of domain is defined as

 ^~VkN,D:=⨁|l|1≤Nl∈Nd0Wkl,N,D,

where This is a subset of the full grid space , and its number of degrees of freedom scales as (the proof is similar to Lemma 2.3 in ), which is larger than that of , but still significantly less than that of with exponential dependence on .

We will now investigate the approximation property of the space We can obtain the following result, which essentially states that the projection onto this newly constructed space has the same order of accuracy as in Lemma LABEL:lem:ltproj.

###### Lemma 2.2 (L2 projection estimate onto ^~VkN,D )

Let be the projection onto the space , then for , , and , , , we have

 |~PDv−v|Hs(ΩN,D)≲{Nd2−N(q+1)|v|Hq+1(Ω)s=0,2−Nq|v|Hq+1(Ω)s=1. \hb@xt@.01(2.6)

Proof. The proof follows same procedure as Appendix A in . We will mainly highlight the difference in the proof (see Steps 1 and 2 below). The main difference lies in the fact that all the hierarchical spaces (and associated projections) have dependence not only on but also on the finest mesh level

Step 1: Decomposition of into tensor products of one-dimensional increment projections. We denote as the standard projection operator from to , and the induced increment projection

 Qkl,N,D:={Pkl,N,D−Pkl−1,N,D,if l=1,…N,Pk0,N,D,if l=0,

and further denote

 ~PkN,D:=∑|l|1≤Nl∈Nd0Qkl1,N,D,x1⊗⋯⊗Qkld,N,D,xd,

where the last subindex of indicates that the increment operator is defined in -direction. We can verify that In fact, for any , it’s clear that . Therefore, we only need

 ∫Ω(~PkN,Dv−v)wdx=0,∀w∈^~VkN,D. \hb@xt@.01(2.7)

It suffices to show (LABEL:eq:portho) for which is a dense subset of . In fact, we have

 v=PkN,Dv+v−PkN,Dv,

where is the projection onto the full grid space Therefore,

 ∫Ω(~PkN,Dv−v)wdx =∫Ω(~PkN,Dv−PkN,Dv)wdx+∫Ω(v−PkN,Dv)wdx =−∫Ω(∑|l|∞≤N,|l|1>Nl∈Nd0Qkl1,N,D,x1⊗⋯⊗Qkld,N,D,xdv)wdx.

The last term in the first row of the equality above vanishes because In addition, for any ,

 ∫[0,1]Qkl,N,Dϕφdx=∫[0,1](I−Pkl−1,N,D)ϕφdx−∫[0,1](I−Pkl,N,D)ϕφdx=0,

Therefore, by properties of the tensor product projections

 ∫Ω(~PkN,Dv−v)wdx=0,∀w∈^~VkN,D,

and the proof for is complete.

Step 2: Estimation of the increment projections. For a function , we have the convergence property of the projection as follows: for any integer with , ,

 |Pkl,N,Dv−v|Hs(Ijl,N,D) ≤ck,s,q(hjl,N)(q+1−s)|v|Hq+1(Ijl,N,D),j=1,⋯,2l−1,

where the mesh size

The estimation above directly applies for . For , by simple algebra, we have

 |Qkl,N,Dv|Hs(Ijl,N,D) ≤~ck,s,q2−l(q+1−s)|v|Hq+1(I⌊j/2⌋l−1,N,D),j=2,⋯,2l−1, |Qkl,N,Dv|Hs(Ijl,N,D) ≤ck,s,q(hl)(q+1−s)|v|Hq+1(Ijl,N,D)+ck,s,q(hl−1−hN/2)(q+1−s)|v|Hq+1(I0l−1,N,D),j=0,1, <~ck,s,q2−l(q+1−s)|v|Hq+1(I0l−1,N,D), |Q