On heterogeneous coupling of multiscale methods for problems with and without scale separation

In this paper, we discuss partial differential equations with multiple scales for which scale resolution is needed in some subregions, while a separation of scale and numerical homogenization is possible in the remaining part of the computational domain. Departing from the classical coupling approach that often relies on artificial boundary conditions computed from some coarse grain simulation, we propose a coupling procedure in which virtual boundary conditions are obtained from the minimization of a coarse grain and a fine-scale model in overlapping regions where both models are valid. We discuss this method with a focus on interface control and a numerical strategy based on non-matching meshes in the overlap. A fully discrete a priori error analysis of the heterogeneous coupled multiscale method is derived, and numerical experiments that illustrate the efficiency and flexibility of the proposed strategy are presented.


Introduction
The past few years have witnessed a growing number of new numerical schemes for multiscale problems. Broadly speaking, the numerical challenge for the approximation of such problems is to avoid scale resolution, i.e., the use of a fine mesh that resolves the smallest scale in the problem. Indeed, such direct approaches are often computationally too expensive for practical applications. In this paper, we consider in a polygonal domain Ω ⊂ R d , d = 1, 2, 3, with boundary Γ = Γ D ∪ Γ N , the model problem −div (a ε ∇u ε ) = f in Ω, where f ∈ L 2 (Ω), g D ∈ H 1/2 (Γ D ), and g N ∈ L 2 (Γ N ), and where a ε ∈ (L ∞ (Ω)) d×d is a highly oscillatory tensor that satisfies, for 0 < α < β, α|ξ | 2 ≤ a ε (x)ξ · ξ , |a ε (x)ξ | ≤ β|ξ | ∀ξ ∈ R d , for a.e. x ∈ R.
For such a problem, one can identify two broad classes of multiscale methods, namely numerical methods that seek to approximate an effective solution of the original problems in which the small scales have been averaged out. The existence of such effective solutions usually relies on homogenization theory [9,21]. The attractivity of such methods is the possibility to obtain numerical approximations that correctly describe the macroscopic behavior of the multiscale problem at a cost that, however, is independent of the smallest scale. Another class of multiscale methods aims at building coarse basis functions that encode the multiscale oscillations in the problem. This class of methods usually comes with a cost that is no longer independent of the small scale, but the construction of the basis functions can be localized and, once constructed, this basis can often be reused in a multi-query context. We refer to [5] for a review and references of the first class of methods and to [15,20,23] for review and references of the second class of methods. The framework that makes the first class of methods efficient is that of a simultaneous coupling of a macro-and a micro-method. In such approach, a separation of the scales is often required. The second class of methods solves the fine-scale problems on overlapping patches, and it has recently been shown that convergences can also be obtained without assuming scale separation [20,23].
In this contribution, we address an intermediate situation between separated and nonseparated scales in the following sense: we assume that in a subset ω 2 of the computational domain Ω the macro-/micro-upscaling strategy can be applied but that in an other part ω of the domain one needs full resolution of the scales. Here we assume that this second domain is sufficiently small so that standard resolved finite element method (FEM) can be used. In the region ω 2 , we chose to use the finite element heterogeneous multiscale method (FE-HMM) [2]. While our method easily generalizes to multiple regions with and without scale separation, we assume here for simplicity that Ω = ω 2 ∪ ω. The main issue for such a coupling strategy is to set adequate boundary conditions at the interfaces of both computational domains. We note that such problems have numerous applications in the sciences; we mention, for example, heterogeneous structures with defects [10,16] or steady flow problems with singularities [17]. Coupling strategies between fine-scale and upscaled models have already been studied in the literature, for example, in [27] where a precomputed global homogenized solution is used to provide the boundary conditions in the fine-scale subregions. More recently, a coupling strategy based on an L 2 projection of the homogenized solution onto harmonic fine-scale functions has been discussed [8].
The aim of this paper is to pursue the study of a new approach that we have proposed in [6,7] relying on an optimization-based coupling strategy. By introducing small overlapping region ω 0 between ω 2 and ω, where both fine-scale and homogenized model are valid, we consider the unknown boundary conditions for both models as (virtual) control and minimize the discrepancy of the solutions from the two models in the overlap. Two possible scenarios are illustrated in Fig. 1. Such ideas have appeared earlier in the literature for coupling different type of partial differential equations [18] or for atomistic-to-continuum methods [28]. We also note, related to our method, the recent work on the coupling of local and nonlocal diffusion models [13].
We briefly describe the main contribution of this paper. First, in [7], the theory and the numerics have been developed for the cost function · L 2 (ω 0 ) , called distributed observation in the classical terminology of optimal control. Here we consider the cost function · L 2 (Γ 1 ∪Γ 2 ) , called boundary or interface control. Such controls can reduce the cost of Two scenarios for the domain decomposition of Ω the iterative method to solve the optimality system compared to the cost of solving the optimization problem with distributed observation [14]. Second, in [7] we used the same mesh in the overlap ω 0 with the consequence of having to use a mesh size that scales with the fine mesh of ω. Here we discuss the use of independent meshes in the overlap through appropriate interpolation techniques. As a result, we have again a significant reduction in the computational cost of the coupling as the macroscopic numerical method in ω 2 = Ω \ ω does not need an increasing number of micro-solvers as the mesh in ω 1 is refined. Finally, numerical examples were carried out in [7] only for the situation where ω Ω (Fig. 1, left); here we discuss also the scenario for which ∂ω∩∂Ω = ∅ (Fig. 1, right).
The paper is organized as follows. In Sect. 2, we describe the model problem, introduce the two minimization costs functions considered in this paper, and give an a priori error analysis between the coupled and the fine-scale solutions. In Sect. 3, we define the multiscale numerical discretization of the optimization problem and perform a fully discrete a priori error analysis. Finally, Sect. 4 contains several numerical experiments that illustrate the theoretical results and the performance of the new coupling strategy.

Problem formulation
Let ω ⊂ Ω be the region without scale separation and ω 0 be the overlap region. Assume that Γ 1 = ∂ω 1 \Γ and Γ 2 = ∂ω 2 \Γ are Lipschitz continuous boundaries. We decompose the tensor a ε of problem (1) into a ε = a ω + a ε 2 , where a ε 2 = a ε 1 ω 2 and a ω = a ε 1 ω are tensors with and without scale separation, respectively. The tensor a ε 2 H-converges toward an homogenized tensor a 0 2 [25]. Further, we set a 1 = a ε 1 ω 1 , u 1 = u ε 1 , and u 2 = u 0 2 . The heterogeneous control restricted to Dirichlet boundary controls is given by the following problem: find u ε 1 ∈ H 1 (ω 1 ) and u 0 2 ∈ H 1 (ω 2 ), such that 1 2 u ε 1 − u 0 2 2 H is minimized under the following constraints, for i = 1, 2, where the boundary conditions θ i , called the virtual controls, are to be determined; (H, · H ) is a Hilbert space specified below. Further, we define the space of admissible Dirichlet controls The strategy is to solve a minimization problem in the space of admissible controls, where we minimize the cost In this paper, two Hilbert spaces (H, · H ) are considered.
The solutions are split into where (v ε 1 , v 0 2 ) are called the state variables and satisfy, for i = 1, 2, where The solutions u ε 1,0 and u 0 2,0 exist and are unique, thank to the Lax-Milgram lemma, and the solutions v ε 1 and v 0 2 can be uniquely determined if the controls θ 1 and θ 2 are known. As u ε 1,0 and u 0 2,0 are independent of the virtual controls (θ 1 , θ 2 ), they can be computed beforehand.
The well posedness of the optimization problem is proved following Lions [22]. The key point consists in proving that the cost function induces a norm over U = (U D 1 , U D 2 ). One consider then the completion of U (still denoted by U) with respect to the cost induced norm, and the minimization problem admits a unique solution (θ 1 , θ 2 ) ∈ U satisfying the Euler-Lagrange formulation for all (μ 1 , μ 2 ) ∈ U, and where O is either ω 0 or Γ 1 ∪ Γ 2 . While the optimization problem with the cost function of case 1 has been considered in [7], we prove that the optimal controls problem is well posed with the cost function of case 2. For the homogenization theory (H-convergence), we consider a family of problems (1) indexed by ε. In what follows, we will often assume ε ≤ ε 0 , where ε 0 is a parameter used in a strong Cauchy-Schwarz inequality (see Lemma 6.3). We assume that θ i ∈ U D i and hence u i (θ i ) is in H 1 (ω i ), for i = 1, 2.

Minimization over Γ 1 ∪ Γ 2
As a fist step, we write the cost in terms of the state variables v ε 1 and v 0 2 , We set and show that π induce a norm over U.

Lemma 2.1
The bilinear form π is a scalar product over U.
Assume now that π(θ 1 , θ 2 ) = 0. It holds that and v 0 2 are H 1 functions on ω 1 and ω 2 , respectively, we obtain We now use H-convergence on the tensor a ε 1 , to obtain an homogenized tensor a 0 1 in ω 1 . It holds that v ε 1 converges weakly in H 1 toward v 0 1 the homogenized solution of Using the compact embedding of H 1 in L 2 , the solution v ε 1 converges strongly, up to a subsequence, toward v 0 1 in L 2 , and it holds that v 0 where C s < 1 is the strong Cauchy-Schwarz constant. The tensors a 0 2 and a 0 1 are equal in the overlapping region ω 0 , due to the locality of H-convergence, and one can bound the H 1 norm over ω 0 by the H 1/2 norm over its boundary Collecting the results leads to v 0 1 = 0 and v 0 2 = 0 a.e. in ω 0 . Further from the Caccioppoli inequality, Theorem 6.1, it holds v 0 1 = 0 a.e. in ω 1 and v 0 2 = 0 a.e. in ω 2 . We can conclude that θ i = 0 a.e. in Γ i , i = 1, 2, by using the trace inequality; i.e., The norm induced from the scalar product π is given by where O is either ω 0 or Γ 1 ∪ Γ 2 .

A priori error analysis
Let u ε be the solution of the heterogeneous problem (1), and let us derive a priori error bounds between u ε and the solution of the couplinḡ where u rec 2 is the reconstructed homogeneous solution u 0 2 , given by problem CITE, with periodic correctors, and ω + is a subdomain of Ω such that ω ⊆ ω + ⊆ ω 1 . The reconstructed homogeneous solution needs to be considered here instead of u 0 2 , as the later is a good approximation of the heterogeneous solution u ε only in the L 2 norm, but fails in the H 1 norm. The term u rec 2 is given by where u 0 2 = u 0 2 (θ 2 ).
We consider the cost function of case 2, and we refer to [7] for an analysis of case 1. For ε fixed, we consider a homogenized problem on ω 2 with boundary conditions on Γ 2 given by the trace of the heterogeneous solution u ε . The problem reads: find u 0 the solution of where γ 2 : H 1 (ω 2 ) → H 1/2 (Γ 2 ) denotes the trace operator on Γ 2 . Similarly, we define the trace operator γ 1 on Γ 1 . Assuming that the tensor a ε 2 is periodic in the fast variable, i.e., with periodic boundary conditions, and where (e i ) d i=1 denotes the canonical basis of R d . Assuming sufficient regularity on u 0 and on χ j , it can be proved that where the constant is independent of ε. For proofs, we refer to [9,21,24].
Estimates for the fine solution Let us define an operator P :

It can be split into
where the state variables v ε 1 and v 0 2 are solutions of (4) for i = 1, 2, respectively, and where U 0 is given by For the cost function of case 1, it has been shown in [7] that the operator Q is bounded in the operator norm, i.e., Here we assume that Q is bounded for the norm in U induced by the scalar product (6) for the cost function of case 2.
Theorem 2.2 Let u ε be solution of (1) andū be given by (7). Assume that u 0 and χ j are smooth enough for (10) to hold, and that Q ≤ C. Then we have where the constant C depends on the constant of the Caccioppoli inequality, the bound Q , and the trace constants associated with the trace operators γ 1 and γ 2 on Γ 1 and Γ 2 , respectively.
To prove Theorem 2.2, we note that the difference u ε −ū is a ε 1 -harmonic in ω 1 ; thus, Caccioppoli inequality (see Theorem 6.1) can be applied, Assuming that Q is bounded and using Lemmas 2.3 and 6.2, we can show that is bounded by Cε, which concludes the proof of Theorem 2.2.

Lemma 2.3
Let u ε and u 0 solve (1) and (9), respectively, and let (θ 1 , θ 2 ) ∈ U be the optimal virtual controls. Then Proof From the definition, it holds We look at the numerator. As the pair (θ 1 , θ 2 ) minimizes the cost function J , the Euler-Lagrange formulation (5) holds and The result follows.
The next lemma gives an upper bound to the norm in Lemma 2.3.

Lemma 2.4
Let u ε and u 0 be solution of (1) and (9), respectively. Assume that u 0 and χ j have enough regularity for (10) to hold. Then where the constant C is independent of ε.
Proof It holds Using the continuity of the traces, the first term can be bounded by This prove the result.
The proof of Theorem 2.2 follows from (11) and Lemmas 2.3 and 2.4.

Estimates for the coarse solution
The a priori error estimates to the coarse-scale solver follow from [7, Theorem 3.6] using Lemma 2.4. We skip the details.
where the constant C is independent of ε, but depends on τ , τ + , and the ellipticity constants of a ε 2 .

Fully discrete coupling method
In this section, we describe the fully discrete overlapping coupling method and perform an a priori error analysis. The fine-scale solver requires a triangulation of sizeh sufficiently small to resolve the multiscale nature of the tensor. In contrast, the coarse-scale solver on ω 2 takes full advantage of the scale separation and allows for a mesh size larger than the fine scale. We use the FEM in ω 1 and the FE-HMM in ω 2 . As the finite elements of the fine and coarse meshes in ω 0 are different, an interpolation between the two meshes should be considered. One can also chose to use the same finite elements in the overlap, leading to a discontinuity at Γ 1 in the mesh over ω 2 . In that latter situation, the discontinuous Galerkin FE-HMM [4] should be used instead of the FE-HMM.
In what follows, we consider for simplicity the problem (1) with homogeneous Dirichlet boundary conditions, i.e., we set g D = 0 and Γ N = ∅. Further, we assume that ε is small enough so that we can use the strong Cauchy-Schwarz lemma (Lemma 6.3) and its discrete version (Lemma 6.5) hold.
Numerical method for the fine-scale problem Let Th be a partition of ω 1 , in simplicial or quadrilateral elements, with mesh sizeh ε whereh = max K ∈Th h K , and h K is the diameter of the element K . In addition, we suppose that the family of partitions {Th} is admissible and shape regular [11], i.e., (T1) admissible: ω 1 = ∪ K ∈T h K and the intersection of two elements is either empty, a vertex, or a common face; (T2) shape regular: there exists σ > 0 such that h K /ρ K ≤ σ , for all K ∈ Th and for all Th ∈ {Th}, where ρ K is the diameter of the largest circle contained in the element K .
For each partition Th of the family {Th}, we define a FE space in ω 1 where R p is the space P p of polynomials of degree at most p on K if K is a triangle, and the space Q p of polynomials of degree at most p in each variable if K is a rectangle. Further, V p 0 (ω 1 , Th) denotes the space of functions in V p D (ω 1 , Th) that vanish on ∂ω 1 . Let u 1,h be the numerical approximation of u ε 1 satisfying problem (3) Th) is obtained by the optimization method and u 1,0,h ∈ V p 0 (ω 1 , Th) satisfies where F 1 is given by Thanks to the Poincaré inequality, the coercivity and boundedness of the bilinear form B 1 can be proved; the existence and uniqueness of u 1,0,h follow. Quadrature formula A macroscopic quadrature formula is given by the pair {x j,K , ω j,K } of quadrature nodes x j,K and weights ω j,K , for j = 1, . . . , J. The sampling domain of size δ around each quadrature point is denoted by We assume that the quadrature formula verifies the necessary assumptions to guarantee that the standard error estimates for a FEM hold [11]. The numerically homogenized tensor a 0,h 2 (x j,K ) is obtained using numerical solutions of micro-problems defined in K δ j . In each sampling domain, we consider a mesh T h in simplicial or quadrilateral elements K with mesh size h = max K ∈T h h K satisfying h < ε. The micro-FE space is where the space W (K δ j ) depends on the boundary conditions in the micro-problems; W (K δ j ) = H 1 0 (K δ j ) for Dirichlet coupling, or W (K δ j ) = W 1 per (K δ j ) for periodic coupling.
The discrete micro-problems read: find ψ i,h The numerically homogenized tensor can be computed by ). We define a macro-bilinear form B 2,H (·, ·) over for all w 1,h ∈ V p 0 (ω 1 , Th) and w 2,H ∈ V p 0 (ω 2 , T H ). We introduce discrete Lagrange multipliers for each of the constraints and obtain a discrete optimality system: find where the unknown vector U is given by U = (v 1,h , v 2,H , λ 1,h , λ 2,H ) , and Fully discrete error estimates The coupling solution, denoted byūh H , is defined as where u rec 2,H (θ 2,H ) is a fine-scale approximation obtained from the coarse-scale solution u 2,H (θ 2,H ) using a post-processing procedure in the following way. We assume that the tensor a ε 2 is Y -periodic in y, and we restrict the FE spaces to piecewise FE spaces. Periodic coupling is then used with sampling domains K ε of size ε. The reconstructed solution u rec 2,H (θ 2,H ) is given by where ψ j,h K ε are the micro-solutions of (12) in the sampling domain K ε . As the numerical solutions might be discontinuous in ω 2 , we consider a broken H 1 semi-norm, We next state our main convergence result for the optimization-based numerical solution.
Let u H ∈ V 1 0 (ω 2 , T H ) be the FE-HMM approximation of the homogenized solution u 0 .
Theorem 3.1 (A priori error analysis in ω + ) Let ε 0 be given by the strong Cauchy-Schwarz lemma, Lemma 6.3, and consider ε ≤ ε 0 . Let u ε and u 0 be the exact solutions of problems (1) and (9), respectively, andūh H be the numerical solution of the coupling (18). Further, let u H ∈ V 1 0 (ω 2 , T H ) be the FE-HMM approximation of u 0 . Assume u ε ∈ H s+1 (Ω), with s ≤ 1, u 0 ∈ H 2 (ω 2 ), and assume that (10) holds, then The analysis of the error e HMM,L 2 is by now standard for the FE-HMM. One decompose the error into [2] e HMM, where u 0 H is a FEM approximation of u 0 with numerical quadrature andū H is a semidiscrete FE-HMM approximation of u 0 , where the micro-functions are in the exact Sobolev space W (K δ ). Under suitable regularity assumption [12], we have where the constant C is independent of ε,h, H, and h.
Next, following [1,2] we can bound the micro-and modeling errors. If we assume the following regularity on ψ i K ε ∈ W (K δ ), the non-discretized micro-solutions of problem (12), we obtain a bound on the micro-error where the constant C is independent of ε,h, H, and h (we recall that for the reconstruction we use periodic boundary conditions in the micro-problems (12) over sampling domains are of size δ = ε). If we collocate (i.e., freeze) the slow variable x to the quadrature point x K in the tensor a ε 2 , i.e., we consider a ε 2 (x K , x/ε) in the macro-and micro-bilinear forms, we obtain an optimal modeling error assuming that δ/ε ∈ N >0 . Without collocation, the modeling error becomes where the constant C is independent of ε,h, H, and h. Remark 3.2 When δ/ε / ∈ N and δ > ε, Dirichlet boundary conditions are used instead of the periodic conditions in the micro-problems (12), and the modeling error becomes where the constants are independent of ε,δ,h, H, and h. Remark 3.3 Higher-order FE macro-and micro-spaces can also be considered, and we refer to [2,3] for details.
Next, we state an error estimates in the coarse-scale region for the optimization-based numerical solution with correctors.

Theorem 3.4 (Error estimates in Ω\ω + ) Let u ε be the exact solution of problem (1)
andūh H be the numerical solution of the coupling (18). Let a ε 2 (x) = a 2 (x, x/ε), where a 2 (x, y) is Y -periodic in y and satisfies a 2 (x, y) ∈ C(ω 2 ; where the constants are independent of H,h, h, and ε.
Proof It follows the lines of [7,Theorem 4.4], where DG-FE-HMM is replaced by FE-HMM.
Remark 3. 5 We note that the above theorem is also valid when using discontinuous Galerkin macro-solver (i.e., the DG-FE-HMM [4]). This has been studied in [7] for the cost function of case case 1. A similar proof applies for the cost function of case 2.

Numerical experiments
In this section, we give three numerical experiments that can be seen as a complement of the ones carried in [7], where we focused on a minimization in L 2 (ω 0 ), with interior subdomains and matching grids in the overlap ω 0 . In the first experiment, we still consider the minimization over L 2 (ω 0 ) and compare matching and non-matching meshes. The second experiment illustrates the coupling with the cost function of case 2 over Γ 1 ∪ Γ 2 and comparisons with the cost function of case 1 over ω 0 . In the last example, we combine non-matching grids and a minimization over the boundary. We observe several order of magnitude of saving in computational cost when compared to the method proposed in [7].
In the experiments, we will use a tensor a ε which is highly heterogeneous non-periodic and oscillate at several non-separated scales in ω, and which has scale separation in ω 2 , with a locally periodic structure.

Comparison of matching and non-matching grids on the overlap
Experiment 1 For this experiment, we use the cost function of case 1 Using FEM and FE-HMM in ω 1 and ω 2 , respectively, leads to two main restrictions: the mesh size in ω 1 should be smaller than the fine scale, whereas the mesh size in ω 2 can be larger than the fine scales, in order to take full advantage of the FE-HMM. Since both methods are defined in ω 0 , we can chose to have the same FE in both meshes on the overlap, or one can impose two different meshes. With the first choice, no interpolations must be considered between Th and T H over ω 0 , but T H is composed of FE with mesh size as small as the fine scales. In that situation, DG-FE-HMM is chosen instead of FE-HMM due to the discontinuity at the interface Γ 1 . The second choice requires interpolation between the meshes in ω 0 , but T H is not restricted by the size of the fine mesh Th. We show that both cases give similar convergence rates, but the computational cost is significantly reduced in the second case.

Let us consider a Dirichlet elliptic boundary value in
with f ≡ 1 and a ε is given by . Let H = 1/8, ε = 1/10, and a micro-mesh size h = ε/L, so that the micro-error is negligible. We initialize the fine mesh toh = 1/16. We use uniform simplicial meshes in ω 1 and ω 2 and assume that Th is obtained from T H using a uniform refinement in ω 0 . This allows simplification in the interpolation between the two meshes in the overlap. We couple the FEM over ω 1 with the mesh Th(ω 1 ) with the FE-HMM over ω 2 with mesh T H (ω 2 ) and compare it with a coupling between FEM over Th(ω 1 ) with DG-FE-HMM over a mesh composed of coarse FE from T H (ω 2 \ ω 0 ) with small FE from the fine mesh Th(ω 0 ). The reference fine-scale solution is computed on a very fine mesh, and we compare the two numerical solutions with the reference one. After three iterations, we plot the numerical approximations of the fine-scale solution u ε 1 and coarse-scale solution u 0 2 (in transparent), for a coupling with minimization of the cost function of case 1 with non-matching grid (Fig. 2a) and with matching grids (Fig. 2b). A zoom of the coarse-scale solutions in the overlap region ω 0 is shown in Fig. 2c for the coupling with non-matching grids and in Fig.  2d with matching grids, where the coupling is performed with the cost function of case 1 (as the fine meshes become too dense after three iterations, we plot for the zoom the solution after one iteration to better visualize the difference in the meshes).
We refine either only in ω 1 for the fine-scale solver (non-matching grids) or in addition in ω 0 for the coarse-scale solver (matching grids). We set δ = ε for the sampling domains and consider a micro-mesh size h = ε/L, so that the micro-error is negligible. Figure  3a shows the H 1 norm in ω with non-matching grids (bullet) and with matching grids (diamond); we see that the errors are similar. We also measured the times, using MATLAB timer, to compute the numerical solutions. We see in Fig. 3b that using non-matching grids is faster as the number of micro-problems that have to be computed with the coarse solver, is smaller and fixed, whereas it increases when matching grids are used, causing a significant time overhead.
The rate of convergence in ω is influenced by H and ε, and whenh is refined, we expect a saturation, depending on H and ε, in the convergence. Let ε = 1/20 and initialize the fine mesh toh = 1/64. We set H = 1/8, 1/16, and 1/32 and refineh in each iteration. In Fig. 4, we plot the H 1 norm between the reference and numerical solutions w.r.t the mesh size in ω. We see indeed that the error saturates at a threshold value that depends on H.  The tensor a ε in ω 2 has scale separation and is Y -periodic in the fast variable, and the homogeneous tensor a 0 2 can be explicitly derived as Let ω 1 = [0, 1/2] × y and ω = [0, 1/4] × y, with y ∈ [0, 1]. An illustration of a numerical solution is given in Fig. 6a. At first, we consider the cost of case 1, Let ε = 1/50, and h/ε = 1/L be small enough to neglect the micro-error. We initialize the fine mesh toh = 1/128. For different macro-mesh sizes H = 1/8, 1/16, 1/32 and 1/64, we refineh and monitor the convergence rates between the numerical solution of the coupling and the reference solution. In Fig. 5a, the H 1 norm is displayed for H = 1/8 (dots), H = 1/16 (dashes-dots), H = 1/32 (dashes) and H = 1/64 (full lines). One can see that the error saturates at a value depending on the macro-mesh size H. Now, we compare the costs of case 1 over ω 0 with the cost of case 2 over Γ 1 ∪ Γ 2 . We fix ε = 1/10, H = 1/16, and h = ε/L small enough in order to neglect the micro-error. We initialize the fine mesh toh = 1/32 and refine the mesh only in ω 1 . The numerical approximations of u ε 1 and u 0 2 are shown in Fig. 6a, for the cost of case 1 over ω 0 , and in Fig.  6b, for the cost of case 2 over Γ 1 ∪ Γ 2 . The H 1 and L 2 errors between u H and a reference solution in ω 0 are shown in Fig. 6c, d, respectively. Computational times are compared as well in Fig. 7, for the cost over ω 0 (diamonds) and the cost over Γ 1 ∪ Γ 2 (bullets). As the number of degrees of freedom of the saddle point problem (17) is reduced when minimizing over the boundaries Γ 1 ∪ Γ 2 , we see that the coupling over ω 0 is more costly than the coupling over Γ 1 ∪ Γ 2 . Considering an interpolation between the two meshes in the interface ω 0 gives similar results as, due to the periodicity of a ε 2 , we need only to resolve one cell problem to compute the homogenized tensor a 0 2 . We next vary the size of the overlap ω 0 and consider ω 1 = [0, 1/4 + mH] × y, for m = 1, 4, 8, where H = 1/32 is the coarse mesh size, and initializeh = 1/64. We minimize over the overlap ω 0 . We observe that both couplings are influenced by the size of τ = dist(Γ 1 ∪ Γ 2 ) and this is shown in the H 1 errors in Fig. 8. The rates deteriorate when τ goes to zero.
Minimization with interface control on non-matching grids For the last experiment, we combine the two previous effects. The fastest coupling is obtained by performing by considering the minimization with of the cost of case 2 with interpolation of the two meshes in the overlap, whereas the slowest coupling is obtained by the minimization with the cost function of case 1 using identical meshes in the overlap. −div (a ε ∇u ε ) = f in Ω, with f ≡ 1 and a ε is given by We set H = 1/16 and ε = 1/10. We initializeh = 1/32. In Fig. 9a, we see the H 1 error for the two settings is similar, whereas the computational cost using minimization over the overlap and non-matching grid in ω 0 dramatically decreases (see Fig. 9b).

Conclusion
In this paper, we aimed at reducing the computational cost of our optimization-based method, proposed in [6,7]. Our focus was to reduce the total number of degrees of freedom of our method by investigating two strategies; i.e., (i) to use a cost function over the boundary of the overlapping region; (ii) to consider independent meshes in the overlap and an interpolation procedure between the fine and coarse meshes in the overlapping region.
We derived a fully discrete a priori error analysis of the coupling method based on the new cost function and the heterogenous meshing strategy. Numerical experiments have been carried out to assess the difference between the two costs functions and the two meshing strategies, with and without matching grids in the overlap. In all numerical experiments, we observe orders of magnitude of saving in the computational cost when we compare the numerical settings used in [7] with a combination of the strategies (i) and (ii).
We note that elliptic problems with a non-null right-hand side and problems where ∂ω ∩ Γ = ∅, can also be considered and we refer to [19] for details. We give next a bound of the L 2 norm over ω by the L 2 norm over the overlap ω 0 .

Lemma 6.2
Let v ε 1 and v 0 2 be solutions of (4), for i = 1, 2,, respectively. The following bounds hold: where τ is the width of the overlap and C is a constant depending on α, β, and the Poincaré constant is associated with ω 1 and ω 2 , respectively.
In the next lemma, we state a strong version of the Cauchy-Schwarz inequality and refer to [7] for the proof. Let us recall the problems for the state variables: where a 1 = a ε 1 and a 2 = a 0 2 .

Lemma 6.3 (Strong Cauchy-Schwarz)
Let v ε 1 ∈ H 1 D (ω 1 ) and v 0 2 ∈ H 1 D (ω 2 ) be solutions of (19) for i = 1, 2. Then, there exist an ε 0 > 0 and a positive constant C s < 1 such that for all ε ≤ ε 0 , it holds Discrete versions of the Caccioppoli and the strong Cauchy-Schwarz inequalities are stated below. v h ∈ V p (ω 1 , T h ) solution of where the constant C is independent of h.
We now give the discrete strong Cauchy-Schwarz inequality, and to simplify the notations, we omit the ε dependency in v 1 .