Sparse grid central discontinuous Galerkin method for linear hyperbolic systems in high dimensions
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 nonperiodic 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 timedependent linear hyperbolic system with variable coefficients
\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 boxshaped domains. The model (LABEL:eq:model0) arises in many contexts [17], 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 [2]. 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 highdimensional numerical simulations. Sparse grid techniques have been incorporated in collocation methods for highdimensional stochastic differential equations [35, 34, 25, 23], finite element methods [36, 3, 28], finite difference methods [9, 11], finite volume methods [14], and spectral methods [10, 8, 29, 30] for highdimensional 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 [4] is a class of finite element methods using discontinuous approximation space for the numerical solution and the test functions. The RungeKutta DG scheme [5] 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 [12] is well suited for timedependent 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 [12], in this paper, we restrict our attention to smooth solutions of (LABEL:eq:model0). It is known that for nonsmooth solutions, adaptivity should be invoked to capture discontinuity like structures. This can be achieved using the idea in [13] and is left for our future work.
Based on the scheme constructed in [12], 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 nonperiodic 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 initialboundary 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 nonperiodic 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 componentwise form as
\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 nonperiodic 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 [33]. Consider a general interval , we define the th level mesh to be a uniform partition of cells with length and , for any Let
be the usual piecewise polynomials of degree at most on . Then, we have the nested structure
Similar to [33], we can now define the multiwavelet subspace , as the orthogonal complement of in with respect to the inner product on , i.e.,
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 multiindices. For a multiindex , where denotes the set of nonnegative integers, the and norms are defined as
The componentwise arithmetic operations and relational operations are defined as
By making use of the multiindex notation, we denote by the mesh level in a multivariate sense. We define the tensorproduct mesh grid and the corresponding mesh size Based on the grid , we denote by as an elementary cell, and
as the standard tensorproduct 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 tensorproduct construction, the multidimensional increment space can be defined as
Therefore, we have The sparse finite element approximation space we consider, is defined by
This is a subset of , and its number of degrees of freedom scales as [33], 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 [12], we can have estimates for projection operator onto the spaces
To facilitate the discussion, below we introduce some notations about norms and seminorms. 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 seminorm, and to denote the broken Sobolev norm and seminorm that are supported on a general grid . For any set , we define to be the complement set of in For a nonnegative integer and set , we define the seminorm on any domain denoted by
and
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 [12].
Lemma 2.1 ( projection estimate)
Let be projections onto the spaces , respectively, then for , , and which is periodic on , , , we have for
\hb@xt@.01(2.2) 
This lemma shows that the norm and seminorm 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
The semidiscrete 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
\hb@xt@.01(2.3)  
\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 [1]. In 1D, the bases of are denoted by
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 multidimensions are defined by tensor products
where we have used the notation and to denote the multiindex for the bases.
As for temporal schemes, we can use the total variation diminishing RungeKutta (TVDRK) methods [32] to solve the ordinary differential equations for the coefficients resulting from the discretization. To calculate the righthandside of (LABEL:eq:DGformulation_p)(LABEL:eq:DGformulation_d), the fast matrixvector 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 matrixvector product which appears at the righthandside of (LABEL:eq:DGformulation_p)(LABEL:eq:DGformulation_d)
where can be the coefficient of the basis in sparse grid space and are the corresponding onedimensional transform of coefficients from basis to basis in the th dimension in our scheme. Note that we have onedimensional bases in each dimension, and we use to denote the th basis. The bases are ordered according to grid increment. Using Algorithm 1 in [31], we should calculate all the onedimensional transform along each direction associated with a block lower triangular matrix, and then calculate all the onedimensional transforms having a block upper triangular structure. The fast matrixvector product on sparse grid with LU decomposition can be proceeded as follows.

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

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

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 variablecoefficient, we use Gaussian quadrature to compute these terms. Since these integrals are multidimensional integrations, we use the socalled unidirectional principle to separate the integration into multiplication of onedimensional integrals. For example, if is separable,
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 [5], CDG [21], sparse grid DG [12], and the sparse grid CDG schemes. We only consider the twodimensional 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 [5]. 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 twodimensional CDG method is larger than the CFL number for onedimensional CDG method in [21]. 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.
DG  CDG  sparse grid DG  sparse grid CDG  

1  2  3  1  2  3  1  2  3  1  2  3  
0.33  –  –  0.48  –  –  0.66  –  –  0.87  –  –  
0.40  0.20  0.13  0.66  0.36  0.24  0.81  0.41  0.25  1.17  0.65  0.44  
0.46  0.23  0.14  0.90  0.52  0.35  0.92  0.46  0.28  1.58  0.94  0.62 
2.4 Nonperiodic problems
Here, we consider nonperiodic 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 nonperiodic boundary condition on the domain the finite element space on dual mesh with cell size is represented by
\hb@xt@.01(2.5) 
where the mesh is partitioned as
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
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.
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,
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
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 [33]), 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 ( projection estimate onto )
Let be the projection onto the space , then for , , and , , , we have
\hb@xt@.01(2.6) 
Proof. The proof follows same procedure as Appendix A in [12]. 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 onedimensional increment projections. We denote as the standard projection operator from to , and the induced increment projection
and further denote
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
\hb@xt@.01(2.7) 
It suffices to show (LABEL:eq:portho) for which is a dense subset of . In fact, we have
where is the projection onto the full grid space Therefore,
The last term in the first row of the equality above vanishes because In addition, for any ,
Therefore, by properties of the tensor product projections
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 , ,
where the mesh size
The estimation above directly applies for . For , by simple algebra, we have