A posteriori error estimates for a multi-scale finite-element method

We are interested in the discretization of a diffusion problem with highly oscillating coefficient, by a multi-scale finite-element method (MsFEM). The objective of this method is to capture the multi-scale structure of the solution via local basis functions which contain the essential information on small scales. In this paper, we perform an a posteriori analysis of this discretization. The main result consists of building error indicators with respect to both small and large meshes used in this method. We present a numerical test in which the experiments are in good coherency with the results of analysis.


Introduction
The a posteriori error analysis is a powerful tool to improve the quality of approximated solutions of a model of partial differential equations. The estimates obtained in this context allow us to perform a self-adaptation of the mesh and to reach a desired accuracy, via a fixed tolerance. The norm of the error is bounded by estimators which depend only on the mesh size, the data, the approximated solution, and which are explicitly computable. The pioneering work of the a posteriori error estimates, in the context of finite-element method, was done by Rheinbolt in 1978 (see Babuška et andRheinboldt 1978a, b).
In this work, we are interested in the methods of multi-scale finite element (MsFEM) treated in Carballal Perdiz (2010), Efendiev and Hou (2009), Hou and Wu (1997), Lozinski et al. (2013), and Ahmed Blal (2014). Generally, this method is based on two main ingredients: building multi-scale local basis functions on a coarse mesh T H , and coupling them with a global variational formulation on a fine mesh T h , providing an accurate approximation of the solution. Multi-scale basis functions are designed to capture the characteristics of the multi-scale solution and contain information on small scales. They are built from those of standard finite element in the coarse mesh, such that they have the same support and satisfy the equation L φ = 0 on each element of T H , where L is the principal operator in the model. There are also variants of MsFEM: MsFVM ("multi-scale finite volume method", see Lee et al. 2002), MsMFEM ("multi-scale mixed finite-element method", see Chen and Hou 2003), MsFVEM ("multi-scale finite volume element method", see Hou 2009), and DG-MsFEM ("discontinuous Galerkin multi-scale finite-element method", see Efendiev and Hou 2009).
In this work, we develop a posteriori error estimates for finite-element multi-scale method for a diffusion problem with a mesh adaptation. For this, as in the construction of the method, we will develop estimates for the overall solution related to the coarse mesh and couple these estimates with those related to the fine mesh. In Aarnes and Efendiev (2006), we find adaptive techniques for a Finite Volume Multi-scale method based on indicators on some physical criteria. In Abdulle and Nonnenmacher (2011), the authors give residual type indicators for homogenization problems. In Henning et al. (2014), we find a residual error estimate for another type of MsFEM, that is a Petrov-Galerkin MsFEM with over-sampling. In this method, the finite-dimensional "coarse-scale" subspace and the continuous fine-scale space are defined through a projection operator from H 1 0 (Ω), and the authors use a reconstruction operator and a corrector operator to define the discrete problem. The a posteriori estimates are on the error between the exact solution and the reconstructed discrete one, and the indicators are expressed with these operators. There are three types of indicators. The first one is related to the coarse mesh and is similar to one of ours. The second is based on the jump of the fluxes through the inter-element of the fine mesh and the last one is retard the over-sampling.

Formulations
Let Ω ⊂ R 2 be a bounded polygonal domain with a Lipschitz-continuous boundary Γ = ∂Ω. We use the standard space L 2 (Ω) equipped with the usual norm · 0 , together with the Sobolev space H 1 (Ω) of functions in L 2 (Ω), such that their first derivatives (in distribution sense) belong to L 2 (Ω), equipped with the norm v 1 := v 2 0 + ∇v 0 1/2 . As usual, is equipped with the norm |v| 1 := ∇v 0 . We consider the diffusion problem associated with the operator Lv := −div (ν∇v) defined by: We assume that f ∈ L 2 (Ω), ν ∈ C 0,1 (Ω), the space of Lipschitz-continuous functions, and there exist two positive constants ν m and ν M satisfying: Let us note that ν ∈ L ∞ (Ω) is sufficient to have the existence and uniqueness of the solution in H 1 0 (Ω), but for the need of the a posteriori error analysis, we assume more regularity for ν.
We denote by |·| ν the energy norm defined by: If D is a subset of Ω, we use the notation: and also the following notations: We remark that: The existence and uniqueness of the solution of (1) are obtained using the Lax-Milgram theorem. In many applications, the coefficient ν may present a highly oscillatory character. To obtain the large-scale solutions accurately and efficiently without resolving the small-scale details, we will use a Multi-scale Finite-Element Method (MsFEM) (see Babuška and Osborn 1983;Carballal Perdiz 2010;Efendiev and Hou 2009;Hou and Wu 1997;Lozinski et al. 2013;Ahmed Blal 2014), an approach that captures the multi-scale structure of the solution. Let T H a regular mesh of Ω: where H = max K H K and H K is the diameter of the element K . We denote by N H and E H , respectively, the set of all internal nodes x i for i = 1, 2, . . . , N H , and the set of all edges E i for i = 1, 2, . . . , N E of T H , excluding the edges on the boundary ∂Ω. The first-order standard finite-element space based on T H and approximating V is denoted by P 1 (T H ) and defined by: where P 1 (K ) and P 1 (E)) are the spaces of polynomial functions with degree ≤ 1 in K and E, respectively. As usual, the basis functions {ψ i H } N i=1 associated with P 1 (T H ) satisfy {ψ i H }(x j ) = δ i j , and the support of {ψ i H }, denoted by ω i , is the union of the β i elements K i,d , such that x i is one of their vertices: The multi-scale basis functions are defined (see Efendiev and Hou 2009;Hou and Wu 1997), for i = 1, . . . , N H , as the function Φ i H ∈ H 1 (ω i ) satisfying the following problems, for d = 1, . . . , β i : with the notation L K := L| K . The Multi-scale Finite-Element (MsFEM) space adapted to the operator L is then defined by: The MsFEM approximation of the problem (1) is now given by: To determine the multi-scale basis functions, we have now to solve the problems (4). To do so, we consider a new finer mesh. For each K ∈ T H , we consider a regular conforming mesh T h (K ) all independently, and define T h as the overall fine mesh: Let us note that each T h (K ) is a conforming mesh, but the global mesh can be nonconforming one.
We denote by E h (K ) the set of all edges of T h (K ) and by E h the set of all edges of T h excluding the edges on the boundary ∂Ω. The edges from E h can be divided into two groups: those forming the boundary of coarse triangles and thus constituting the coarse mesh edges (denoted by E H ) and those internal for some triangles K ∈ T H forming the set E 0 h (K ). In the sequel, we need also the following notations. For any E ∈ E H , we denote by K 1 (E) and K 2 (E), the two coarse adjacent triangles, such that E = K 1 (E) ∩ K 2 (E). Let, moreover, h (E) be the set of the triangles of both fine meshes T h (K 1 ) and T h (K 2 ) that touch E. We define also E b h (E) as the set of internal edges of both fine meshes T h (K 1 ) and T h (K 2 ) that touch E and are not in E. For i = 1, . . . , N H , the sets ω i are also the supports of the multi-scale basis functions Φ i H and we put: with the notations a K (u, v) = K ν∇u · ∇v, and for K ∈ T H : We define the MsFEM discrete space by: and the discrete MsFEM approximation of the problem (1) is defined by: The multi-scale solution u h H of the global discrete problem (10) is therefore written as a linear combination of multi-scale basis functions: (2010) and Lozinski (2010), we can prove the following a priori error estimate.

As in Carballal Perdiz
Proposition 1 Let u be the solution of (1) and u h H be the solution of (10). We have the following error estimate: where V L := v ∈ H 1 0 (Ω) such that ∀E ∈ E H v | E is linear , and C 0 be a constant which depends only on the geometry of K .
In the other hand, since V H is also a subspace of H 1 0 (Ω), by Cea's Lemma, we have: , satisfies the problem: and by applying the Poincaré's inequality and (3): where C(K ) is the Poincaré's constant in K which can be written as C(K ) = C 0 H K , where C 0 depends only on the geometry of K and is independent of H K , its diameter.
where we use the shorthand notation: for x ≤ Cy with positive constant independent of x, y, and meshes. To derive the lower bound of the error, we need special functions introduced by Verfürth (1996) and called bubble functions. Using the barycentric coordinates, λ K i , for i = 1, 2, 3, associated with an element K ∈ T H , these bubble functions are defined by: where E = K 1 ∩ K 2 . We have the following lemma (see Verfürth 1996, lemma 3.3, page 66): v be a polynomial function in K , and σ be a polynomial function on E. We have:

Interpolation operators
We consider here operators of Clément type, related, respectively, to the spaces V H and where Φ i H and Φ i h,H , i = 1, . . . N H are the multi-scale basis functions defined, respectively, by (4) and (7). These operators satisfy the following properties.
where ω K (resp. ω E ) is the union of all the triangles of T H that share at least one node with K (resp. E).
The proof is given in the Annex. We shall also need a Clément interpolation operator Π h : H 1 (K ) → P 1 (T h (K )), defined relatively to the P 1 finite elements on the fine mesh T h (K ) for any coarse triangle K . We take a version of such an operator constructed by (see Scott and Zhang (1990) where T is any triangle from T h (K ), e any internal edge from E 0 h (K ),ẽ any boundary edge from the mesh T h (K ), ω T (resp. ω e ) is the union of all the triangles of T h (K ) that share at least one node with T (resp. e) and γẽ is the union of all the boundary edges that share at least one node withẽ. We define now another operator. For v ∈ H 1 (K ), let Π h v, such that: (26) satisfies the following estimates: where α T and α e are boundary switches, i.e., α T = 1 (resp. α e = 1) if T (resp. e) touches ∂ K ; otherwise, α T = α e = 0, and where ω T (resp. ω e ) is the union of all the triangles of T h (K ) that share at least one node with T (resp. e).
The proof is in the Annex.

A posteriori error estimates
In this section, we will give different a posteriori error estimates. The first ones are related to the approximation of the problem (1) in the space V H spanned by the multi-scale functions, and given by the problem (5). Obviously, the function u H will not be computed directly, but over the approximation of the multi-scale functions. Here and below, we use the notation (14) and the notation (2) for D equal K , E, T , e, etc. . .. For E ∈ E h , we denote by n E the unit, normal, outward-pointing vector field and by ν∇ n E u H E the jump of ν∇u H · n E across the element E. In the following, the notation x y means that x ≤ Cy with positive constant independent of x, y, ν, and meshes.

Estimates related to multi-scale basis functions
Theorem 1 Let f ∈ L 2 (Ω), u be the solution of (1), and u H be the solution of (5), and let η K , for K ∈ T H defined by: Then, we have: where f K is an approximation of f .
Proof 1. We begin by proving the upper bound of the error (30). By taking w = u − u H and for any w H ∈ V H , the Galerkin orthogonality and an integration by parts in each K ∈ T H give: By Cauchy-Schwarz inequality and since div(ν∇u H ) = 0 in K , we obtain: (21) and use the estimates of Proposition 2 and (3) to obtain: which gives (30) since w = u − u H . 2. To prove (31), the upper bound of the indicator, we use the technique of bubble functions which is now standard. Let w K = b K f K , where b K is the bubble function in K defined by (15). By (17), we have: Since L K u H = 0 in K , and K w K f dx − K ν∇u · ∇w K dx = 0, and then, using Green's formula, we have: The Cauchy-Schwarz and inverse inequalities, the definition of w K , and (3) give: , and then: Let now w = ν∇(u H − u). First of all, we remark that for all element K ∈ T H and all E ∈ ∂ K , we have divw = f in K and ν∇ n E u E = 0. In the other hand, since L K u H = 0 in K , u H | ∂ K is regular and K is convex, we deduce by regularity theorem (Cf. Grisvard 1985) that u H | K ∈ H 2 (K ), and so, w · n is in L 2 (E). By Lemma 1, we therefore have: and therefore:

Estimates for the fully approximated problem
In this section, we give a posteriori error estimates for the discrete MsFEM approximation problem given by (10). The a posteriori indicators are defined, for all K ∈ T H , all T ∈ T h (K ), and all E ∈ E H , by: We will prove the upper and lower bounds of the error by these indicators.

An upper bound for the error
Theorem 2 Let f ∈ L 2 (Ω) and u h H be the solution of (10). We have:

Proof of Theorem 2
Let u h H be the solution of (10) and w = u − u h H . Then, by the Galerkin orthogonality, we have for any w h H ∈ V h H : Using (1) and an integration by parts over all the triangles T ∈ T h (K ), for all K ∈ T H , yields: Since the edges from E h is divided into two groups: E H , the set of the boundary of coarse triangles, and E 0 h (K ), the set of internal edges from E h (K ) for some triangle K ∈ T H , the last sum can be rewritten as: (20), and use Proposition 2 and the inequality (3) to estimate the terms in T 1 : To estimate the terms in T 2 , we shall use the definition of the space V h H . Since u h H is an element of V h H , it satisfies on any K ∈ T H [see (7), (8), and (9)]: An integration by parts over all the triangles Thus, one can write for any (26). With the help of Proposition 3, one gets: We make the estimates taking into account that v h | ∂ K = 0 and Π h (w) = v h + w h H , and using Proposition 3, we obtain: which give: and since v h | ∂ K = 0, we obtain: and α e are defined in Proposition 3. We now gather the terms in the sums above and use Cauchy-Schwarz, and then, (24) to get: Using again the inequality (3), we obtain: Combining together the bounds for T 1 and T 2 and recalling that w = u − u h H , this leads to (36)

A lower bound for the error
To have an optimal a posteriori error estimate of the error, we have to prove that the indicators η K ,h , δ K , and ξ E defined in (33), (34), and (34), respectively, are locally lower bounds of the error.
Theorem 3 Let f ∈ L 2 (Ω), u h H the solution of (10) and η K ,h , δ K , and ξ E defined by (33), (34), and (34), respectively. The following estimates hold for all K ∈ T H and all E ∈ E H :

Remark 2
In Theorem 3, the error indicators are increased by a quantity related to the H /h ratio; this means that we do not need to refine the fine mesh, so that this ratio is not big enough.

Implementation
In this section, we will show a numerical test which confirm the reliability and effectiveness of the a posteriori error estimates for the MsFEM method developed in the previous section. All the calculations presented are done with the software FreeFem++ (Cf. Hecht and Pironneau 2021).
To establish our MsFEM method, we need to build two meshes of Ω, T H the coarse one and T h the fine one. We first build the conforming coarse mesh T H = K with H = max K H K and H K is the diameter of K . The fine mesh is then constructed in the following manner. Each element K of T H will be meshed with a step equal to H K /n K , where n K is an integer fixed in advance, by eventually taking into account the a priori information on the problem. We thus obtain a conforming mesh T h (K ) = T ⊂K T . We denote the diameter of the element T by h T and h = max T h T . The global fine mesh of the domain is obtained by catering these last meshes and will not be necessary conforming: In the first subsection, we present the resolution algorithm for the MsFEM scheme developed in this work and then the mesh adaptation algorithm. In the second subsection, we present the numerical calculations for the problem (45).

Algorithms
The resolution algorithm has the following structure:

Algorithm 1 (T H , n, f ) (Resolution algorithm)
-a coarse mesh T H = K , representing as well as possible the geometry of the problem; -{S i , i = 1, . . . , N } the set of vertices of T H ; -n a non-negative integer; -the data f ; 1. for each K ∈ T H , build the meshes T h (K ) with h T ≈ H K /n; 2. for each i, 1 ≤ i ≤ N , calculate the basis function Φ i in each K (a part of the support of this function) as the solution of (7); 3. solve the discrete problem (10) (which is a N × N system), to calculate the solution u h H .
The indicators developed in Theorem 1 and Theorem 2 are robust tools to find the elements K and T , where the error between the exact solution and the approximated solution is too large, and then must be refined to improve the quality of the numerical solution. Let η K , δ K ,h , and ξ E,h defined respectively by (29), (34), and (34). First, we refine T H , using only the indicators η K , and ξ E,h . Indeed, δ K ,h is much smaller than the others and can be used to adapt the meshes related to the basis functions (fine mesh).
The mesh adaptation process has the following general structure.
H , the error between the reference solution and the solution obtained by the standard Finite-Element Method with the coarse mesh. We will consider different norms and put:

Numerical test
We consider the initial coarse mesh T H given in Fig. 1(c) and the corresponding fine mesh given in Fig. 1(d), where we have chosen n = 4. By Algorithm 1(T H , n, f ), we obtain the solution represented by its isovalues in Fig. 1(b). We note that with this initial meshes, the MsFEM solution is less regular than the reference one. The first adaptation of the coarse mesh T H is made by cutting by 4 the elements K to obtain a finer conforming mesh given in Fig. 3a and the corresponding fine mesh given in Fig.  3b. The solution related to the new meshes is given in Fig. 3f. We note that the modifications are made near the two refined elements. Figure 3 (c), (d), and (e) shows that the isovalues of indicators are smaller near the refined elements.
This procedure of Algorithm 2 is continued, until the estimators reach T ol H . This is done at the 6 th iteration for the MsFEM solution. In Fig.4, we find all the information about the solution MsFEM at this iteration. The solution FEM presented in this figure is obtained by starting with the initial coarse mesh and using adaptivity via the estimators I K , until these estimators reach T ol H . This is done at the 8 th iteration and the results are also presented in Fig. 4.
We notice that MsFEM captures the oscillations much better than the FEM, with a lower cost, since the size of the MsFEM system is smaller than the FEM system (see Table 1). We see in Fig. 4 that the error given by FEM is large in the oscillation zone, while MsFEM gives a better approximation in this zone.   Table 1, we give a summary of the numerical results for this test. The notation I T E means "iteration". The suffix "max" or "min" means the maximum value or the minimum value, respectively.
At first sight, we note that the error decreases with the iterations of the mesh adaptation as well as the indicators and their gaps between the maximum and the minimum values. These indicators detect the location where the MsFEM solution is not accurate enough. One of the gaps can sometimes increase, but globally, it decreases. For the same tolerance of the error, MsFEM needs less iterations than FEM, and the solution is better in the different norms, as shown in Table 1. Furthermore, this strategy leads to an equidistributivity of the error, since the gap between the indicators η K ,h and ξ E,h decreases from the first iteration to the last one as well as the others gaps (see Table 1).
In Fig. 5a, we plotted the log of the maximum of the indicator η K as a function of the number of degree of freedom, and in Fig. 5b, we find a comparison between the error MsFEM and the error FEM as functions of the number of iterations.
We tested the robustness of our MsFEM adaptive algorithm when the parameter ε goes to 0. In Fig. 5c, we see that the error decreases in a similar way for the three values ε = 10 −2 , ε = 10 −3 or ε = 10 −4 . Furthermore, we do not need the have the fine mesh very fine. Indeed, in Fig. 5c, we can see that before adaptation procedure, the error is the same for the three fine meshes related to h = H /4, h = H /8, and h = H /16, and it is not recommended to take the mesh T h very fine when we adapt T H .

Conclusion
In this work, we have presented an a posteriori analysis of a multi-scale finite-element method (MsFEM), for a diffusion problem with highly oscillating coefficients. We derived upper and lower bounds for the approximation error and presented a numerical test confirming their performance in regard to their efficiency and reliability. The indicators obtained are of residual type. Those related to the fine mesh are, in some how, standard and represent the residual equation and the jump of the normal derivative of the solution through the interfaces of the fine mesh. The others are also of residual type and take into account the linearity constraint on the interfaces of the coarse mesh.
It was noticed in the work Lozinski et al. (2013) that the error decreases with the step H of the coarse mesh and a good precision is reached for a relatively small number of basis functions, whereas the difference between the errors associated with the different steps h of the fine mesh is minimal. As a consequence, in the numerical test, we opted to refine only the