A discontinuous Galerkin multiscale method for convection-diffusion problems

# A discontinuous Galerkin multiscale method for convection-diffusion problems

Daniel Elfverson
###### Abstract

We propose an discontinuous Galerkin local orthogonal decomposition multiscale method for convection-diffusion problems with rough, heterogeneous, and highly varying coefficients. The properties of the multiscale method and the discontinuous Galerkin method allows us to better cope with multiscale features as well as interior/boundary layers in the solution. In the proposed method the trail and test spaces are spanned by a corrected basis computed on localized patches of size , where is the mesh size. We prove convergence rates independent of the variation in the coefficients and present numerical experiments which verify the analytical findings.

## 1 Introduction

In this paper we consider numerical approximation of convection-diffusion problems with possible strong convection and with rough, heterogeneous, and highly varying coefficients, without assumption on scale separation or periodicity. This class of problems, normally refereed to as multiscale problem, are know to be very computational demanding and arise in many different areas of the engineering sciences, e.g., in porous media flow and composite materials. More precisely, we consider the following convection-diffusion equation: given any we seek such that

 −∇⋅A∇u+b⋅∇u =fin Ω, (1)

is fulfilled in a weak sense, where is the computational domain with boundary . The multiscale coefficients will be specified later. There are two key issues which make classical conforming finite element methods perform badly for these kind of problems,

• the multiscale features of the coefficient need to be resolved by the finite element mesh and

• strong convection leads to boundary and interior layers in the solution which need to be resolved.

To overcome the lack of performance using classical finite element methods in the case of multiscale features in the coefficient many different so called multiscale methods have been proposed, see [25, 26, 23, 7, 13, 12, 10, 11] among others, which perform localized fine scale computations to construct a different basis or a modified coarse scale operator. Common to the aforementioned approaches is that the performance of the method rely strongly on scale separation or periodicity of the diffusion coefficients. There is also approaches which perform well without scale separation or periodicity in the diffusion coefficient but to high computational cost by either having to solve eigenvalue problems [2] or where the support of the localized patches is large [37, 4]. See also [38].

In the variational multiscale method (VMS) framework [25, 26] the solution space is split into coarse and fine scale contribution. This idea was employed for multiscale problems in a adaptive setting for classical finite element in [31, 34, 32] and to the discontinuous Galerkin (DG) method in [14]. A further development is the local orthogonal decomposition (LOD) method, see [36, 20, 19] for classical finite element and [15] for DG methods. The LOD operates in linear complexity without any assumptions on scale separation or periodicity and the trail and test spaces are spanned by a corrected basis function computed on patches of size . The LOD has e.g. been applied to eigenvalue problems [35], non-linear elliptic problems [21], non-linear Schrödinger equation [17], and in Petrov-Galerkin formulation [16].

There is a vast literature on numerical methods for convection dominated problems, we reefer to [28, 24, 27], among others. There has also been a lot of work on DG methods, we refer to [39, 33, 3, 29] for some early work and to [8, 22, 40, 9] and references therein for recent development and a literature review. DG methods exhibit attractive properties for convection dominated problems, e.g., they have enhanced stability properties, good conservation property of the state variable, and the use of complex and/or irregular meshes are admissible. For multiscale methods for convection-diffusion problems, see e.g. [1, 41, 18].

In this paper we extended the analysis of the discontinuous Galerkin local orthogonal decomposition (DG-LOD) [15] to convection-diffusion problems. For problems with strong convection using the standard LOD won’t suffice, since convergence can no longer be guarantied. Instead we propose to include the convective term in the computations of the corrected basis functions. We prove convergence results under some assumptions of the magnitude of the convection and present a series of numerical experiment to verify the analytic findings. For problems with weak convection it is not necessary to include the convective part [21].

The outline of this paper is as follows. In section 2 the discrete setting and underlying DG method is presented. In section 3 the multiscale decomposition, the DG-LOD, and the corresponding convergence result are stated. In Section 4 numerical experiments are presented. Finally, the proofs for some of the theoretical results are given in Section 5.

## 2 Preliminaries

In this section we present some notations and properties frequently used in the paper.

### 2.1 Setting

Let for be a polygonal domain with Lipschitz boundary . We assume that: the diffusion coefficients, , has uniform spectral bounds , defined by

 (2)

and the convective coefficient, , is divergence free

 ∇⋅b(x)=0 a.e. x∈Ω. (3)

We denote .

We will consider a coarse and a fine mesh, with mesh function and respectively. Furthermore, we assume that the fine mesh resolve and that the coarse mesh do not resolve the fine scale features in the coefficients. Let , for , denote a shape-regular subdivision of into (closed) regular simplexes or into quadrilaterals/hexahedras (), given a mesh function defined as for all . Also, let denote the -broken gradient defined as for all . For simplicity we will also assume that is conforming in the sense that no hanging nodes are allowed, but the analysis can easily be extend to non-conforming meshes with a finite number of hanging nodes on each edge. Let be the reference simplex or (hyper)cube. We define to be the space of polynomials of degree less than or equal to if is a simplex, or the space of polynomials of degree less than or equal to , in each variable, if is a (hyper)cube. The space of discontinuous piecewise polynomial function is defined by

 Pp(Tk):={v:Ω→R∣∀T∈Tk,v|T∘FT∈Pp(^T)}, (4)

where , is a family of element maps. We will work with the spaces . Let denote the -projection onto . Also, let denote the set of all edges in where and denote the set of interior and boundary edges, respectively. Given that and are two adjacent elements in sharing an edge , let be the the unit normal vector pointing from to , and for let be outward unit normal of . For any we denote the value on edge as . The jump and average of is defined as, and respectively for , and for . For a real number we define its negative part as .

Let denote any generic constant that neither depends on the mesh size or the variables and ; then abbreviates the inequality .

### 2.2 Discontinuous Galerkin discretization

For simplicity let the bilinear form , given any mesh function , be split into two parts

where represents the diffusion part and represents the convection part. The diffusion part is approximated using a symmetric interior penalty method

where is a constant, depending on the diffusion, large enough to make coercive. The convective part is approximated by

 ach(u,v):=(b⋅∇hu,v)L2(Ω)+∑e∈Eh(Ω)(be[u],[v])L2(e) (7) −∑e∈Eh(Ω)(νe⋅b[u],{v})L2(e)+∑e∈Eh(Γ)((νe⋅b)⊖u,v)L2(e),

where upwind is imposed choosing the stabilization term as [5].

The following definitions and results are needed both on the fine and coarse scale, for this sake let . The energy norm on is given by

 |||v|||2k,d =∥A1/2∇kv∥2L2(Ω)+∑e∈Ekσek∥[v]∥2L2(e), (8) |||v|||2k,c =∑e∈Ek∥b1/2e[v]∥2L2(e), |||v|||2k =|||v|||2k,d+|||v|||2k,c.

From Theorem 2.2 in [30] we have that for each , there exist an averaging operator with the following property

 ∥∇k(v−Ickv)∥L2(Ω)+∥k−1(v−Ickv)∥L2(Ω)≲∑e∈Ek1k∥[v]∥2L2(e). (9)

In the error analysis we will also need a localized energy norm, defined in a domain (aligned with the mesh ) as

 |||v|||2k,d,ω =∥A1/2∇kv∥2L2(ω)+∑e∈Eke∩¯ω≠0σek∥[v]∥2L2(e), (10) |||v|||2k,c,ω =∑e∈Eke∩¯ω≠0∥b1/2e[v]∥2L2(e), |||v|||2k,ω =|||v|||2k,d,ω+|||v|||2k,c,ω.

## 3 Multiscale method

In this section we preset the multiscale decomposition and extend the results in [15] to convection-diffusion problems. For the constants in the convergence results to be stable we assume the following relation of the convective term

 O(∥Hb∥L∞(Ω)α)≤1 (11)

How the magnitude of (11) affects the convergence of the method is investigated in the numerical experiments.

### 3.1 Multiscale decomposition

In order to do the multiscale decomposition the problem is divided into a coarse and a fine scale. To this end let and , with the respective mesh function and , denote the two different subdivisions, where is constructed by some (possible adaptive) refinements of .

The aim of this section is to construct a coarse finite element space based on , which takes the fine scale behavior of the data into account. We assume that the mesh resolves the variation in the data, i.e., the solution to: find such that

 ah(uh,v)=F(v)for all v∈Vh, (12)

gives a sufficiently good approximation of the weak solution to (1). Note however that never have to be computed in practice, it only acts as a reference solution. We introduce a coarse projection operator and let the fine scale reminder space be defined by the kernel of , i.e.,

 Vf:={v∈Vh∣ΠHv=0}⊂Vh. (13)

The coarse projection operator has the following approximation and stability properties.

###### Lemma 1.

For any and , the approximation property

 H|−1T∥v−ΠHv∥L2(T)≲α−1/2|||v|||h,T, (14)

and stability estimate

 |||ΠHv|||H≲Cs|||v|||h, (15)

is satisfied, with

 Cs=(C2A+∥Hb∥L∞(Ω)α)1/2. (16)
###### Proof.

The approximation property follows directly from [15, Lemma 5]. Let be a Clément type interpolation operator proposed in [6, Section 6] which satisfy

 ∥∇CHu∥L2(T)+∥H−1(u−CHu)∥L2(T)≲∥∇u∥L2(ω1T), (17)

where are the union of all elements that share a edge with . We define the conforming function using averaging operator in (9). We obtain

 |||ΠHv|||2H=∑T∈TH∥A1/2∇(ΠHv−Π0v)∥2L2(T) (18) +∑e∈Eh(Ω∪ΓD)(σH∥[vc−ΠHv]∥2L2(e)+∥b1/2e[vc−ΠHv]∥2) ≲∑T∈THβ(1H2∥v−Π0v∥2L2(T)+(1H2+∥b∥L∞(T)H)∥vc−v∥2L2(T))

using that is the -projection onto constants, a trace inequality, and stability of . Next, using that

 ∥CHIchv−v∥L2(Ω) ≤∥CHIchv−Ichv∥L2(Ω)+∥Ichv−v∥L2(Ω) (19) ≲H∥∇Ichv∥L2(Ω)+∥Ichv−v∥L2(Ω) ≲α−1/2H|||v|||h

in (18) concludes the proof. ∎

The following lemma shows that for every there exist a (non-unique) in the preimage of which is conforming.

###### Lemma 2.

For each , there exist a such that , , and .

###### Proof.

Follows directly from [15, Lemma 6], since . ∎

The next step is to split any into some coarse part based on , such that the fine scale reminder in the space is sufficiently small. A naive way to do this splitting is to use a -orthogonal split. An alternative definition of the coarse space is . A set of basis functions that span is the element-wise Lagrange basis functions where for simplexes or for quadrilaterals/hexahedra. The space is known to give poor approximation properties if does not resolve the variable coefficients in (1). We will use another choice, see [36, 15], based on , to construct a space of corrected basis functions. To this end, we define a fine scale projection operator by

 ah(Fv,w)=ah(v,w)for all w∈Vf, (20)

and let the corrected coarse space be defined as

 VmsH:=(1−F)VH. (21)

The corrected space are spanned by corrected basis functions which can be computed as: for all find such that

 ah(ϕT,j,v)=ah(λT,j,v)for all v∈Vf. (22)

Note that, . From (21) we have that any can be decomposed into a coarse and a fine scale contribution, .

###### Lemma 3 (Stability of the corrected basis function).

For all , the following estimate

 |||ϕT,h−λT,j|||h≲Cϕβ1/2∥H−1λT,j∥L2(Ω) (23)

holds, where .

###### Proof.

Let , where , , from Lemma 2. We have

Using that the diffusion part in (24) of the bilinear form is continuous in with the constant , Lemma 2, and a inverse inequality, we get

For the convection part in (24), we have

 (b⋅∇h(ϕT,h−λT,j),bT,j)L2(Ω) (26) ≲∥Hb⋅∇h(ϕT,h−λT,j)∥L2(Ω)∥H−1bT,j∥L2(Ω) ≲∥Hb∥L∞(Ω)∥∇h(ϕT,h−λT,j)∥L2(Ω)∥H−1λT,j∥L2(Ω),

and obtain

 |||ϕT,h−λT,j|||h≤Cϕβ1/2∥H−1λT,j∥L2(Ω). (27)

with . ∎

### 3.2 Ideal discontinuous Galerkin multiscale method

An ideal multiscale method seeks such that

 ah(umsH,v)=F(v)for all v∈VmsH. (28)

Note that, to construct in the space a variational problem has to be solved on the whole domain for each basis function, which is not feasible for real computations. The following theorem shows the convergence of the ideal (non-realistic) multiscale method.

###### Theorem 4.

Let be the solution to (12), and be the solution to (28), then

 |||uh−umsH|||≲C1α−1/2||H(f−ΠHf)||L2(Ω) (29)

holds, with

See Section 5. ∎

### 3.3 Discontinuous Galerkin multiscale method

The fast decay of the corrected basis functions (Lemma 6), motivates us to solve the corrector functions on localized patches. This introduces a localization error, but choosing the patch size as (Theorem 7) the localization error has the same convergence rate as the ideal multiscale method in Theorem 4. The corrector functions are solved on element patches, defined as follows.

###### Definition 5.

For all , let be a patch centered around element with size , defined as

 ω0T :=int(T), (30) ωLT :=int(∪{T′∈TH∣T∩¯ωL−1T≠0}),L=1,2,….

See Figure 1 for an illustration.

The localized corrector functions are calculated as follows: for all find such that

 ah(ϕLT,j,v)=ah(λT,j,v),for all v∈Vf(ωLT). (31)

The decay of the corrected basis function is given in the following lemma.

###### Lemma 6.

For all , where is the solution to (22) and is the solution to (36), the following estimate

 ∣∣∣∣∣∣ϕT,j−ϕLT,j∣∣∣∣∣∣h≲C2γL|||λT,j−ϕLT,j|||h (32)

holds, where is the size of the patch, , , and , where is a generic constants neither depending on the mesh size, the size of the patches, or the problem data.

###### Proof.

See Section 5. ∎

The space of localized corrected basis function is defined by . The DG multiscale method reads: find such that

 ah(ums,LH,v)=F(v)for all v∈Vms,LH. (33)

An error bound for the DG multiscale method using a localized corrected basis is given in Theorem 7. Note that it is only the first term in Theorem 7 that depends on the regularity of .

###### Theorem 7.

Let and be the solutions to (12) and (33), respectively. Then

 |||u−ums,LH|||h≤ |||u−uh|||h+Ccα−1/2∥H(f−ΠHf)∥L2(Ω) (34) +C5∥H−1∥L∞(Ω)Ld/2γL∥f∥L2(Ω)

holds, where is the size of the patches, is a constant defined in Theorem 4, and , where is defined in Lemma 13, and and are defined in Lemma 6.

See Section 5. ∎

###### Remark 8.

Theorem 7 is simplified to,

 |||u−ums,LH|||h≤|||u−uh|||h+C1∥H∥L∞(Ω). (35)

given that the patch size is chosen as with an appropriate and . In the numerical experiments we choose .

###### Remark 9.

If the convective term is small it is not necessary to include it in the computation of the correctors [21]. Instead the following correctors can be used: for all find such that

This gives the right convergence results if

 O(∥b∥L∞(Ω)α)=1 (37)

compared to (11) if the convective term is included.

## 4 Numerical experiment

We consider the domain and the forcing function . The localization parameter which determine the size of the patches is chosen as , i.e., the size of the patches are . Consider a coarse quadrilateral mesh, , of size , . The corrector functions are solved on sub-grids of the quadrilateral mesh, , where . We consider three different permeabilities: , which is piecewise constant with respect to a Cartesian grid of width in y-direction taking the values or , and which is piecewise constant with respect to a Cartesian grid of width both in the x- and y-directions, bounded below by and has a maximum ratio . The permeability is taken from the layer in the SPE 10 benchmark problem, see http:www.spe.org/web/csp/. The diffusion coefficients and are illustrated in Figure 2.

For the convection term we consider: , for different values of .

To investigate how the error in relative energy-norm, , depends on the magnitude of the convection we consider: and with . Figure 3 shows the convergence in energy-norm as a function of the coarse mesh size for the different values of .

Also, to see the effect of heterogeneous diffusion of the error in the relative energy-norm, , we consider: Figure 4 which shows the error in relative energy-norm using and and Figure 5 which shows the error in relative energy-norm using and .

We obtain convergence of the DG multiscale method to a reference solution in the relative energy-norm, , independent of the variation in the coefficients or regularity of the underlying solution.

## 5 Proofs from Section 3

In this section we state the proofs of the main results which was postponed from in section 3. To this end we start by proving some technical lemmas in Section 5.1 which we use to prove the main results in Section 5.2.

### 5.1 Technical lemmas

In the proofs of the main results, Theorem 4, Lemma 6, and Theorem 7, we will need some definitions and technical lemmas stated below.

Continuity of the DG bilinear form for convective problems can be proven on a orthogonal subset of . The space is an orthogonal subset of but on a coarse scale.

###### Lemma 10 (Continuity in (Vh×Vf) and (Vf×Vh)).

For all, or in , it holds

 a(v,w)≲Cc|||v|||h|||w|||h (38)

where

 Cc=CA+∥Hb∥L∞(Ω)α−1. (39)
###### Proof.

Since is continuous in with the constant , continuity in follows from . For the convective part , we have

 ac(v,w) =∑T∈Th(b⋅∇v,w)L2(T)+∑e∈Ek(be[v],[w])L2(e) (40) −∑e∈Ek(Ω)(νe⋅b[v],{w})L2(e)+∑e∈Ek(Γ)((νe⋅b)⊖v,w)L2(e) ≲∑T∈Th(∥b∥L∞(T)∥∇v∥L2(T)∥w−ΠHw∥L2(T)) +∑e∈Ek(∥b∥L∞(e)h−1/2∥[v]∥L2(e)∥w∥L2(S+∪S−)).

where and . Using a discrete Cauchy-Schwartz inequality and summing over the coarse elements, we get

 ac(v,w) ≲α−1/2∥Hb∥L∞(Ω)|||v|||h∥H−1(w−ΠHw)∥L2(Ω), (41) ≲∥Hb∥L∞(Ω)α−1|||v|||h|||w|||h,

which concludes the proof for . The proof of is obtained by first integrating by parts. ∎

The following cut-off function will be frequently used in the proof of the main results.

###### Definition 11.

The function , for , is a cut off function fulfilling the following condition

 ζd,DT|ωdT =1, (42) ζd,DT|Ω∖ωDT =0, ∥[ζd,DT]∥L∞(Eh(T)) ≲∥h∥L∞(T)(D−d)H|T,

and , for all .

For the cut off function has the following stability property.

###### Lemma 12.

For any and from Definition 11, the estimate,

 |||ζd,DTv|||h≲Cζ|||v|||h,ωDT, (43)

holds, where .

###### Proof.

For the diffusion part we use the following result from [15],

 |||(1−ζd,DT)v|||h