Computational aspects of the multiscale discontinuous Galerkin method for convection-diffusion-reaction problems

1.   Introduction

The finite element method (FEM) is a numerical technique for finding approximate solutions to boundary value problems for partial differential equations. It uses the subdivision of the whole domain into simpler parts, called finite elements, and variational methods from the calculus of variations to solve the problem by minimizing an associated error function. The discontinuous Galerkin finite element method is a variant of the classical finite element method. Its main difference with classical finite element methods is the continuity of the solution across element interfaces. The discontinuous Galerkin (DG) method does not require the continuity of the solution along edges. Since this leads to ambiguities at element interfaces, the technique from finite volume methodology (FVM), namely the choice of numerical fluxes, has been introduced. From this point of view, DG methods combine features of the finite element methods and finite volume methods. Thus DG methods have several advantages. For instance, DG methods are highly parallelizable and very well suited for handling adaptive strategies. However, DG methods have more degrees of freedom than classical finite element methods.

On the other hand, the solution of the convection dominated problem has typically singularity and to resolve it one requires very fine meshes in the domain or very high order polynomials in the approximate spaces, which produces very large degrees of freedom especially in the DG method.

Over the decades, several variants of DG method such as hybridizable DG (HDG), DG with Lagrange multiplier (DGLM), multiscale DG (MDG), have been developed to reduce the degrees of freedom of DG. Concerning MDG, it was first introduced by Hughes et. al. and investigated for advection-diffusion equation [3,6]. They have introduced extra streamline diffusion term in the setting of the MDG to treat the advection term. In [2], MDG was also introduced for elliptic problem without theoretical analysis. In [10], one of the authors has presented the MDG for convection-diffusion-reaction equations without artificial viscosity.

Relating other multiscale methods for convection diffusion equations, we refer the readers to, for example, [4] and the references cited therein.

MDG has the computational structure of the continuous Galerkin (CG) method based on the variational multiscale idea (see [2]), which is indeed a DG method. Storage and computational efforts are reduced significantly in high order approximation.

In this paper, we study computational aspects of the MDG [10]. Especially, we investigate the matrix structure of the MDG. The MDG solution is obtained by composition of the DG and the inter-scale operator. We show that the composition results to the matrix product of the DG matrix and the inter-scale matrix corresponding to the local problem on the element. We apply ILU preconditioned GMRES to effectively solve the resulting global system.

This paper is organized as follows. In Section 2, we introduce the model problem called the convection-diffusion-reaction equation. In Section 3, we introduce finite element spaces for the MDG method. In Section 4, we define the parameters and describe the MDG method for the model problem. In Section 5, we form the computational structure of the MDG. Finally, in Section 6, we show the numerical results with convection dominated problems.

2.   Convection-diffusion-reaction problems

Let $\Omega$ be a bounded open domain in $\mathbb{R}^m$, $m = 2,3$ with a polyhedral boundary $\partial \Omega$. We consider the convection-diffusion-reaction equation

where $f \in L^{2}(\Omega)$ and $c \in L^{\infty}(\Omega)$ are real valued, ${\bf b}$ is the solenoidal vector field defined on $\overline{\Omega}$, and ${\bf A}$ is a symmetric matrix whose entries are bounded, piecewise continuous real valued functions defined on $\overline{\Omega}$, with

By $\vec{n}(x)$ we denote the unit outward normal vector to $\partial \Omega$ at $x \in \partial \Omega$. We divide the boundary as follows: Let

The sets $\partial \Omega_{-}$ and $\partial \Omega_{+}$ will be referred to as the inflow and outflow boundaries, respectively. Clearly we see that $\partial \Omega = \partial \Omega_{o} \cup \partial \Omega_{-} \cup \partial \Omega_{+}$. If $\partial \Omega_{0}$ is nonempty, we shall further divide it into disjoint subsets $\partial \Omega_{D}$ and $\partial \Omega_{N}$ whose union is $\partial \Omega_{o}$ with $\partial \Omega_{D}$ is nonempty and relatively open in $\partial \Omega$. We supplement (2.1) with the boundary conditions:

and adopt the (physically reasonable) hypothesis that ${\bf b} \cdot \vec{n} \ge 0 \mbox { on } \partial \Omega_{N}$ whenever $\partial \Omega_{N}$ is nonempty. We assume that the values of the diffusion coefficient ${\bf A}$ and the velocity field ${\bf b}$ ensure well-posedness of (2.1) and (2.4). Additional assumptions on these coefficients will be set later on (see [10]).

3.   Finite element spaces

In this section, we introduce the finite element spaces for the MDG method. We recall the space $L^{2}(\Omega)$ of all square integrable functions and the space $H^{1}(\Omega)$ of all functions in $L^{2}(\Omega)$ that have square integrable derivatives. To define approximation spaces, we consider a uniformly regular partition $\mathcal{T}_h$ of $\Omega$ into finite elements $E$.

Let $\mathcal{T}_h$ be a regular family of triangulations of $\Omega$ in the sense that there exists a $\kappa>0$ such that $\frac{h}{h_{\min}}<\kappa$, where $h_{\min} = \min(h_{E})|_{E\in\mathcal{T}_h}$, $h_{E}$ is the diameter of $E \in \mathcal{T}_h$, and $h = \max(h_{E})|_{E\in\mathcal{T}_h}$. We assume that $\mathcal T_{h}$ contains only regular nodes, that is, each element vertex is also a vertex to all adjacent elements and there are no hanging nodes. The elements $E \in \mathcal{T}_h$ are either triangles and/or quadrilaterals in two dimensional space or tetrahedra and/or hexahedra in three dimensional space. The vector $\vec{n}_E(x)$ is also used as the outward unit normal to $\partial E$ at $x \in \partial E$, which coincides with the one on $\partial \Omega$. We denote by $\mathcal {E}_h$ the set of all edges of $\mathcal{T}_h$, by $\mathcal {E}_h^i$ the set of all interior edges, and by $\mathcal{E}_h^o = \mathcal {E}_h \backslash \mathcal{E}_h^i$ the set of all outer edges.

In consistent DG methods, the solution values are coupled by generalized flux functions across the edges and they appear by the jumps and averages. Let $e$ be an interior edge shared by two elements $E_{1}$ and $E_{2}$ and let $\vec{n}_1$ and $\vec{n}_2$ be the unit normal vectors on $e$ pointing exterior to $E_{1}$ and $E_{2}$, respectively. With $\phi_{j} : = \phi_{\mid_{E_{j}}}$, following [1,2,10], we define

For a vector valued function $\vec{\tau}$ which is piecewise smooth on $\mathcal {T}$, with analogous meaning for $\vec{\tau}_{1}$ and $\vec{\tau}_{2}$, we define

Notice that the jump $[[ \phi]]$ of the scalar function $\phi$ across $e \in \mathcal {E}_h^i$ is a vector parallel to the normal to $e$, and the jump $[[\vec{\tau}]]$ of the vector function $\vec{\tau}$ is a scalar quantity. For $e \in \mathcal {E}_h^o$, the set of boundary edges, we let

We do not require either of the quantities $\{ \phi \}$ or $[[\vec{\tau}]]$ on boundary edges, and leave them undefined there.

We now assign to $\mathcal {T}_h$ the broken Sobolev space of order $p$,

equipped with the broken Sobolev norm $\left ( \sum_{E \in \mathcal {T}_h} \| w \|^{2}_{H^{p}(E)} \right )^{1/2}$. When $p = 0$, we denote by $\| \cdot \|_{E}$ the norm in $L^{2}(E)$ and similarly by $\| \cdot \|_{\partial E}$ the norm in $L^{2} (\partial E)$. For the construction of the MDG, we first consider the continuous finite element space

and then associate with it the discontinuous approximation spaces:

$W_{h}(\Omega)$ is then viewed as being obtained from $\overline{W}_{h}(\Omega)$ by releasing interelement continuity constraints. Also $W_{h}(\Omega)$ is a formal union of the local spaces $W_{h}(E)$.

4.   Multiscale discontinuous Galerkin method

In this section, we introduce the MDG method (see [6,10]). To do that, we start with a DG method for (2.1) and (2.4) given as follows: Find $u_{h} \in W_{h}(\Omega)$ such that

where $B(\cdot,\cdot)$ is a bilinear form defined on $H^{2}(\mathcal{T}_h) \times H^{2}(\mathcal{T}_h)$ and $L(g_{D},g_{N},f;\cdot)$ is a linear functional defined on $H^{2}(\mathcal{T}_h)$. We assume that the values of $B(\cdot,\cdot)$ on the interior edges are determined by averages and jumps. That is,

where $B_{E}(\cdot,\cdot)$ is a bilinear element form defined on $E \in \mathcal{T}_h$, $B_{e}^{o} (\cdot,\cdot)$ is a boundary edge bilinear form defined on the boundary edges, $B_{e}^{i}(\cdot,\cdot)$ is an interior edge bilinear form that depends only on the values from the two elements that share interior edge $e$ in $\mathcal{E}_{h}^{i}$, and $(\widetilde{u}_{h},\widetilde{v}) = (\langle [[u_{h}]],[[{\bf A}\nabla u_{h}]], {\bf b}\cdot [[u_{h}]],\{ u_{h} \},\{ {\bf A} \nabla u_{h} \} \rangle, \langle [[v]], [[a\nabla v]], {\bf b} \cdot [[v]], [[v]], \{ v \}, \{ {\bf A}\nabla v \}\rangle)$. We here note that a large class of consistent DGs satisfies (4.2) (see [1,6]).

Now, to define the local problem of the MDG, we decompose $u_h$ of $\overline{u}_h$ and $u'_h$. That is, let

Then, (4.1) takes the below forms:

Coarse scale equation:

Fine scale equation:

By treating the function $\overline{u}_{h}$ as data, we rewrite (4.4) as follows: Find $u'_{h} \in W_{h}(\Omega)$ such that

Take $v' \in W_{h}(E)$ and, assume, without loss of generality, that $v'$ is extended to be zero outside of $E$. Using (4.2), we then see that (4.5) is equivalent to the following: For each $E \in \mathcal{T}_h$, find $u'_{h} \in W_{h}(E)$ such that, $\forall v' \in W_{h}(E)$,

(4.6) relates fine scales to the coarse scales, but remains coupled to the continuous elements through the numerical flux terms in (4.6). MDG defines the local problem in a way that the fine scales are expressed in terms of the coarse scales that will uncouple (4.6) on inner boundaries.

We now note that, for $\tau \in W_h(E)$,

where $v^i$ and $v^o$ are the values of $v$ from the inside and the outside of $E$, respectively. We then see that, $u'_h \in W_h(E)$,

Similarly, since

we have that

Substituting (4.7)-(4.8) into (4.6), local problem of the MDG is ended up as follows: $\forall \overline{u}_{h} \in \overline{W}_{h}$, find $u'_{h} \in \cup_{E\in \mathcal{T}_h} W_{h}(E)$ such that

where $b_{E}(\cdot, \cdot)$ and {$l_{E}(\cdot, \cdot\,;\cdot)$} are linear with respect to its first argument and affine with respected to its second and third arguments (see [7]):

We note that (4.9) is a DG formulation for the local problem $-\nabla \cdot({\bf A}(x)\nabla u_h)$ $+{\bf b}(x)\cdot \nabla u_h+c(x)u_h = f$ on $E$ having weakly imposed boundary condition on $\partial E$, given that $g_{D} = \overline{u}_h$ on $\partial \Omega_{D}$ and $g_{N} = ({\bf A}\nabla \overline{u}_h)\cdot n$ on $\partial \Omega_{N}$.

We now denote by $\mathfrak{I}:\overline{W}_{h}(\Omega) \times L^{2}(\Omega) \rightarrow W_{h}(\Omega)$ the operator which associates to each $(\overline{u}_h,f)\in \overline{W}_{h}(\Omega) \times L^{2}(\Omega)$ the solution $u_h'$ of the local problem on each $E\in\mathcal{T}_h$. Let

and

Global MDG method is then given to find $u^{MDG}_h \in \mathfrak{I}(\overline{W}_h,f)$ such that

Remark. Concerning the analysis of the stability and the error estimates of the MDG method, we refer the reader to [10].

5.   Computational aspects of MDG

In this section we form the matrix equation of the MDG method. To make the presentation simple, we consider the case of standard NIPG (Nonsymmetric Interior Penalty Galerkin) DG method [1,10].

Let $\mathcal{T}_h = \{E_k\}_{k = 1}^L$. Let $\{ \phi_i^{E_k} \}_{i = 1} ^{N_k}$ be the nodal basis of $W_h(E_k)$ and let $\{ \overline{\phi}_i \}_{i = 1} ^{N}$ be the nodal basis functions of $\overline W_h(\Omega)$. We then see that

where $\overline{N}_k$ is the number of degrees of freedom of the continuous finite element on $E_k$.

By following the process of the previous section, one can obtain the local problem (4.9), which is now given as follows (see [7,10]): Find $u_h^E\in W_h(E)$, such that, for $v_h^\prime\in W_h(E)$, $\forall E\in{\mathcal{T}}_h,$

where

Now, let $u_h^{MDG}\in\mathfrak{I}(\overline{u}_h;f)$ be the solution to the MDG method (4.12) for (2.1) and (2.3) and let

where $\overline U = (\overline U_1,\overline U_2,\,\cdots,\,\overline U_N)^T,$ $U^k = (U^k_1, U^k_2,\,\cdots,\,U^k_{N_k})^T$, $F^k = (F^k_1, F^k_2,\,\cdots,\, F^k_{N_k})^T$, for $k = 1,\cdots, L$, and $\chi_E$ is the characteristic function having value one on $\overline E$ and zero otherwise. Let $\overline U^k = (\overline U_1^k,\overline U_2^k,\,\cdots,\,\overline U_{\overline {N}_k}^k)^T$ be the vector of the continuous nodal values on ${E_k}$.

Now, let $\widetilde S^k$ be the $N_k\times\overline N_k$ inter-scale matrix of the local problem (5.2), that is,

where $\widetilde S^k = ({\mathcal{L}}^{k,m})^{-1}\overline{\mathcal{R}}^m,$ $({\mathcal{L}}^{k,m})_{ji} = b_E(\phi^{E_k}_i,\phi^{E_m}_j)$, and $(\overline{\mathcal{R}}^m)_{jr} = b_{\partial E}(\overline \phi_r, \phi^{E_m}_j)$. We denote by $(\Lambda)_i$ and $(\Lambda)_{ij}$ the $i$-th row vector and the $ij$ component of matrix $\Lambda$, respectively.

Now, let $S^k$ be the $N_k\times{N}$ zero extension of $\widetilde S^k$ and let ${\breve{S}}^k$ be the $N\times{N}$ zero extension of $S^k$. We denote by ${\breve{\phi}}^{E_k}_i$ the zero extension of $\overline \phi^{E_k}_i$ into $\Omega$. We then see that

Noting

and from (5.4), we express $u_h^{MDG}$ in terms of $\overline U$ and we have that

Similarly, we consider $v_h\in\mathfrak{I}(\overline{u}_h;0)$ such that, for $1\le l\le L$,

where $\overline I^l$ is the $\overline N_l\times 1$ vector such that $\overline I^l_p = \delta^l_p$ and $\overline I$ is the $N\times 1$ zero extension of $\overline I^l$. Substituting (5.8)-(5.9) into (4.12) and noting (5.1) and (5.7), we find that

Let $A^{l,k}$ be the $N_l\times N_k$ matrix such that

and $\widetilde A^{l,k}$ be the $N\times N$ zero extension of $A^{l,k}$ such that

We then by noting (5.7), that (5.10) is rewritten as follows:

Here we have used the fact that ${\breve{S}}^l$ and $\widetilde A^{l,k}$ are the $N\times N$ matrices and $\overline I$ and $\overline U$ are the $N\times 1$ vectors. We now note, by the definition of ${\breve{S}}^l, S^l$, and $A$, that

Considering $l = 1,\cdots, L$ and also from (5.12)-(5.13), we see that the MDG matrix of (4.12), $A^{MDG}$, is given by

Remark. We see, by (5.13), that

where

Here $S$ is the $M\times N$ inter-scale matrix corresponding to the inter-scale operator $\mathfrak{I}(\overline{u}_h;f)$, where $M = \sum_{k = 1}^LN_k$, and $A$ is the $M\times M$ DG matrix. We note that the MDG matrix is a product of the interscale matrices of local problems and the DG matrix.

Remark. We further see, from (5.14)-(5.15), that

In the next section, we compute the MDG solution in the following way:

$\mathbf{1.}\,\,\,$ Do for all $E\in{\mathcal T}_h$:

● construct $A^{k,l}$: compute the local mass and stiffness matrices of usual DG.

● construct $\widetilde S^k$: solve the local problem on $E^k$.

● calculate $(\widetilde S^l)^TA^{k,l}{\widetilde S^k}$.

$\mathbf{2.}\,\,\,$ Sum up the computations following the global numbering.

$\mathbf{3.}\,\,\,$ Solve the global system by ILU preconditioned GMRES.

6.   Numerical results

In this section, we test the MDG with convection dominated problems. When the convection strongly dominates the diffusion ($k = 10^{-6}$), a fast linear solver is needed to solve the large system. We apply ILU preconditioned GMRES method to solve the matrix equation effectively. We compute $L^2$-error and convergence order defined by

where $u(\cdot,\cdot)$ denotes the exact solution and $u^h$ is an approximate solution.

6.1.   Diffusion coefficient $k = 10^{-4}$

We solve the problem (2.1) and (2.4) in the case of the diffusion coefficient $k = 10^{-4}$ with the following data: ${\bf A} = k{\bf I}$ (here ${\bf I}$ is the 2 by 2 identity matrix), ${\bf b} = (2,3)$, $c = 1$, $\Omega = (0,1) \times (0,1)$ and $g = 0$ (see [9]).

We take $f$ such that the exact solution is given by

where

The graph of the exact solution and the approximate solution of DG with uniform mesh $h = 1/32$ are given in Figure 2 and Figure 3, respectively. The DG solution shows spurious oscillation on the inner layer and in the direction of convection. To reduce these errors, we take smaller meshes or use high order degree of the approximate polynomials. We have tested the MDG and DG to compare the results. The results are given in Table 1. When $P_1$ elements are used, the DG solution requires the degrees-of-freedom six times more than the one of MDG solution. The graphs of the approximate solutions with $h = 1/128$ are given in Figure 4. Similarly, when $P_4$ elements are used, the DG solution requires the degrees-of-freedom three times more than the one of MDG solution. The graphs of the approximate solutions with $h = 1/32$ are given in Figure 5.

6.2.   Diffusion coefficient $k = 10^{-6}$

In this subsection, we consider the case of $k = 10^{-6}$ in diffusion coefficient.

We take $h = 1/512$ and $P_1$ element. The DG solution shows spurious oscillation as seen in Figure 6. It is numerically shown that $P_1$ approximation requires very small $h$ to resolve the spurious oscillation and the computation is not efficient. Figure 7 shows the approximate solutions without oscillation when $P_4$ and $h = 1/256$. DG solution requires degrees-of-freedom three times more than the one of MDG to produce similar results.

In this test, we have also applied the ILU preconditioner to the matrix equation. In Tables 2-3, we compare the results with and without ILU preconditioning for DG and MDG, respectively. We see that ILU preconditioned GMRES effectively solve the large system. As seen in Tables 2-3, DG requires degrees-of-freedom five times more than the ones of MDG to obtain the qualitatively similar solution. CPU time of the MDG solution was also significantly reduced.

6.3.   Diffusion coefficient $k = 10^{-6}$ with polynomials of mixed degrees

If the convection strongly dominates, then the linear approximation is not efficient to resolve the spurious oscillations. Figure 9 shows the graphs of the solutions of MDG with $h = 1/32$ and $1/256$.

In this subsection, we apply high order approximation only in the region where detailed information is needed. By observing the oscillations occur along the convected direction, we simply apply high order polynomial in the convected direction. We show that MDG is very flexible in increasing the polynomial degree $p$ where it is needed.

As seen in the previous tests, oscillation starts to develop at the inner layer and propagates to the gray area in Figure 8. We apply different orders of polynomials in different regions (see [5,8,10,13]). We apply $P_1$ in the grid part and $P_4$ in the shadow part, respectively, in Figure 10.

Matrix structure and $P_1$ and $P_4$ elements are given in Figure 11. $3\times 3$, $15\times 15$ block matrices corresponding to $P_1$ and $P_4$ elements, respectively, are placed on the diagonal.

We compare the result with the one of DG with $P_4$ elements. MDG solution exhibits similar behavior of DG with much less degrees of freedom. The figures are given in Figure 12. We compute the $L^2$-error for the MDG solution with polynomials of mixed degrees. The results are given in Table 4.

As a conclusion, the MDG method numerically shows an advantage in adaptive $p$ refinements. ILU GMRES is applied to the resulting system. Combination of $P_1$ element in the smooth region and higher order element in region where more detailed information are needed, can be a reasonable candidate in refinement. Further research on the adaptive refinement will be done in the future.