A numerical study of superconvergence of the discontinuous Galerkin method by patch reconstruction

1.   Introduction

The discontinuous Galerkin method by patch reconstruction (PRDG) was firstly introduced in [12] by Li et al. to solve elliptic equations. The method has been successfully applied into the biharmonic problem [11], the Stokes problem [13,14], the eigenvalue problem [15], the linear elasticity problem [16], the convection-diffusion problem [19] and the sequential least square method [17]. This method was originally motivated by the patch reconstruction technique in the finite volume scheme for the hydrodynamic solver [20]. Based on the least square patch reconstruction, a novel piecewise polynomial discontinuous finite element space is constructed. The space is a subspace of the general discontinuous Galerkin(DG) space, which enables the method to enjoy the well-developed theories and schemes for the DG method. The propose of this article is numerically investigate the superconvergence property of the PRDG method.

Superconvergence properties of the finite element method are classical topics which have been well studied [10,3,21,18]. The properties can be employed to design the a posterior estimator for $h$-refinement [1] and the post-processing strategies [26]. In the last decades, there are numerous works focusing on the superconvergence behavior of the DG method [24,6,9,8,25,4]. In particular, Cockburn et al. [7] studied the superconvergence of the LDG method for elliptic equations on the Cartesian grid. Castillo analyzed the superconvergence properties for the DG method with conservative numerical flux in [5]. The weak Galerkin(WG) method was introduced by Wang and Ye in [22], and also observed a superconvergence properties numerically. In [23], Wang et al. analyzed the superconvergence for the polynomial preserving recovered gradient of the WG method.

In this article, we perform a numerical study of the superconvergence of the PRDG method and present the detailed algorithms. The PRDG method only takes one degree of freedom (DOF) per element and achieves arbitrary order with identical numbers of unknowns. The method utilizes the DOFs with very high efficiency. In some particular cases, the method even can attain higher accuracy than FEMs with the same quantity of unknowns. The superconvergence property of the PRDG method is also conducive to efficiency, which means the method can attain a desirable digit precision with only a few DOFs. Another advantage of the PRDG method is that it has great flexibility on meshes, e.g. the general polygonal mesh is allowed. However, we only consider the uniform meshes in this paper for the purpose of superconvergence investigation.

The article is organized as follows. The approximation space and the details of the algorithm are introduced in Section 2. We first describe the principle to choose the element patches and the sampling nodes, and then elaborate the process to calculate the global stiff-matrix. Section 3 states the approximation properties of such space and the superconvergence properties of the proposed method. Numerical experiments are presented in Section 4 to demonstrate the properties of the resulting linear system and to verify the theoretical results in all 3 dimensions.

2.   Numerical implementation of the discontinuous Galerkin method by patch reconstruction

We consider an open bounded Lipschitz domain ${\Omega}$ in $\mathbb{{R}}^D$, $D = 1,2,3$, such as a convex polygonal domain in $\mathbb{{R}}^D$. Let $\mathcal{T}_h$ be a uniform partition of domain ${\Omega}$. Let $h_K$ and $\left\lvert{K}\right\rvert$ denote its diameter and area/volume, respectively. With the uniform partition, $h{:} = h_K$. We focus our discussion on uniform meshes for the purpose of the superconvergence investigation although an arbitrary polygonal/polyhedron mesh is allowed in general. The shape regularity constraints will be claimed in the next section, and we just focus on the numerical implementation in this section.

2.1.   Reconstruction and interpolation

With qualified partitions $\mathcal{T}_h$, we first define a finite element space $V_h$ and an interpolation operator $\mathcal{{R}}$ on the mesh at a reasonable cost. If the high order polynomial approximation is needed, a common approach is to define numerous DOFs locally on each element. To avoid abusing the DOFs, we define only one DOF per element. It is denoted with $x_K$ for the sampling node. Let $U_h$ be the piecewise constant space associated with $\mathcal{T}_h$, i.e.,

Now we construct the finite element space $V_h$ by the space $U_h$. More precisely, the finite element space with piecewise polynomials $V_h$ can be regarded as $U_h$ embedded by the reconstruction operator $\mathcal{{R}}$, which can be expressed as follows,

Since there is only one DOF per element and it requires more DOFs to construct higher order polynomial, the patch reconstruction technique is used. The high order polynomial is constructed from the element patch $S(K)$. $S(K)$ is a subset of $\mathcal{T}_h$ that includes $K$ itself and neighboring elements around $K$. The sampling node $x_{K}$ is assigned for each element $K$, which is a point located inside the element. The details of how to set up the sampling nodes and the element patch will be discussed later. $\mathcal{{I}}_K$ denote the set of sampling nodes corresponding to $S(K)$. Let $\# \mathcal{{I}}_K$ denote the number of elements belonging to $\mathcal{{I}}_K$, and $\#S(K)$ denote the number of elements belonging to $S(K)$.

Now we introduce how to obtain the local approximation polynomial with the sampling nodes and the element patch. For $\forall v\in U_h$ and $\forall K\in \mathcal{T}_h$, a local polynomial $\mathcal{{R}}_{K} v$ of degree $m$ can be found by seeking the best local least-square approximation.

$\mathcal{{R}}_K v$ gives an interpolation polynomial on the patch $S(K)$. We limit the interpolation on element $K$. The global interpolation operator $\mathcal{{R}}$ is given by

Above is the general process to construct the finite element space $V_h$, and some practical issues should be addressed here. Firstly, the barycenters of element $K$ are preferred as the sampling nodes $x_K$ when we investigate the superconvergence of the PRDG method, although they can be chosen inside $K$. Secondly, the number $\# S(K)$ should be larger than the DOFs that the polynomial of degree $m$ needs. The number is

where $D$ is the dimension of space. A threshold value $c_0$ is assigned to $\# S(K)$, which is reached recursively by adding the nearest Von Neumann neighbors to $S(K)$. The recursive process is terminated when the size of $S(K)$ meets the requirement. Figure 2.1 demonstrate the appropriate patch with proper regularity and the inappropriate patch with an isolated element. The regularity of the element patch might influent the approximation convergence result. And this is the reason why we use the above principle to choose the element patch.

The least square problem (2.1) can be easily solved by the generalized inverse of matrix $M = (W^TW)^{-1}W^T$, where $W$ is the Vandermonde-type matrix. The details of examples in 1D can be found in [11], 2D in [13] and 3D in [15]. Here we demonstrate an one dimensional example to explain the relation between the matrix $M_K$ and the basis function. Assume the element $K_1$ has the element patch $S(K_1) = \{K_1,K_2,K_3\}$, and the sampling nodes are $x_{K_1},x_{K_2},x_{K_3}$, the linear approximation is obtained by solving the local least-square problem:

The above problem can be directly solved by the generalized inverse,

The Vandermonde-type matrix $W$ is specified as,

$q$ is a vector

Thus the matrix $M_{K_1} = (W^TW)^{-1}W^T$ can be expressed explicitly,

The linear approximation polynomial $\mathcal R_{K_1}\; v$ is

where

We reorganize that

Obvious, the basis functions are as follows,

The linear polynomial can be rewritten as

More details can be found in [11] for 1D case.

The matrix $M_{K}$ stores the information of basis functions whose values are not zero in element $K$. The matrix $M_{K}$ only depends on the sampling nodes set $\mathcal{{I}}_K$ and it is irrelevant to the vector $v\in U_h$. $M_{K}$ is calculated offline and stored in memory, which is distinct from the FEM and DG methods. The basis function $\lambda_{K}$ is the interpolation function of the characteristic function corresponding to $K$. And the basis functions are not explicitly used in the calculation. The interpolation polynomial with order $m$ of the vector $v\in U_h$ is specified as

where $L$ is the basis vector of a polynomial with order $m$. $L = [1,x,x^2,...,x^m]$ for 1D, $L = [1,x,y,x^2,xy,y^2,...,y^m]$ for 2D, and $L = [1,x,y,z,x^2,xy,y^2,xz,yz,z^2,..., z^m]$ for 3D, respectively.

2.2.   Computation of the stiffness matrix

We consider the following Poisson problem and solve it with the PRDG method,

Let $\mathcal{{E}}_h$ denote the collection of all edges of $\mathcal{T}_h$, $\mathcal{{E}}_h^{0}$ denote the collection of all the interior edges and $\mathcal{{E}}_h^{b}$ denote the collection of all boundary edges. The variational form of the Poisson problem is to seek $u_h\in U_h$ such that

The IPDG scheme [2] can be directly applied, where the bilinear term $a_h$ and linear operator $(f,\mathcal{{R}}v_h)_h$ are defined as: for any $v_h,w_h\in U_h$

and

where $\eta$ is a positive constant. The average $\left\{\,{v}\,\right\}$ and the jump operator $[ [ v ] ]$ is defined as follows. Assume $e$ is a common edge shared by elements $K_1$ and $K_2$, and let $n_1$ and $n_2$ be the outward unit normal at $e$ of $K_1$ and $K_2$, respectively. Given $v_i{:} = \left. v \right|_{\partial K_i}$, we define

For a vector-valued function $\varphi$, we define $\varphi_1$ and $\varphi_2$ similarly,

For $e \in\mathcal{{E}}_h^{b}$, let

The linear algebra system is obtained from the weak form (2.3),

where $A$ is the global stiff matrix whose size is $N\times N$, $u_h\in U_h$ is a $N\times 1$ vector, and the right hand side (RHS) $F$ is also a $N\times 1$ vector. $N$ is the number of elements in $\mathcal{T}_h$.

We now illuminate the computation of the right hand side $F$. The $i-$th component of $F$ is the numerical integration of the inner product between $f$ and $\lambda_{K_i}$, where $\lambda_{K_i}$ is the basis function corresponding to the element $K_i$, i.e.

As we mentioned, the basis function $\lambda_{K_i}$ is not explicitly used in the computation. Instead, we calculate the local right hand side $\psi_{K_i}$ element by element, and then assemble the local RHS $\psi_{K_i}$ to the global RHS $F$.

The corresponding relation between $\psi_{K_i}$ and $F$ is determined by $s_{K_i}$. $s_{K_i}$ is the element index of $S(K_i)$.

Similarly, the elements in the global stiff matrix $A$ is computed from $a_h(\lambda_{K_i}, \lambda_{K_j})$. It is not explicitly calculated. Instead, the local element stiff matrix $\kappa_{K_{i}}$ and local trace stiff matrix $\kappa_{e}$ are calculated per element in $\mathcal{T}_h$ and per edge in $\mathcal{{E}}_h$, respectively. Assembling the locate stiff matrixes together gives $A$.

The size of the local element stiff matrix is $\# S(K_{i}) \times \# S(K_{i})$, and $\kappa_{K_{i}}$ is computed as follows

where ${\nabla}_x L$ is the gradient operator on $L$, ${\nabla}_x L = [0,1,0,0,2x,y,0,z,0,0,..., mx^{m-1}]$. The vector ${\nabla}_y L$ and ${\nabla}_z L$ are given in the same way. Because the $M_{K_{i}}$ is calculated offline, the Jacobian matrix and the affine transformation are no longer needed. Numerical integration is conducted on the real element $K_{i}$ instead of the reference element. Furthermore, the corresponding matrix is constituted by $s_{K_{i}}^{T}s_{K_{i}}$.

Consider an interior edge $e\in \mathcal{{E}}_h^{0}$, which is shared by elements $K_{i}$ and $K_{j}$. The local trace stiff matrix $\kappa_{e}$ is calculated on edge $e$ which is relevant with two element patches $S(K_{i})$ and $S(K_{j})$, and its size is $\left (\# S(K_{i})+\# S(K_{j}) \right) \times \left (\# S(K_{i})+\# S(K_{j}) \right)$. In practice, the matrix $\kappa_{e}$ is decomposed to four submatrix $\kappa_{e_{1}}$, $\kappa_{e_{2}}$, $\kappa_{e_{3}}$, $\kappa_{e_{4}}$ for simplicity, which reads

The sizes of four submartrices are $\# S(K_{i}) \times \# S(K_{i})$, $\# S(K_{i}) \times \# S(K_{j})$, $\# S(K_{j}) \times \# S(K_{i})$ and $\# S(K_{j}) \times \# S(K_{j})$, respectively. Numerical integration of $\kappa_{e_{1}}$ is calculated by

where $n_{i} = (n_{x_{i}}, n_{y_{i}} ,n_{z_{i}} )$ is the unit outer norm of $K_i$ on $e$. The other submatrix and the corresponding relation are computed analogously.

We end this section by some comments about the linear system (2.3). The number of unknowns in $u_h$ always equals to the number of elements $N$. The scale of stiff matrix $A$ always maintains $N\times N$ regardless of the approximation order $m$ which may vary. The sparsity pattern of $A$ may vary when the size of element patch $S(K)$ varies.

3.   Superconvergence results for elliptic problems

In this section, we discuss the convergence property of the PRDG method for elliptic problems. The mesh partition $\mathcal{T}_h$ is uniform in all dimensions. So it satisfies the shape regularity conditions naturally [12]. For the element patch $S(K)$, we define $d_K: = \text{diam} S(K)$ and $d = \max_{K\in \mathcal{T}_h}d_K$. Moreover, we assume the following.

Assumption 1. For $\forall K\in \mathcal{T}_h$, there exist constants $C$ and $c$ that are independent of $K$, such that $B_c\subset S(K)\subset B_C$ with $C\ge 2c$, and $S(K)$ is star-shaped with respect to $B_c$, where $B_c$ is a disk with radius $c$.

The Assumption 1 is the geometric constraint for the element patch, which is also the reason why we must obey the principle to choose the candidate of element patch. Next, we claim the assumption on the sampling node set.

Assumption 2. For any $K\in \mathcal{T}_h$ and $p\in\mathbb{{P}}_m(S(K))$,

This assumption leads to the uniqueness of the least squares problem (2.1), which implies that the number $\# \mathcal{{I}}_K$ must be large enough. However, the size of the element patch is a subtle issue. When the element patch is too large, the diameter $d$ of $S(K)$ increases and might change the numerical performance. There is a quantitative estimate of this assumption,

with

The constant $\Lambda(m, \mathcal{{I}}(K))$ is analogous to the Lebesgue constant in approximation theory. The uniform upper bound of $\Lambda(m, \mathcal{{I}}(K))$ can be found if the element patch satisfies some constraints, we refer to [12] for details. Now, we are ready to state the approximation property of the local reconstruction operator,

Lemma 3.1. [12,Lemma 3] There exists a unique solution to (2.1) while Assumption 2 holds. Furthermore, $\mathcal{{R}}_K$ satisfies

For any $K\in \mathcal{T}_h$ and $g\in C^0(S(K))$, the stability property holds

and the quasi-optimal approximation property is valid

where $\Lambda_m{:} = \max_{K\in \mathcal{T}_h} \{1+\Lambda(m,\mathcal{{I}}(K))\sqrt{\#</italic> <italic>\mathcal{{I}}(K)}\}$.

The nearly optimal approximation property naturally exists.

Lemma 3.2. [12,Lemma 4] If Assumption 2 holds, then there exists $C$ for $u\in C^0(\Omega) \cup H^{m+1}(\Omega)$ such that

The approximation estimates of the global reconstruction operator are as follows,

Lemma 3.3. [12,Equation (3.4)] For $u\in H^{m+1}(\Omega)$, together with the Agmon inequality and the local approximation estimates (3.6) and (3.7), there exist a positive constant $C$, such that

3.1.   General convergence property

We now introduce the convergence property of the PRDG method for elliptic problems. The coercivity and the boundedness of $a_h$ are obvious. For sufficiently large $\eta$, there exist $\alpha$ such that

where the energy norm is defined as

And the boundedness is satisfied, there exist $\beta$ such that

the sum space $V(h)$ is defined as follow, $V(h) = \mathcal{{R}}(U_h)+H^{2}(\Omega)\cap H^{1}_{0}(\Omega)\subset H^{2}( \mathcal{T}_h)$.

This immediately gives the existence and the uniqueness of the weak form (2.3). The error estimate for the Poisson problem is given by the following theorem.

Theorem 3.4. [12,Theorem 1] Let $u\in H^{m+1}(\Omega)$ and $u_h$ be the solution of the Poisson problem and (2.3), respectively. Then there exists a positive constant $C$, such that

If Assumption 1 holds, then we can simplify(3.10) into

The above estimates are for general meshes. In cases of uniform meshes, we can obtain a better convergence rate at the barycenters or the face centers of the mesh.

3.2.   Superconvergence property

In one dimensional space, the meshes are equally distributed. Then refine the mesh by split one grid into two identical grid, which guarantees each grid points in the coarse mesh will also be the grid points in the fine mesh, see Figure 3.1.

A convergence rate $2(m+1)$ can be observed at some special points. In one dimensional case, the grid points and the barycenters of the grid have the superconvergence property. We define two norms for the estimates, which are related to the discrete point values. There are $N$ grid points on the uniform partition $\mathcal{T}_h$, which are the barycenters of elements, denoted by $x_{1},....,x_{N}$. We define the norm $\|\,{\cdot}\,\|_{{h}}$ that depends on the mesh as follows,

Obviously, $\|\,{\cdot}\,\|_{{h}}$ satisfies the triangle equality and the absolutely homogeneous. The positive definite property is also satisfied in the reconstruction space $V_h$ due to the barycenters are chosen as the sampling nodes.

Since there are two values on the grid points and the discontinuity, we investigate the average $\left\{\,{\cdot}\,\right\}$ at the grid points. Assume there are $M$ grid points on the uniform partition $\mathcal{T}_h$, denoted by $x_{1},....,x_{M}$, the mesh depended quantity $|\cdot|_{h}$ is defined as follows,

For one dimensional case, the following superconvergence result is observed numerically.

Proposition 1. If the partition $\mathcal{T}_h$ are uniform, $u$ and $u_h$ are the solutions in Theorem 3.4, and additionally $u\in H^{2m+1}(\Omega)$, then the superconvergence property with order $2m+1$ at the grid points and the barycenters are achieved, i.e.

We consider two types of uniform partitions in two dimensional spaces: uniform trianglar mesh and uniform square mesh. These meshes are demonstrated in Figure 2.1 and Figure 3.2, respectively. The special points which give the superconvergence property are the barycenters of each element. However, we also numerically observe the superconvergence phenomenon at the center of the element faces in uniform square mesh, where the convergence rate is $m+2$.

Proposition 2. If the partition $\mathcal{T}_h$ are uniformly distributed in two dimensional space, $u$ and $u_h$ are the solution in Theorem 3.4, and additionally $u\in H^{m+2}(\Omega)$, such that the superconvergence property with $m+2$ order at the barycenters of each elements are satisfied, i.e.

Meanwhile, if the partition $\mathcal{T}_h$ is uniform square mesh, the superconvergence also can be observed at center of element faces.

The superconvergence results in three dimensional space do not coincide with that in two dimensional spaces. We cannot observe the superconvergence property at the barycenters of elements. The domain $\Omega$ is partitioned into uniform tetrahedron or hexahedron meshes which are shown in Figure3.3. The superconvergence property with rate $m+2$ can be observed numerically at the centers of each face of the elements.

Proposition 3. If the partition $\mathcal{T}_h$ are uniformly distributed in three dimensional space, $u$ and $u_h$ are the solution in Theorem 3.4, and additionally $u\in H^{m+2}(\Omega)$, then the superconvergence property with rate $m+2$ at the center of each face of the elements are satisfied, i.e.

4.   Numerical experiments

In this section, we present numerical results for the superconvergence property of the PRDG method. All the meshes are uniformly distributed as was presented. We first demonstrate the relation between the sparsity pattern of the resulting linear system and the size of the element patch. We employ the uniform tetrahedron mesh to partition the cubic domain with $3072$ elements. The linear reconstruction is conducted with the patch size $7$ and $16$, and the quadratic reconstruction is conducted with the patch size $16$. Figure 4.1 shows the three sparsity patterns with different reconstructions and patch sizes.

Firstly, we observe that the size of all the linear system is $3072\times 3072$, regardless of the reconstruction. Then we compare the left two subfigures, which are the linear reconstruction with different patch sizes. The sparsity pattern changes when the element patch and the number of nonzero elements increases. However, the right two subfigures demonstrate the linear reconstruction and the quadratic reconstruction with the same patch size. The sparsity patterns and the number of nonzero terms coincide with each other. Generally, the sparsity of the resulting linear system is mainly determined by the size of the element patch. When the reconstruction order increases, larger element patch is requires, which will lead to slightly worse sparsity. In the numerical experiments, we usually take the element patch sightly larger than the DOFs of high order polynomial reconstruction requirement.

4.1.   One dimensional results

In one dimensional occasion, the exact solution is taken as $u(x) = \sin(2\pi x)$, the corresponding right hand side is $f = 4 \pi^2\sin(2\pi x)$ and the Dirichlet boundary condition $u(0) = u(1) = 0$ is weakly satisfied.

In Figure 4.2 and Table 4.1, we present the convergence rate for this Poisson equation which is solved by the PRDG method with different norms. At the mesh points and the barycenters of the grids, we obtain convergence rate $2m+1^{\text{th}}$ while the order of $L^2$ norm error is $m+1$. The results agree with the Proposition 3.1.

4.2.   Two dimensional results

We consider the Poisson equation with the Dirichlet boundary in the square domain. The exact solution is $u(x) = \sin(2\pi x)\sin(2\pi y)$, and the superconvergence property is slightly different from that in one dimensional case. Figure 4.3 and Table 4.2 shows the numerical results in uniform triangle meshes. The superconvergence appears at the barycenters of each element with rate $m+2$. However, the convergence rate at the centers of each element faces is consistent with the $L^2$ norm error.

Figure 4.4 and Table 4.3 present the numerical results for square meshes. The superconvergence property can be achieved at the grid points and the element face centers at the same time. The convergence rate is $m+2$ which is also coincides with Proposition 3.2.

4.3.   Three dimensional results

Finally, we present a three dimensional example. The Poisson equation with exact solution $u = \sin(2\pi x)\sin(2\pi y)\sin(2\pi z)$ is considered. Figure 4.5 and Table 4.4 shows the linear reconstruction numerical results. The superconvergence can only be achieved at the centers of element faces. The convergence rate is $m+2$ which agrees with Proposition 3.3 well.

5.   Conclusions

The superconvergence property of the PRDG method for elliptic problems is investigated numerically. Details of numerical implementations are presented and the sparsity patterns of resulting linear systems are demonstrated. The superconvergence results are achieved at the face centers or the barycenters of the element in function values with the uniform partitions.

Acknowledgments

The authors would like to thank the anonymous referees. They have very constructively helped to improve the original version of this paper.

The research is supported by the National Natural Science Foundation of China (Grant No. 11671312, 91630313), the Natural Science Foundation of Hubei Province (Grant No. 2019CFA007), and China Postdoctoral Science Foundation (Grant No. 2019M660558). The numerical calculations have been done on the supercomputing system in the Supercomputing Center of Wuhan University.