FINITE ELEMENT MULTISCALE METHODS FOR POISSON’S EQUATION WITH RAPIDLY VARYING HETEROGENEOUS COEFFICIENTS

An abstract framework for constructing finite element multiscale methods is presented. Using this framework we propose and compare two different multiscale methods, one based on the continuous Galerkin finite element method and one on the discontinuous Galerkin finite element method. In these multiscale methods the solution is split into coarse and fine scale contributions. The fine scale contribution is obtained by solving localized constituent problems on patches and is used to obtain a modified coarse scale equation. The localized constituent problems are completely parallelizable i.e, no communication between the different problems are needed. The modified coarse scale equation has considerably less degrees of freedom than the original problem. Numerical experiments are presented where the effect of the patch size of the local constituent problems as well as the convergence of the multiscale methods are investigated and compared for the proposed multiscale methods. We conclude that for a given accuracy and a fixed number of patches, smaller patches can be used for the discontinuous Galerkin multiscale method compared to the continuous Galerkin multiscale method.


INTRODUCTION
There are numerous applications which involves solutions that varies over several different scales, for example flow in porous media such as oil reservoir simulations and CO 2 storage.These, so called multiscale problems, are often impossible to solve with standard single mesh methods since the finest scale needs to be resolved to get a reliable result, see e.g.[5].
To resolve this problem several multiscale methods have been developed during the last two decades e.g., the Multiscale Finite Element Method (MsFEM) by Hou and Wu [9] and the Variational Multiscale Method (VMS) by Hughes [10].See also [8,7,12] and references therein for recent development and exposition.Using the framework of the Variational Multiscale Methods Larson and Målqvist introduced the Adaptive Variational Multiscale method [11].This method has further been developed in [12], where the framework for constructing multiscale methods used in this paper is presented and further discussed.
Lately, there have been a lot of interest in discontinuous Galerkin multiscale methods.Discontinuous Galerkin (DG) methods appeared in the 1970s; see [3,6] for some early work for elliptic problem and [4,14,15] for a literature review.A desired property with DG methods is that they admits good conservation properties of the state variable and are ideally suited for application to complex and irregular meshes.Conservation is a crucial property for multiscale problems.Recently proposed multiscale discontinuous Galerkin methods include e.g., [1] based on the MsFEM, and [2] based on the Heterogeneous Multiscale Method.
In this paper a continuous Galerkin multiscale method and a discontinuous Galerkin multiscale method for solving Poisson's equation with rapidly variable heterogeneous coefficients are studied.The continuous Galerkin version was first presented in [11], while the DG version is new.In the proposed multiscale method the solution is split into coarse and fine scale contributions.The fine scale contribution is obtained by solving localized constituent problems on patches and is used to obtain a modified coarse scale equation.Both a symmetric and a non-symmetric version of the modified coarse scale equation are presented.Numerical experiments are presented, where the size needed for the constituent problems to get a sufficient approximation as well as the convergence of the different multiscale methods, are investigated.We conclude that for a given accuracy and a fixed number of patches, smaller patches can be used for the discontinuous Galerkin multiscale compared to the continuous Galerkin multiscale method.On the coarse scale the discontinuous Galerkin multiscale method is approximating the L 2 -projection, rather than the nodal values, which is the case for continuous Galerkin multiscale method.The property of approximating the L 2 -projection is preferable in a multiscale setting.Also, DG has better conservation properties than CG.
The precise setting of the paper is the following.We consider the following model problem: where Ω ⊂ R d for d = 1, 2, 3, is a polygonal domain and α ∈ L ∞ (Ω), such that α > β > 0, β ∈ R has multiscale structure.Equation (1) has a unique solution u ∈ H 1 (Ω) up to a constant for each f ∈ L 2 (Ω) provided that Ω f dx = 0 is satisfied.Defining the L 2 -scalar product as (•, •) L 2 (ω) on a domain ω ⊆ Ω, the weak formulation of (1) reads: The rest of the paper is organized as follows.In Section 2, we present the different finite element methods needed to construct the multiscale methods.In Section 3 an abstract framework for constructing multiscale methods as well as the specific multiscale methods used in the numerical examples are proposed.Section 4, is devoted to some implementation details.Finally, in Section 5 numerical experiment are presented.

FINITE ELEMENT METHODS
Let K = {K} be a shape-regular mesh and let Γ denote the set of all edges (or faces in 3D) of the mesh K.The set Γ is the union of two disjoint subsets Γ = Γ I ∩ Γ B , where Γ I is the union of the interior edges and Γ B the union of the boundary edges.Given an interior edge e = ∂K + ∩ ∂K − ⊂ Γ I for K + , K − ∈ K, denote K + the element with the higher index and n as the outward unit normal of K + on e. Defining v + := v| ∂K + and v − := v| ∂K − , we set the average and jump operator as, for e ∈ Γ I and for e ∈ Γ B .Also, for a non negative integer p, we denote by P p (K), the set of all polynomials on K of total degree at most p.

Continuous Galerkin method
In the continuous Galerkin (CG) finite element discretization we are using a conforming approximation of the test space i.e.,

Discontinuous Galerkin method
In the discontinuous Galerkin method discretization we use a non-conforming approximation i.e., S h = {v ∈ L 2 (Ω) : v| K ∈ P r (K), K ∈ K} ⊂ V.The discontinuous Galerkin method reads: find u h ∈ S h such that where the bilinear form B dg : S h × S h → R and the linear functional F dg : S h → R are given by respectively; here h e := diam(e), and σ e ∈ R is a positive constant, depending on the variable α, large enough to make the bilinear form (7) coercive with respect to the natural energy norm.We refer, e.g., to [14,4] and references therein for details on the analysis of DG methods for elliptic problems.

ABSTRACT MULTISCALE METHOD
In the VMS framework, the fine scale finite element space, W h , is decoupled into coarse and fine scale contributions W h = W c ⊕ W f , where W c is associated with a coarse mesh K c .The split between the coarse and the fine scales is determined by an inclusion operator I c : W h → W c .The coarse and fine scale contributions are defined as, W c := I c W h and W f := (I − I c )W = {v ∈ W : I c v = 0}.There are several different chooses of I c e.g. the L 2 -projection onto W c or the nodal interpolant onto the coarse mesh.Let B : W ×W → R be a bilinear form, we can then define a multiscale map T : W c → W f from the coarse to the fine scale as The reference solution and the test function can be decomposed into a coarse and fine-scale contribution; The multiscale problem reads: find The fine scale component u f can be computed by letting v c = 0 in (10) and using the multiscale map (9).We arrive to the problem: find The coarse scale solution is obtained by letting v f = 0 in (10): find u c ∈ W c such that In ( 12), T u c and u f are unknown and obtained by solving ( 9) and (11).Note that Then a symmetric formulation of the coarse scale problem is obtained by considering The linear systems ( 12) and ( 13) has dim(W c ) unknowns, but (9) and (11) are equally hard to solve as the original problem and need to be approximated.

Localization of the multiscale method
Let N be the index set of all nodes, {x i }, in the mesh K c .Further, given that the coarse space is spanned by basis functions W c = span{φ j }, let M i be the index set of all φ j such that φ j (x i ) = 1, in the continuous setting M i = {i} and in the discontinuous case M i have several entries.For each basis function φ j we solve: find T φ j ∈ W f such that where φ j + T φ j can be viewed as a modified basis function.Because the fast decay of φ j + T φ j away from supp(φ j ), see [13] for the conforming case, we can solve (9) on small overlapping patches ω i ⊂ Ω for each basis function φ j where j ∈ M i .Defining W f (ω i ) to be W f restricted to the patch ω i , ( 9) is transformed to: for each i ∈ N and j ∈ M i find The term (11) can be handled in i similar fashion by splitting the right hand into local contributions using a partition of unity.The size of the patches is determined by adding a superscript L, ω L i , as in Definition 1.
Definition 1 Let {φ j : j = 1, . . ., dim(W c )} be the Lagrange basis (continuous or discontinuous) of W c .The sum Φ i := j∈M i φ j constructs a standard continuous Lagrangian basis function.We say that ω 1 i is an 1-layer patch, if ω 1 i = supp(Φ i ).Further, we say that ω L i is an L-layer patch if Finally, the set ω L i \ω L−1 i will be referred to as an L-ring.This is illustrated in Figure 2.
Example of a 1 layer patch ω 1 i and 2 layer patch ω 2 i around node i.

Continuous Galerkin multiscale method
The split between the coarse and fine scale spaces, V h = V c ⊕ V f , is realized by choosing the inclusion operator to be the nodal interpolant; I c = Π c .To keep the conformity of the method the fine scale problem is solved on patches using Dirichlet boundary condition.The multiscale problem reads: for all i ∈ N find The modified coarse scale equations is then formulated as: find U c ∈ V c such that for the non-symmetric formulation and as for the symmetric formulation.The solution to the multiscale problem is

Discontinuous Galerkin multiscale method
Exploiting the discontinuous nature of S h the split between the coarse and fine spaces, S h = S c ⊕ S f , is realized by choosing the inclusion operator to be the element wise L 2projection onto S c ; I c = P c .This is more natural in a multiscale setting since the coarse scale solution approximate the average on each coarse element rather than the nodal values.The discontinuous nature of S h also allows for using Neumann boundary conditions on the fine scale problems.With V c = span{φ j }, we need to solve the fine scale problem: for all i ∈ N and j ∈ M i where The modified coarse scale equations are formulated as: find U c ∈ S c such that for the non-symmetric formulation or for the symmetric formulation.The solution to the multiscale problem is

IMPLEMENTATION
In the proposed multiscale method, the fine scale problem is perfectly parallelizable i.e., no communication between different fine scale problems are required.Algorithm 1 shows how the multiscale methods can be implemented.Note that the outer for-loop is perfectly parallel.An schematic overview is given in Figure 2 where the lines between the boxes represent communication.Also, note that the assembly of the coarse stiffness matrix and load vector is also done in parallel, in the fine scale problems.The extra constraints on the fine scale problems are realized using Lagrange multipliers Algorithm 1 Multiscale Method 1: Initialize the coarse mesh with mesh size H. 2: Let the fine mesh size be h = H/2 n and the size of the patches L. Determine the patch ω L i .

5:
for j ∈ M i do 6: Compute the fine scale contribution for the modified basis functions T φ j .

7:
end for 8: Compute the right hand side correction U f i .9: end for 10: Solve the modified coarse scale problem to obtain U c .

a,f,Ω
Split into local problems, transfer data

Decay of modified basis functions
Consider the domain ω L i , for L = 1, . . ., N .On ω N i for N = 8, let K c be a coarse mesh consisting of 16 × 16 elements and K f be a fine mesh consisting of 128 × 128 elements.Let T L φ j ∈ W f (ω L i ) be the solution of computed on ω L i and extended by 0 in Ω \ ω L i , where φ j ∈ M i , is a basis function on the coarse scale.Three types of permeabilities, called Ones, Period, and SPE, are used.For One, a = 1, for Period, α = 1 or α = 0.1 with a period of 1/64 in x-direction, and SPE, data is taken from the 31st layer permeability data in the tenth SPE comparative solution project1 and illustrated in Figure 3.The aspect ration is a max /a min = 5.9823 • 10 5 .The decay of the coarse modified basis function φ j + T L φ j is illustrated by computing T L φ j for L = 1, . . ., N − 1 using T N φ j as a reference solution.The space W f and the bilinear form B(•, •), are defined as V f and B cg (•, •) for the continuous Galerkin multiscale method, and as S f and B dg (•, •) for the discontinuous Galerkin multiscale method.Exponential decay, in the broken energy norm for L = 1, . . ., N when N = 4, is observed in Figure 4.The fast decay motivates us to solve the constituent problems on patches ω L i ⊂ Ω using a small number of L-rings.This, in turn, means less computational work and a smaller overlap between the localized problems.The DG method converges faster than CG to to the reference solution in the relative broken energy   23) for different permeability using continuous Galerkin (solid line) and discontinuous Galerkin (dashed line).norm (24).Hence, smaller patches are needed for solving the local problems using DG than CG to achieve the same accuracy.

Comparison of the continuous and discontinuous Galerkin multiscale methods
Consider the model problem (1) on the unit square Ω = (0, 1) × (0, 1).Let K be a reference mesh with M N × M N elements, and K c a coarse mesh of N × N elements i.e., each coarse elements is further subdivided into M ×M elements.In the numerical experiment N = 16 and M = 8.Let, f (x, y) = −1 for {0 < x, y < 1/128}, f (x, y) = 1 for {127/128 < x, y < 1}, and f = 0 otherwise, be the forcing function.The same permeabilities, Ones, Rand and SPE, as in Section 5.1 are used.In the numerical experiments all patches, ω L i , are of the same size, L, and for each iteration L is increased by one.The continuous Galerkin multiscale method and the discontinuous Galerkin multiscale method are compared, see Figure 5.We conclude: • To obtain a given accuracy, in the relative broken energy norm (24), the discontinuous Galerkin multiscale method requires approximately one layer less than the continuous Galerkin multiscale method.For a comparison of the degrees of freedom required for the fine scale problems, see Table 1.
• This is a bit unfair comparison since the reference solution is the DG respectively CG solution computed on the fine scale.DG has a more enriched test and trial space and may give a better approximation than CG because of the discontinuous permeability coefficients.
• On the coarse scale the discontinuous Galerkin multiscale method is approximating the L 2 -projection rather than the nodal values, which is the case for continuous Galerkin multiscale method.This is preferable in a multiscale setting.
• The DG method has better conservation properties which is an important property in many multiscale applications.

Figure 2 .
Figure 2. Implementation scheme of the discontinuous Galerkin multiscale method.

Figure 3 .
Figure 3. Permeability structure for SPE (c) in log scale.

Table 1 .
Degree of freedom for the fine scale problems d (16n) d