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

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 $[a,b]$, we define the $n$-th level mesh $\Omega_{n}([a,b])$ to be a uniform partition of $2^{n}$ cells with length $h_{n}=2^{-n}(b-a)$ and $I_{n}^{j}=[a+jh_{n},a+(j+1)h_{n}]$, $j=0,\ldots,2^{n}-1,$ for any $n\geq 0.$ Let

be the usual piecewise polynomials of degree at most $k$ on $\Omega_{n}$. Then, we have the nested structure

Similar to [33], we can now define the multiwavelet subspace $W_{n}^{k}([a,b])$, $n=1,2,\ldots$ as the orthogonal complement of $V^{k}_{n-1}([a,b])$ in $V^{k}_{n}([a,b])$ with respect to the $L^{2}$ inner product on $[a,b]$, i.e.,

For notational convenience, we let $W_{0}^{k}([a,b]):=V_{0}^{k}([a,b])$, which is the standard piecewise polynomial space of degree $k$ on $[a,b]$. This gives the hierarchical decomposition $V_{n}^{k}([a,b])$ on $\Omega_{n}$ as $V_{n}^{k}([a,b])=\bigoplus_{0\leq l\leq n}W_{l}^{k}([a,b])$.

For a $d$ dimensional domain $[a,b]^{d},$ we recall some basic notations about multi-indices. For a multi-index $\mathbf{\alpha}=(\alpha_{1},\cdots,\alpha_{d})\in\mathbb{N}_{0}^{d}$, where $\mathbb{N}_{0}$ denotes the set of nonnegative integers, the $l^{1}$ and $l^{\infty}$ norms are defined as

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

By making use of the multi-index notation, we denote by $\mathbf{l}=(l_{1},\cdots,l_{d})\in\mathbb{N}_{0}^{d}$ the mesh level in a multivariate sense. We define the tensor-product mesh grid $\Omega_{\mathbf{l}}([a,b]^{d})=\Omega_{l_{1}}([a,b])\otimes\cdots\otimes\Omega_{l_{d}}([a,b])$ and the corresponding mesh size $h_{\mathbf{l}}=(h_{l_{1}},\cdots,h_{l_{d}}).$ Based on the grid $\Omega_{\mathbf{l}}$, we denote by $I_{\mathbf{l}}^{\mathbf{j}}=\{\mathbf{x}:x_{i}\in I_{l_{i}}^{j_{i}},i=1,\cdots,d\}$ as an elementary cell, and

as the standard tensor-product piecewise polynomial space on this mesh, where $Q^{k}(I^{\mathbf{j}}_{\mathbf{l}})$ denotes the collection of polynomials of degree up to $k$ in each dimension on cell $I^{\mathbf{j}}_{\mathbf{l}}$. If $\mathbf{l}=(N,\cdots,N)$, the grid and space will be further denoted by $\Omega_{N}([a,b]^{d})$ and ${\bf V}_{N}^{k}([a,b]^{d})$, respectively.

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

Therefore, we have ${\bf V}_{N}^{k}([a,b]^{d})=\bigoplus_{\begin{subarray}{c}|\mathbf{l}|_{\infty}\leq N\\ \mathbf{l}\in\mathbb{N}_{0}^{d}\end{subarray}}\mathbf{W}_{\mathbf{l}}^{k}([a,b]^{d}).$ The sparse finite element approximation space we consider, is defined by

This is a subset of ${\bf V}_{N}^{k}([a,b]^{d})$, and its number of degrees of freedom scales as $O((k+1)^{d}2^{N}N^{d-1})$ [33], which is significantly less than that of ${\bf V}_{N}^{k}([a,b]^{d})$ with exponential dependence on $Nd$. 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 $\Omega=[0,1]^{d},$ we let $\Omega_{N,P}:=\Omega_{N}([0,1]^{d})$ be the primal mesh and $\Omega_{N,D}$, which is the periodic extension of $\Omega_{N}([-h_{N}/2,1-h_{N}/2]^{d})$ restricted to $[0,1]^{d}$, be the dual mesh. Similarly, we let $\hat{\bf V}_{N,P}^{k}:=\hat{\bf V}_{N}^{k}([0,1]^{d})$ and $\hat{\bf V}_{N,D}^{k}$ to be the periodic extension of $\hat{\bf V}_{N}^{k}([-h_{N}/2,1-h_{N}/2]^{d})$ restricted to $[0,1]^{d}$. Here and below, the subscripts $P$ and $D$ 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 $L^{2}$ projection operator onto the spaces $\hat{\bf V}_{N,P}^{k},\hat{\bf V}_{N,D}^{k}.$

To facilitate the discussion, below we introduce some notations about norms and semi-norms. Let $G=P,D$, on primal or dual mesh $\Omega_{N,G}$, we use $\|\cdot\|_{H^{s}(\Omega_{N,G})}$ to denote the standard broken Sobolev norm, i.e. $\|v\|^{2}_{H^{s}(\Omega_{N,G})}=\sum_{\mathbf{0}\leq\mathbf{j}\leq 2^{\mathbf{N}}-\mathbf{1}}\|v\|^{2}_{H^{s}(I_{N,G}^{\mathbf{j}})},$ where $\|v\|_{H^{s}(I_{N,G}^{\mathbf{j}})}$ is the standard Sobolev norm on $I_{N,G}^{\mathbf{j}},$ (and $s=0$ is used to denote the $L^{2}$ norm). Similarly, we use $|\cdot|_{H^{s}(\Omega_{N,G})}$ to denote the broken Sobolev semi-norm, and $\|\cdot\|_{H^{s}(\Omega_{\mathbf{l},G})},|\cdot|_{H^{s}(\Omega_{\mathbf{l},G})}$ to denote the broken Sobolev norm and semi-norm that are supported on a general grid $\Omega_{\mathbf{l},G}$. For any set $L=\{i_{1},\ldots i_{r}\}\subset\{1,\ldots d\}$, we define $L^{c}$ to be the complement set of $L$ in $\{1,\ldots d\}.$ For a non-negative integer $\alpha$ and set $L$, we define the semi-norm on any domain denoted by $\Omega^{\prime}$

and

which is the norm for the mixed derivative of $v$ of at most degree $q+1$ in each direction. In this paper, we use the notation $A\lesssim B$ to represent $A\leq\text{constant}\times B$, where the constant is independent of $N$ and the mesh level considered. The following results are obtained from Lemma 3.2 in [12].

This lemma shows that the $L^{2}$ norm and $H^{1}$ semi-norm of the projection error scale like $O(N^{d}2^{-N(k+1)})$ and $O(2^{-Nk})$ with respect to $N$ when the function $v$ has bounded mixed derivatives up to enough degrees. This lemma will be used in Theorem 3.2 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 $G=P$ or $D$, let $T_{h,G}$ be the collection of all elementary cell $I^{\mathbf{j}}_{N,G}$, $\Gamma_{N,G}:=\bigcup_{T\in\Omega_{N,G}}\partial T$ be the union of the interfaces for all the elements in $\Omega_{N,G}$ (here we have taken into account the periodic boundary condition when defining $\Gamma_{N,G}$) and $S(\Gamma_{G}):=\Pi_{T\in\Omega_{N,G}}L^{2}(\partial T)$ be the set of $L^{2}$ functions defined on $\Gamma_{N,G}$. For any $q\in S(\Gamma_{N,G})$ and $\mathbf{q}\in[S(\Gamma_{N,G})]^{d}$, we define their averages $\{q\},\{\mathbf{q}\}$ and jumps $[q],[\mathbf{q}]$ on the interior edges as follows. Suppose $e$ is an interior edge shared by elements $T_{+}$ and $T_{-}$, either on primal or dual mesh, we define the unit normal vectors $\bm{n}^{+}$ and $\bm{n}^{-}$ on $e$ pointing exterior of $T_{+}$ and $T_{-}$, respectively, then

The semi-discrete sparse grid CDG scheme for (2), based on the weak formulation introduced in [18, 20], is defined as follows: we find $u^{l}_{h}\in\hat{{\bf V}}_{N,P}^{k}$ and $v^{l}_{h}\in\hat{{\bf V}}_{N,D}^{k}$, such that $\forall\,l=1,\cdots,m$

for any $\varphi_{h}\in\hat{{\bf V}}_{N,P}^{k}$ and $\psi_{h}\in\hat{{\bf V}}_{N,D}^{k}$, where ${\bf u}_{h}=(u_{h}^{1},\cdots,u_{h}^{m}),{\bf v}_{h}=(v_{h}^{1},\cdots,v_{h}^{m})$ and $\tau_{\max}$ is an upper bound for the time step due to the CFL restriction (see Section 2.3 for detailed discussions).

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 $W_{l}^{k}([0,1])$ are denoted by

and they satisfy $\int_{a}^{b}v_{p,l}^{j}(x)v_{p^{\prime},l^{\prime}}^{j^{\prime}}(x)dx=\delta_{pp^{\prime}}\delta_{jj^{\prime}}\delta_{ll^{\prime}}$. Figures 1(a) and 2(a) provide illustrations of the basis functions for $k=0,1$ and $l=0,1,2.$ The bases in $W_{\mathbf{l}}^{k}$ in multi-dimensions are defined by tensor products

where we have used the notation ${\bf s}=(\mathbf{l},\mathbf{j},\mathbf{p})$ and $s_{i}=(l_{i},j_{i},p_{i})$ to denote the multi-index for the bases.

As for temporal schemes, we can use the total variation diminishing Runge-Kutta (TVD-RK) methods [32] to solve the ordinary differential equations for the coefficients resulting from the discretization. To calculate the right-hand-side of (4)-(5), 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 (4)-(5)

where $f_{\bf s}$ can be the coefficient of the basis in sparse grid space and $t_{s_{i},j_{i}}^{i},\,i=1,\cdots,d,$ are the corresponding one-dimensional transform of coefficients from basis $v_{s_{i}}$ to basis $v_{j_{i}}$ in the $i$-th dimension in our scheme. Note that we have $n=2^{N}(k+1)$ one-dimensional bases in each dimension, and we use $v_{s_{i}}$ to denote the $s_{i}$-th basis. The bases are ordered according to grid increment. Using Algorithm 1 in [31], 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 $f_{\bf s}\rightarrow b_{\bf j}$ on sparse grid with LU decomposition can be proceeded as follows.

1. 1.

   Calculate (block) LU decomposition $t_{s,j}^{i}=\sum_{m=1}^{n}(Pl)^{i}_{s,m}(uQ)^{i}_{m,j},\,s,j=1,\cdots,n,$ for $i=1,\cdots,d,$ where $P^{i},Q^{i}$ are the permutation matrices, $l^{i},u^{i}$ are lower and upper triangular matrices.

2. 2.

   Compute the transform with a (block) lower triangular matrix for $i=1,\cdots,d,$

   $b_{s_{1},\cdots,s_{i-1},s_{i}^{{}^{\prime}},s_{i+1},\cdots,s_{d}}\leftarrow\sum_{s_{i}:l_{1}+\cdots+l_{d}\leq N}f_{\bf s}(Pl)_{s_{i},s_{i}^{{}^{\prime}}}^{i}.$

3. 3.

   Compute the transform with a (block) upper triangular matrix for $i=1,\cdots,d,$

   $b_{\bf s}\leftarrow\sum_{s_{i}^{{}^{\prime}}:l_{1}+\cdots+l_{i-1}+l_{i}^{{}^{\prime}}+l_{i+1}+\cdots+l_{d}\leq N}b_{s_{1},\cdots,s_{i-1},s_{i}^{{}^{\prime}},s_{i+1},\cdots,s_{d}}(uQ)_{s_{i}^{{}^{\prime}},s_{i}}^{i}.$

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 $\phi(x)=\phi_{1}(x_{1})\cdots\phi_{d}(x_{d})$ is separable,

When the variable coefficient $A_{i}(t,x)$ is separable, we can use unidirectional principle directly. If it is not separable, we can find $A_{i}^{h}(t,x)$ as the $L^{2}$ projection of $A_{i}(t,x)$ onto the sparse grid finite element space, and then use $A_{i}^{h}(t,x)$ to compute the integrals.

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 two-dimensional case solving constant coefficient equation $u_{t}+u_{x_{1}}+u_{x_{2}}=0$ for now. The results are listed in Table 1. 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 two-dimensional CDG method is larger than the CFL number for one-dimensional 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 $u_{t}+c_{1}u_{x_{1}}+c_{2}u_{x_{2}}=0,$ the CFL number for sparse grid DG and sparse grid CDG method will change with the value of the coefficients $c_{1},c_{2}.$ Results in higher dimensions are yet to be studied. A preliminary calculation shows that for equation $u_{t}+u_{x_{1}}+u_{x_{2}}+u_{x_{3}}=0$ 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 1. 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.

Here, we consider non-periodic problems, where equation (1) or (2) 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 $[0,1],$ the finite element space on dual mesh with cell size $h_{n}=1/2^{n}$ is represented by

where the mesh is partitioned as

which consists of $2^{n}-1$ cells of size $h_{n},$ and two cells at the left and right ends of size $h_{n}/2.$ It is easy to see that this space does not have nested structures, i.e. $V_{n-1,D}^{k}\not\subset V_{n,D}^{k}.$ Therefore, we need a new hierarchy to define the increment polynomial spaces.

For a fixed refined mesh level $N,$ we define the following grid $\Omega_{l,N,D}$ on level $l,\,l=0\ldots N,$ by a collection of cells as

which consists of $2^{l}-1$ cells of size $h_{l},$ and a cell at the left end of size $h_{l}-\frac{h_{N}}{2},$ and a cell at the right end of size $\frac{h_{N}}{2}.$ This grid structure is naturally nested, and therefore $V^{k}_{l,N,D}$ which consists of piecewise polynomials of degree $k$ defined on $\Omega_{l,N,D}$ are also nested, and $V^{k}_{N,N,D}=V_{N,D}^{k}$ as defined in (6).

Then the definitions of sparse finite element space in Section 2.1 can be naturally extended here. We let $W_{l,N,D}^{k}$, $l=1,\ldots N$ be a complement set of $V_{l-1,N,D}^{k}$ in $V_{l,N,D}^{k},$ i.e.

However, we no longer require $W_{l,N,D}^{k}$ to be $L^{2}$ orthogonal to $V_{l-1,N,D}^{k}$, because such definition will be difficult to implement in practice. Instead, we define $W_{l,N,D}^{k}$ to be a span of basis functions that are shifted basis functions of $W_{l}^{k}$ space defined in Section 2.1, namely,

By denoting $W_{0,N,D}^{k}=V_{0,N,D}^{k},$ we have decomposed $V^{k}_{N,D}=\bigoplus_{0\leq l\leq N}W_{l,N,D}^{k}.$ Illustration of basis functions by such definitions for $k=0,1$ and $l=0,1,2$ can be found in Figures 1(b) and 2(b). The dimension of $W_{0,N,D}^{k}$ is $2(k+1),$ while the dimensions of $W_{l,N,D}^{k},\,l=1,\ldots N$ are $2^{l-1}(k+1).$

Finally, the sparse finite element space on the dual mesh of domain $[0,1]^{d}$ is defined as

where $\mathbf{W}_{\mathbf{l},N,D}^{k}=W_{l_{1},N,D,x_{1}}^{k}\times\cdots\times W_{l_{d},N,D,x_{d}}^{k}.$ This is a subset of the full grid space ${\bf V}_{N,D}^{k}=\bigoplus_{\begin{subarray}{c}|\mathbf{l}|_{\infty}\leq N\\ \mathbf{l}\in\mathbb{N}_{0}^{d}\end{subarray}}\mathbf{W}_{\mathbf{l},N,D}^{k}$, and its number of degrees of freedom scales as $O(2^{d-1}(k+1)^{d}2^{N}N^{d-1})$ (the proof is similar to Lemma 2.3 in [33]), which is larger than that of $\hat{\bf V}_{N,P}^{k}$, but still significantly less than that of ${\bf V}_{N,D}^{k}$ with exponential dependence on $Nd$.

We will now investigate the approximation property of the space $\hat{\tilde{{\bf V}}}_{N,D}^{k}.$ We can obtain the following result, which essentially states that the $L^{2}$ projection onto this newly constructed space has the same order of accuracy as $\mathbf{P}_{P},\mathbf{P}_{D}$ in Lemma 2.1.