Goal-oriented error estimation and adaptivity in MsFEM computations

TL;DR

This paper presents a goal-oriented adaptive error estimation strategy for MsFEM that combines Constitutive Relation Error with an adjoint problem to improve accuracy and efficiency in multiscale simulations.

Contribution

It introduces a novel goal-oriented approach coupling CRE with adjoint problems for adaptive MsFEM error control and proposes a local enrichment technique for better bounds.

Findings

Error estimates are highly accurate and reliable.

The adaptive procedure significantly reduces computational costs.

The method performs well across various numerical test cases.

Abstract

We introduce a goal-oriented strategy for multiscale computations performed using the Multiscale Finite Element Method (MsFEM). In a previous work, we have shown how to use, in the MsFEM framework, the concept of Constitutive Relation Error (CRE) to obtain a guaranteed and fully computable a posteriori error estimate in the energy norm (as well as error indicators on various error sources). Here, the CRE concept is coupled with the solution of an adjoint problem to control the error and drive an adaptive procedure with respect to a given output of interest. Furthermore, a local and non-intrusive enrichment technique is proposed to enhance the accuracy of error bounds. The overall strategy, which is fully automatic and robust, enables to reach an appropriate trade-off between certified reliability and computational cost in the MsFEM context. The performances of the proposed method are investigated on several illustrative numerical test cases. In particular, the error estimation is observed to be very accurate, yielding a very efficient adaptive procedure.

Equations

$$-\nabla \cdot ({\mathbb{A}}^{\epsilon }\nabla u^{\epsilon })=f\text{in Ω},u^{\epsilon }=0\text{on ΓD},{\mathbb{A}}^{\epsilon }\nabla u^{\epsilon }\cdot \mathbf{n}=g\text{on ΓN}.$$

$$\forall \epsilon \ge 0,\forall \mathbf{\xi }\in {\mathbb{R}}^d,\alpha |\mathbf{\xi }{|}^2\le {\mathbb{A}}^{\epsilon }(\mathbf{x})\mathbf{\xi }\cdot \mathbf{\xi }\le \beta |\mathbf{\xi }{|}^2\text{a.e. in Ω.}$$

$$\forall v\in V,B^{\epsilon }(u^{\epsilon },v)=F(v)$$

$$B^{\epsilon }(u,v)={\int }_{\Omega }{\mathbb{A}}^{\epsilon }\nabla u\cdot \nabla v,F(v)={\int }_{\Omega }fv+{\int }_{{\Gamma }_N}gv.$$

$$\forall v\in V_H^0,B^{\epsilon }(u_H,v)=F(v).$$

$$\nabla \cdot ({\mathbb{A}}^{\epsilon }\nabla {\varphi }_i^{\epsilon })=0\text{in K},$$

$$\forall v\in V_H^{\epsilon },B^{\epsilon }(u_H^{\epsilon },v)=F(v).$$

$$\nabla \cdot ({\mathbb{A}}^{\epsilon }\nabla {\varphi }_i^{\epsilon })=0\text{in K},{\varphi }_i^{\epsilon }={\varphi }_i^0\text{on ∂K}.$$

$$\Vert u^{\epsilon }-u_H^{\epsilon }{\Vert }_1\le C\left(H+\sqrt{\frac{\epsilon }{H}}+\sqrt{\epsilon }\right),$$

$$\nabla \cdot ({\mathbb{A}}^{\epsilon }\nabla {\psi }_i^{\epsilon })=0\text{in SK},{\psi }_i^{\epsilon }\text{is affine on ∂SK},{\psi }_i^{\epsilon }({\mathbf{s}}_j)={\delta }_{ij},$$

$$\forall v\in V_H^{\epsilon },B_H^{\epsilon }(u_H^{\epsilon },v)=F(v).$$

$$\Vert u^{\epsilon }-u_H^{\epsilon }{\Vert }_{1,H}\le C\left(H+\frac{\epsilon }{H}+\sqrt{\epsilon }\right),$$

$$u_H^{\epsilon }\in V_H^{\epsilon }$$

$$\nabla \cdot ({\mathbb{A}}^{\epsilon }\nabla w_j^{\epsilon ,K})=0\text{in K},w_j^{\epsilon ,K}=\mathbf{x}\cdot {\mathbf{e}}_j\text{on ∂K}.$$

$${\varphi }_i^{\epsilon }={\varphi }_i^0\circ {\mathbf{w}}^{\epsilon }.$$

$$\Vert u^{\epsilon }-u_H^{\epsilon }{\Vert }_1\le C\left(H^k+\sqrt{\frac{\epsilon }{H}}+\sqrt{\epsilon }\right),$$

$$\mathcal{S}=\left\{\mathbf{p}\in \mathcal{H}(div,\Omega );\forall v\in V,{\int }_{\Omega }\mathbf{p}\cdot \nabla v={\int }_{\Omega }fv+{\int }_{{\Gamma }_N}gv\right\}$$

$$E_{CRE}(\hat{u},\hat{\mathbf{p}})={\left({\int }_{\Omega }({\mathbb{A}}^{\epsilon }{)}^{-1}(\hat{\mathbf{p}}-{\mathbb{A}}^{\epsilon }\nabla \hat{u})\cdot (\hat{\mathbf{p}}-{\mathbb{A}}^{\epsilon }\nabla \hat{u})\right)}^{1/2}={\left|\left|\left|\hat{\mathbf{p}}-{\mathbb{A}}^{\epsilon }\nabla \hat{u}\right|\right|\right|}_{\mathcal{F}},$$

$$\forall (\hat{u},\hat{\mathbf{p}})\in V\times \mathcal{S},{\left(E_{CRE}(\hat{u},\hat{\mathbf{p}})\right)}^2={\left|\left|\left|{\mathbf{q}}^{\epsilon }-\hat{\mathbf{p}}\right|\right|\right|}_{\mathcal{F}}^2+{\left|\left|\left|u^{\epsilon }-\hat{u}\right|\right|\right|}^2,$$

$$E_{CRE}(\hat{u},\hat{\mathbf{p}})=2{\left|\left|\left|{\mathbf{q}}^{\epsilon }-{\hat{\mathbf{p}}}^{\ast }\right|\right|\right|}_{\mathcal{F}}\text{where}{\hat{\mathbf{p}}}^{\ast }=\frac{1}{2}(\hat{\mathbf{p}}+{\mathbb{A}}^{\epsilon }\nabla \hat{u}).$$

$$\left|\left|\left|u^{\epsilon }-u_H^{\epsilon }\right|\right|\right|\le E_{CRE}(u_H^{\epsilon },\hat{\mathbf{p}}).$$

$${\left|\left|\left|u^{\epsilon }-u_H^{\epsilon }\right|\right|\right|}_H\le E_{CRE}({\hat{u}}_H^{\epsilon },\hat{\mathbf{p}})+{\left|\left|\left|{\hat{u}}_H^{\epsilon }-u_H^{\epsilon }\right|\right|\right|}_H,$$

$$\forall K,{\int }_{K}f+{\int }_{\partial K}{\hat{g}}_K=0,$$

$${\hat{g}}_{K|\Gamma }=\frac{1}{2}\left({\mathbf{q}}_{H|K}^{\epsilon }+{\mathbf{q}}_{H|K^{\prime }}^{\epsilon }\right)\cdot {\mathbf{n}}_K+\sum _i{\alpha }_i{\varphi }_{i|\Gamma }^{\ast },$$

$$\forall v\in {\mathcal{H}}^1(K),{\int }_K{\hat{\mathbf{q}}}_{H|K}^{\epsilon }\cdot \nabla v={\int }_{K}fv+{\int }_{\partial K}{\hat{g}}_{K}v.$$

$$\forall v\in {\mathcal{H}}^1(K),{\int }_K{\mathbb{A}}^{\epsilon }\nabla w^{\epsilon }\cdot \nabla v={\int }_{K}fv+{\int }_{\partial K}{\hat{g}}_{K}v$$

$$\text{Find uε∈V such that, for any v∈V,}{B}^{\epsilon \ast }({\tilde{u}}^{\epsilon },v)=Q(v),$$

$$\tilde{\mathcal{S}}=\left\{\tilde{\mathbf{p}}\in \mathcal{H}(div,\Omega );\forall v\in V,{\int }_{\Omega }\tilde{\mathbf{p}}\cdot \nabla v=Q(v)\right\}.$$

$$Q(v)={\int }_{\Omega }\left({\tilde{f}}_{Q}v+{\tilde{\mathbf{q}}}_Q\cdot \nabla v\right)+{\int }_{{\Gamma }_N}{\tilde{g}}_{Q}v,$$

$$Q(u^{\epsilon })=B^{\epsilon \ast }({\tilde{u}}^{\epsilon },u^{\epsilon })=B^{\epsilon }(u^{\epsilon },{\tilde{u}}^{\epsilon })=F({\tilde{u}}^{\epsilon }).$$

Full text

1 Introduction and objectives

Highly heterogeneous materials are ubiquitous in Material Sciences and Computational Mechanics. Examples include porous or perforated media, composite materials, …Developing multiscale numerical methods to simulate the (thermal, mechanical, …) response of such materials is an increasingly active research field. As an alternative to solving a full fine-scale problem (the cost of which is usually prohibitively expensive), multiscale modeling aims at linking the different scales by appropriately capturing the effect of the smaller scales on the larger scales. This enables an accurate description of complex physical phenomena and/or the prediction of macroscopic properties at a limited computational cost. In the recent years, significant progress has been witnessed in the design of multiscale methods for elliptic PDEs that take into account scale separation between macroscopic and microscopic features. We can distinguish two families of multiscale methods. The first family includes approaches relying on multilevel finite elements, such as HMM or FE2 methods [22, 26]. They perform a numerical homogenization of the medium in order to find local effective material properties, which are next inserted in a global macroscopic problem which can be handled with a reduced number of degrees of freedom. Such approaches can be considered as numerical generalizations of classical homogenization approaches (see e.g. [9]) to non-periodic media. The second family includes approaches based on superposition of scales, such as the GFEM [64], the MsFEM [33], the LOD [51] and the VMS [36] methods, to name but a few. In this work, we focus on the Multiscale Finite Element Method (MsFEM) which has been the object of numerous research works (see [24] for an overview). This method is based on the classical FEM framework, defining an approximate solution in a finite dimension space related to a macroscopic mesh, but uses specific basis functions which encode fine-scale details. The basis functions are obtained, in a preliminary offline stage, by solving fine-scale equations over local subdomains. The MsFEM approximation is next computed during an inexpensive online stage, and provides a numerical approximation that describes details of the exact solution at different scales, including the small ones.

A recurring concern when using multiscale methods is the control of the accuracy of the numerical solution. In particular, it is critical to understand, quantify and control how errors generated by different sources and at different scales affect the overall error, in order to build relevant adaptive strategies and focus the computational efforts where they are required. Following the large literature on a posteriori error estimation in the monoscale framework of the Finite Element Method (FEM) (see e.g. [67, 3] for classical textbooks and [18] for a recent review), several verification tools have been developed for multiscale methods [65, 44, 45, 52, 1, 37, 2, 31, 20]. However, many open questions remain, such as quantitative assessment of error propagation across scales, relevant adaptive strategies, …Recently, we have proposed a robust a posteriori error estimation tool for MsFEM computations [19]. It uses the concept of Constitutive Relation Error (CRE), which has been developed in the Computational Mechanics community for more than thirty years [42]. It enables both (i) to derive fully computable (i.e. without any unknown multiplicative constant) and guaranteed bounds on the global error measured in the energy norm, which may further be used as stopping criteria in adaptive algorithms, from an extension of the CRE to the multiscale context, and (ii) to assess the various MsFEM error sources and drive an adaptive algorithm from local error maps in order to reach a prescribed error tolerance.

In the present article, we extend the methodology introduced in [19] to a posteriori goal-oriented error estimation, namely the control of the error on specific scalar-valued outputs of interest for problems discretized with MsFEM. In this specific context, material heterogeneities may indeed require standard adaptive methods to add degrees of freedom in regions with limited influence on the output of interest, or not to add some in regions which significantly affect this output. Consequently, alternative a posteriori error estimation and adaptivity methods need to be introduced, taking into account the specific quantities of interest of the problem. For that purpose, we introduce classical duality arguments along with an adjoint problem [6, 59, 63, 21, 62, 53, 29, 16], which enable to target the error estimation and adaptivity procedure towards the approximation of the considered quantity of interest. We also use strategies to bound the error that are similar to those employed for monoscale applications with CRE [38, 17, 39, 58, 68, 43]. In particular, we introduce a non-intrusive approach for the approximate solution of the adjoint solution (using so-called handbook functions), which enables to obtain both guaranteed and sharp error bounds. We also define local indicators on various error sources (coarse mesh size, size of the local fine meshes, …) in order to drive the adaptive strategy, and automatically and iteratively select optimal MsFEM parameters to be used locally for reducing the error in the goal functional. In order to comply with the MsFEM philosophy and the multi-query paradigm, all additional costly computations associated with the proposed procedures are performed in the offline stage. The overall verification tool enables to perform MsFEM computations in which: (i) a given precision on specified outputs of interest is attained at the end of the simulation; (ii) the computational work which is needed is as small as possible. The proposed approach, devoted to the Multiscale Finite Element Method and which fulfills the above two objectives, is an alternative to the goal-oriented adaptive strategy recently proposed in [20] within the Generalized Multiscale Finite Element framework, which uses residual-based techniques and which mainly focuses on goal-oriented adaptivity alone. We also mention the approach [44] proposed in the VMS framework, and the approach [50] proposed in the HMM framework where the overall error is also split between different contributions each related to a discretization parameter.

The performances of our approach are evaluated on several numerical examples. It is interesting to note that it leads in practice to MsFEM computations with guaranteed prescribed accuracy, and in which fine-scale features of the solution are required (and captured using fine discretization parameters) only in very localized parts of the physical domain. A limited accuracy of the fine scale computations is thus sufficient. Consequently, it naturally and automatically shares some resemblance with alternative multiscale computing strategies based on local enrichments or couplings between fine-scale and coarser-scale (upscaled) models [48, 8, 15, 28, 49], these strategies being associated to modeling error estimation [57, 56, 66, 54, 61, 4].

This article is organized as follows. In Section 2, we introduce the model problem we consider along with basic notions on the MsFEM discretization approach. The CRE concept applied to the MsFEM framework is recalled in Section 3. In Section 4, we extend this concept to goal-oriented error estimation within MsFEM computations. This is the core of this work. An associated adaptive strategy, enabling to control the accuracy of outputs of interest, is introduced in Section 5. We discuss several illustrative numerical experiments in Section 6. Eventually, conclusions and perspectives are collected in Section 7.

2 Model problem and MsFEM approximation

To start with, we describe our model problem in Section 2.1. In Section 2.2, we next recall some basic elements on several variants of MsFEM: the initial conforming variant, the oversampling variant (see Section 2.2.3) and the higher-order variant introduced by G. Allaire and R. Brizzi (see Section 2.2.4). The reader familiar with MsFEM may easily skip Section 2.2 and directly proceed to the a posteriori error estimation strategy described in Section 3.

2.1 Model problem

We consider mathematical models described by partial differential equations (PDEs) and characterized by multiscale material heterogeneities. For the sake of simplicity, we focus on the following elliptic problem posed on some open and bounded domain $\Omega\subset\mathbb{R}^{d}$:

[TABLE]

The boundary $\Gamma=\partial\Omega$ is split into a subset $\Gamma_{D}\subset\partial\Omega$ (with measure $|\Gamma_{D}|>0$) where homogeneous Dirichlet boundary conditions are applied, and a subset $\Gamma_{N}=\partial\Omega\setminus\Gamma_{D}\subset\partial\Omega$ where Neumann boundary conditions are applied. We assume that $f\in L^{2}(\Omega)$ and $g\in L^{2}(\Gamma_{N})$, corresponding to body and traction loadings respectively, are slowly varying functions. The quantity $\boldsymbol{q}^{\varepsilon}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}u^{\varepsilon}$ is the flux associated to $u^{\varepsilon}$ through the constitutive relation.

The second-order diffusion tensor $\mathbb{A}^{\varepsilon}\in(L^{\infty}(\Omega))^{d\times d}$ is a heterogeneous object, with rapid oscillations and possibly high contrast. The parameter $\varepsilon$ refers to the typical (presumably small) length scale of the heterogeneities. We assume that $\mathbb{A}^{\varepsilon}$ is a symmetric matrix that is uniformly elliptic and bounded, in the sense that there exist $\alpha>0$ and $\beta>0$ such that

[TABLE]

Introducing the functional space V=\{v\in\mathcal{H}^{1}(\Omega),\ \ \text{v=0$on$\Gamma_{D}}\}, the weak formulation of (1) consists in finding $u^{\varepsilon}\in V$ such that

[TABLE]

with

[TABLE]

The bilinear form $B^{\varepsilon}$ is symmetric, continuous and coercive on $V$. It thus defines an inner product and induces the energy norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\sqrt{B^{\varepsilon}(v,v)}$ on $V$. The $\mathcal{H}^{1}$-norm is denoted $\|\cdot\|_{1}$, and it is equivalent to the norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ on $V$.

The model problem (2) may be numerically solved using classical and direct numerical techniques, such as the finite element method (FEM). In this method, a conforming partition $\mathcal{T}_{H}$ of $\Omega$ is introduced ($H$ denotes the characteristic size of the elements of $\mathcal{T}_{H}$) and a finite dimensional space $V_{H}^{0}\subset V$ is constructed from Lagrange basis functions $\phi_{i}^{0}$ associated with each node $i$ of $\mathcal{T}_{H}$. The solution to (2) is thus approximated by $\displaystyle u_{H}=\sum_{i}u_{i}\,\phi_{i}^{0}\in V_{H}^{0}$, solution to

[TABLE]

A very fine mesh is necessary to represent the multiscale nature of the solution $u^{\varepsilon}$. This leads to choosing $H<\varepsilon$, and thus a prohibitively large number of degrees of freedom. The resulting method has a too high computational complexity.

As an alternative, multiscale model reduction techniques based on a scale separation assumption are typically used to solve (2) accurately and efficiently. On the one hand, homogenization strategies can be employed when the structure of $\mathbb{A}^{\varepsilon}$ is periodic (or locally periodic, or random and stationary, …; see [9]), or when one is only interested in macroscopic/upscaled properties of the model (1) (we refer e.g. to [22, 26]). On the other hand, multiscale numerical approaches have the capability to approximate, efficiently and for general structures of $\mathbb{A}^{\varepsilon}$, the details of $u^{\varepsilon}$ at all scales. In the following, we focus on the Multiscale Finite Element Method (MsFEM) [33, 24, 46], which is one of these latter approaches.

2.2 MsFEM approximation

2.2.1 General strategy

We define two partitions of $\Omega$: (i) a coarse grid partition $\mathcal{T}_{H}$, referred to as the macro mesh; (ii) a fine grid partition $\mathcal{T}_{h}$, subordinate to $\mathcal{T}_{H}$, referred to as the micro mesh (the mesh size $h$ is chosen smaller, and often much smaller, than $H$). The mesh $\mathcal{T}_{h}$ is in practice obtained by refining each coarse element of $\mathcal{T}_{H}$. It should be chosen so that oscillations of the data can be accurately captured on $\mathcal{T}_{h}$.

The main idea in the MsFEM approach is to construct a set $\{\phi_{i}^{\varepsilon}\}$ of multiscale basis functions, associated with each node $i$ of the macro mesh $\mathcal{T}_{H}$, which encode fine scale information and local heterogeneities of the solution space. These are constructed locally, in an offline phase, as the solutions, over each element $K$ of the coarse mesh (or over a local region $S_{K}\supset K$), to problems of the form

[TABLE]

complemented with boundary conditions discussed below. We notice that the local problems (3) are decoupled one from each other. They are thus well-suited to parallel computations. Furthermore, they do not depend on the loads $f$ and $g$, which is a main advantage in a multi-query context. These local problems are in practice solved numerically using the micro mesh $\mathcal{T}_{h}$.

Remark 1

A variant in the construction of multiscale basis functions, referred to as the Generalized Multiscale Finite Element Method (GMsFEM), was recently introduced in [23]. It consists in defining a larger number of basis functions over each element $K$ of $\mathcal{T}_{H}$, which are obtained using local spectral problems.

Once the multiscale basis functions are computed, a classical Galerkin approximation of (2) is performed on the finite dimensional space $V^{\varepsilon}_{H}=\text{Span}\{\phi_{i}^{\varepsilon}\}$. The solution to (2) is thus approximated by some $\displaystyle u^{\varepsilon}_{H}=\sum_{i}u_{i}\,\phi_{i}^{\varepsilon}\in V_{H}^{\varepsilon}$, defined as the solution (computed during the online phase) to the following global problem:

[TABLE]

We note that $V_{H}^{\varepsilon}$ and $V_{H}^{0}$ share the same dimension. The above problem hence involves a limited number of degrees of freedom. The small scale information incorporated in the basis functions $\phi_{i}^{\varepsilon}$ is thus brought to the large scales through couplings in the global stiffness matrix. The assembly of this matrix is inexpensive since it reuses local matrices computed and stored in the offline stage.

The MsFEM approach leads to a considerable decrease in memory requirements and CPU cost as an accurate approximation may be obtained for $H\gg\varepsilon$. Nevertheless, a posteriori error estimates and adaptive techniques are required to ensure a sufficiently high-fidelity approximation at a limited cost. These are addressed in Sections 3 and 4.

2.2.2 Conforming MsFEM

The choice of the boundary conditions for the local problems (3) (when constructing the multiscale basis functions) is of critical importance, since it can significantly affect the accuracy of the MsFEM solution. A first choice, leading to a conforming numerical discretization, is to impose an affine evolution of $\phi_{i}^{\varepsilon}$ along $\partial K$. Introducing the classical first-order and piecewise affine FE basis functions $\phi_{i}^{0}$, we thus consider the local problems

[TABLE]

In the regime of interest $H\geq\varepsilon$, assuming that $\mathbb{A}^{\varepsilon}$ has a periodic structure (i.e. $\mathbb{A}^{\varepsilon}=\mathbb{A}^{\rm per}(\cdot/\varepsilon)$ for some fixed matrix $\mathbb{A}^{\rm per}$) and that the local problems (5) are solved exactly (and under some technical regularity assumptions), the following a priori error estimate holds [34]:

[TABLE]

where the constant $C$ is independent of $\varepsilon$ and $H$. This estimate is derived using the periodic homogenization theory that provides (when $\varepsilon\ll 1$) a detailed description of $u^{\varepsilon}$ and $\phi_{i}^{\varepsilon}$ in terms of a two-scale expansion [9]. It is the sum of three terms:

the term in $O(\sqrt{\varepsilon})$ comes from the homogenization theory, and quantifies the difference between $u^{\varepsilon}$ and its two-scale expansion (this error is essentially located close to $\partial\Omega$; recall in particular that $u^{\varepsilon}$ satisfies homogeneous Dirichlet boundary conditions on $\Gamma_{D}$, whereas its two-scale expansion does not);

• 2.

the term in $O(\sqrt{\varepsilon/H})$ again comes from the homogenization theory, and quantifies the difference between any basis function $\phi_{i}^{\varepsilon}$ and its two-scale expansion (this error is essentially located close to $\partial K$);

• 3.

the term in $O(H)$ is a classical a priori error estimate contribution. It appears when approximating the homogenized solution (associated to $u^{\varepsilon}$) by means of the piecewise affine basis functions $\phi^{0}_{i}$.

2.2.3 Oversampling technique

The basis functions defined by (5) do not oscillate on the coarse mesh edges, which leads to a poor approximation of $u^{\varepsilon}$ close to these edges. To improve on that, a possible strategy is to impose affine Dirichlet boundary conditions on the boundary of a domain $S_{K}\supset K$ which is slightly larger than $K$ (and often homothetic to $K$, for practical reasons), and to only use the interior information, on $K$, to construct the basis functions. The motivation stems from the fact that the information inside $K$ depends only mildly on the (incorrect) boundary conditions imposed on $\partial S_{K}$. This so-called oversampling technique thus consists in first solving the local problems

[TABLE]

where $\boldsymbol{s}_{j}$ are the coordinates of the vertices of $S_{K}$, before defining the MsFEM basis functions in $K$ as $\displaystyle\phi_{i}^{\varepsilon}=\psi_{i|K}^{\varepsilon}$. In general, the obtained basis functions $\phi_{i}^{\varepsilon}$ are discontinuous at the boundaries of the macro elements of $\mathcal{T}_{H}$. This thus leads to a non-conforming MsFEM approximation: $V^{\varepsilon}_{H}=\text{Span }\{\phi_{i}^{\varepsilon}\}\not\subset V$. In this framework, the bilinear form in (4) is replaced by $\displaystyle B^{\varepsilon}_{H}(u,v)=\sum_{K\in\mathcal{T}_{H}}\int_{K}\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}u\cdot\boldsymbol{\nabla}v$, and the approximation $u^{\varepsilon}_{H}\in V^{\varepsilon}_{H}$ is defined as the solution to

[TABLE]

Remark 2

It is also possible to define $\phi_{i}^{\varepsilon}$ in $K$ as a linear combination of the functions $\psi^{\varepsilon}_{j|K}$. We refer e.g. to [47, Eq. (2.8)] for more details.

An a priori error estimate is given in [25] for the oversampling variant of MsFEM, under the same assumptions as for (6) (in particular, $\mathbb{A}^{\varepsilon}$ is again assumed to be periodic). It reads

[TABLE]

where the constant $C$ is independent of $\varepsilon$ and $H$ and where $\|\cdot\|_{1,H}$ is the broken $\mathcal{H}^{1}$-norm. Furthermore, a Petrov-Galerkin formulation is proposed in [35], with the aim to further decrease the cell resonance error (which occurs when $H\approx\varepsilon$). The formulation, in which the above non-conforming multiscale basis functions are used for trial functions while conforming linear Lagrange test functions are considered, reads

[TABLE]

2.2.4 Higher-order MsFEM

In the original version of MsFEM, the multiscale basis functions $\phi_{i}^{\varepsilon}$ are constructed from linear coarse-scale functions $\phi_{i}^{0}$. The construction of higher-order MsFEM basis functions was introduced in [5] using local harmonic coordinates. For each element $K\in\mathcal{T}_{H}$, these coordinates $w_{j}^{\varepsilon,K}$ ($j=1,\dots,d$) are defined as the solutions to

[TABLE]

We next define $w_{j}^{\varepsilon}\in\mathcal{H}^{1}(\Omega)$ by $w_{j|K}^{\varepsilon}=w_{j}^{\varepsilon,K}$, and set $\boldsymbol{w}^{\varepsilon}=(w_{1}^{\varepsilon},\dots,w_{d}^{\varepsilon})^{T}$. Introducing the set $\{\phi_{i}^{0}\}$ of the Lagrange basis functions of order $k$, the higher-order MsFEM basis $\{\phi_{i}^{\varepsilon}\}$ is constructed from the composition

[TABLE]

Choosing Lagrange functions $\{\phi_{i}^{0}\}$ of order $1$ leads to the classical MsFEM basis defined by (5). The approach is easy to implement since the computation of the oscillating functions $w_{j}^{\varepsilon}$ is independent of the order $k$ of the coarse scale basis. This approach is advantageously used in the following sections to both compute a posteriori error estimates and drive the adaptive strategy.

In the case when $\mathbb{A}^{\varepsilon}$ is periodic, an a priori error estimate is given in [5] when considering MsFEM basis functions of order $k$. It reads

[TABLE]

where the constant $C$ is independent of $\varepsilon$ and $H$.

Remark 3

As in the oversampling technique described in Section 2.2.3, the local problems (8) can be posed on a larger domain $S_{K}\supset K$ before using the restriction to $K$ of the solutions. In this case, the support of $\phi_{i}^{\varepsilon}$ may be different from that of $\phi_{i}^{0}$, and a non-conforming MsFEM approximation is obtained.

Remark 4

An alternative higher-order MsFEM approach was recently proposed in [32].

3 The CRE concept applied to MsFEM

In this section, we review the method we have introduced in [19] to obtain a rigorous and fully computable upper bound on the MsFEM error measured in the energy norm, i.e. an upper bound on ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u^{\varepsilon}-u^{\varepsilon}_{H}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$. It uses the concept of Constitutive Relation Error (CRE) based on duality [27, 55]. A detailed overview on CRE (in a single scale framework) can be read in [42].

3.1 Definition of an a posteriori error estimate

Introduce the space of equilibrated fluxes

[TABLE]

and define, for any pair $(\widehat{u},\widehat{\boldsymbol{p}})\in V\times\mathcal{S}$, the CRE functional $E_{CRE}$ as

[TABLE]

where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{F}}$ denotes the energy norm on flux fields. The equality

[TABLE]

referred to as the Prager-Synge equality, is a direct consequence of the property $\displaystyle\int_{\Omega}(\boldsymbol{q}^{\varepsilon}-\widehat{\boldsymbol{p}})\cdot\boldsymbol{\nabla}(u^{\varepsilon}-\widehat{u})=0$ (we recall that $\boldsymbol{q}^{\varepsilon}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}u^{\varepsilon}$). Another equality, derived from (9), reads

[TABLE]

With a given $\widehat{\boldsymbol{p}}\in\mathcal{S}$ at hand, a guaranteed and computable bound on the global MsFEM error ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u^{\varepsilon}-u^{\varepsilon}_{H}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ can be derived from (9), as we now explain. When using a conforming MsFEM, we have $V^{\varepsilon}_{H}\subset V$. We can hence choose $\widehat{u}=u^{\varepsilon}_{H}$ in (9), which directly yields

[TABLE]

When considering a nonconforming version of MsFEM, we may recover a field $\widehat{u}^{\varepsilon}_{H}\in V$ from $u^{\varepsilon}_{H}$ by a local averaging on the edges of the partition $\mathcal{T}_{H}$. Using the triangle inequality, we get

[TABLE]

where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{H}$ is the broken energy norm. The right-hand sides of (11) and (12) are fully computable.

Of course, the accuracy of the CRE estimate depends on the choice of $\widehat{\boldsymbol{p}}$. In view of (9), we see that the difference ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{q}^{\varepsilon}-\widehat{\boldsymbol{p}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathcal{F}}$ should be small in comparison with ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u^{\varepsilon}-\widehat{u}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$.

3.2 Computation of an equilibrated flux field

The technical point in the CRE concept is the construction of a relevant flux field $\widehat{\boldsymbol{p}}\in\mathcal{S}$. A convenient strategy, that we have introduced for all variants of MsFEM in [19] as an extension to that usually employed in the monoscale framework (see e.g. [40, 41, 60]), consists in post-processing the approximate MsFEM flux $\boldsymbol{q}^{\varepsilon}_{H}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}u^{\varepsilon}_{H}$ at hand in order to build a flux $\widehat{\boldsymbol{q}}^{\varepsilon}_{H}$ which belongs to $\mathcal{S}$. The strategy consists in two steps:

in Step 1, some tractions $\widehat{g}_{K}(\boldsymbol{x})$ are defined along the edges of each element $K\in\mathcal{T}_{H}$, with $\widehat{g}_{K}=g$ on $\Gamma_{N}$ (see Fig. 1). They satisfy an equilibrium property at the element level,

[TABLE]

as well as continuity between neighboring elements. Such tractions $\widehat{g}_{K}$ are constructed by solving local low-dimensional linear systems associated to each vertex in $\mathcal{T}_{H}$, and using the partition of unity property verified by the test functions $\phi_{i}^{\ast}$ used in the Galerkin MsFEM formulation, i.e. $\displaystyle\sum_{i}\phi_{i}^{\ast}=1$ ($\phi_{i}^{\ast}=\phi_{i}^{\varepsilon}$ in Galerkin formulations, while $\phi_{i}^{\ast}=\phi_{i}^{0}$ in Petrov-Galerkin formulations). We refer to [19] for details. The tractions are represented as affine combinations of the test functions $\phi_{i}^{\ast}$:

[TABLE]

where $K$ and $K^{\prime}$ are the elements such that $\Gamma=\overline{K}\cap\overline{K^{\prime}}$;

• 2.

in Step 2, the flux $\widehat{\boldsymbol{q}}^{\varepsilon}_{H}$ is constructed over each element $K$ as a solution to the local problem

[TABLE]

In practice, this problem is solved using a dual formulation. We hence search for some $w^{\varepsilon}\in\mathcal{H}^{1}(K)$ satisfying the Neumann problem

[TABLE]

and next set $\widehat{\boldsymbol{q}}^{\varepsilon}_{H|K}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}w^{\varepsilon}$. The Neumann problem (14) is well-posed (up to the addition of a constant) in view of (13). In practice, the problem (14) is numerically approximated by a Galerkin procedure on $\text{Span }\{\phi_{i}^{0}\circ\boldsymbol{w}^{\varepsilon}\}$ (see Section 2.2.4), where $\phi_{i}^{0}$ are Lagrange function of high-order (in practice, of order 4). Since $\widehat{g}_{K|\Gamma}$ is a linear combination of the functions $\phi_{i|\Gamma}^{\ast}$ and $\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}\phi_{i|\Gamma}^{\varepsilon}\cdot\boldsymbol{n}_{K}$, we can use the superposition principle and obtain $w^{\varepsilon}$ as a linear combination of pre-computed elementary solutions.

All technicalities in the construction of $\widehat{g}_{K}$ and $\widehat{\boldsymbol{q}}^{\varepsilon}_{H|K}$ are given in [19]. It is fruitful to notice that the overall strategy for constructing $\widehat{\boldsymbol{q}}^{\varepsilon}_{H}\in\mathcal{S}$ is particularly suited to the MsFEM framework. Indeed, it is based on fine-scale computations at the macro element level only. These can be performed in the offline stage of MsFEM while constructing the multiscale basis functions.

4 Goal-oriented error estimation

A posteriori error estimation on some outputs of the computation (which are typically needed for decision-making) is often very interesting for engineering applications. In many practical problems, one may indeed be interested in approximating some functional $Q(u^{\varepsilon})$, referred to as the quantity of interest, rather than the solution $u^{\varepsilon}$ itself. In this section, we address goal-oriented error estimation in the MsFEM framework, again using the CRE concept introduced in Section 3.

4.1 Adjoint problem and classical error indicator

We consider a linear and continuous functional $Q:V\to\mathbb{R}$, and assume that our quantity of interest is $Q(u^{\varepsilon})$. Typical examples include the average of the multiscale solution $u^{\varepsilon}$, or of the exact flux $\boldsymbol{q}^{\varepsilon}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}u^{\varepsilon}$, in some region of interest $\omega_{Q}\subset\Omega$. We then introduce the following adjoint problem [7, 29]:

[TABLE]

where $B^{\varepsilon\ast}(u,v)$ is the formal adjoint of the primal form $B^{\varepsilon}(u,v)$, satisfying $B^{\varepsilon\ast}(u,v)=B^{\varepsilon}(v,u)$. In the current symmetric case, the forms $B^{\varepsilon}$ and $B^{\varepsilon\ast}$ are identical. We also define the adjoint flux $\widetilde{\boldsymbol{q}}^{\varepsilon}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}\widetilde{u}^{\varepsilon}$ and the space of equilibrated fluxes (with regard to the adjoint problem):

[TABLE]

Remark 5

In practice, the functional $Q$ is often defined by

[TABLE]

where $\left(\widetilde{f}_{Q},\widetilde{\boldsymbol{q}}_{Q},\widetilde{g}_{Q}\right)$ is a set of extraction functions which can be mechanically interpreted in the adjoint problem as body, pre-flux and traction loadings, respectively.

Remark 6

To handle some particular cases (such as quantities of interest related to the normal flux on $\Gamma_{D}$, which correspond to a functional $Q$ defined on a strict subset of $V$), a more general formulation of the adjoint problem can be found in [29].

Remark 7

Pointwise quantities of interest, such as $u^{\varepsilon}(\boldsymbol{x}_{0})$ for some $\boldsymbol{x}_{0}\in\Omega$, cannot be directly considered in this framework since the solution $u^{\varepsilon}\in\mathcal{H}^{1}(\Omega)$ may not be continuous (for $d\geq 2$). To circumvent this issue, a regularization technique was introduced in [62] for (monoscale) FEM. We do not consider such pointwise quantities of interest here.

The primal-dual equivalence follows for $u^{\varepsilon}$ solution to the primal problem (2) and $\widetilde{u}^{\varepsilon}$ solution to the adjoint problem (15):

[TABLE]

Let $u^{\varepsilon}_{H}$ be the MsFEM approximation of $u^{\varepsilon}$, and let $\widehat{u}^{\varepsilon}_{H}$ be a function in $V$ close to $u^{\varepsilon}_{H}$ (if we work with a conforming MsFEM approach, then $\widehat{u}^{\varepsilon}_{H}=u^{\varepsilon}_{H}$). From the adjoint solution $\widetilde{u}^{\varepsilon}$, an exact representation of the MsFEM error on the quantity $Q(u^{\varepsilon})$ is available:

[TABLE]

with $\Delta Q=Q(\widehat{u}^{\varepsilon}_{H}-u^{\varepsilon}_{H})$, and where $R^{\varepsilon}(v)=F(v)-B^{\varepsilon}(\widehat{u}^{\varepsilon}_{H},v)$ is the residual functional associated with $\widehat{u}^{\varepsilon}_{H}\in V$. The representation (17) shows that the adjoint solution $\widetilde{u}^{\varepsilon}$, often termed generalized Green’s function, provides for the sensitivity of the error on $Q$ to the local sources of the MsFEM error in the whole domain $\Omega$.

A convenient error indicator may be derived from (17) using the dual weighted residual (DWR) method proposed in [6] for FEM and extended to multiscale computations in recent works (see for instance [2, 20]). The DWR method is based on hierarchical a posteriori error analysis. We compute an approximate adjoint solution $\widetilde{u}^{\varepsilon}_{+}$ from a hierarchically refined MsFEM space $V^{\varepsilon}_{+}\supset V^{\varepsilon}_{H}$ (higher-order basis functions, refined mesh) and we assume that the error between $\widetilde{u}^{\varepsilon}$ and $\widetilde{u}^{\varepsilon}_{+}$ is small. We then write

[TABLE]

where the right-hand side and $\Delta Q$ are fully computable.

Remark 8

Due to the Galerkin orthogonality, choosing the same MsFEM space $V^{\varepsilon}_{H}$ to compute $u^{\varepsilon}_{H}$ and $\widetilde{u}^{\varepsilon}_{+}$ would lead to a meaningless indicator (in that case, we would have $R^{\varepsilon}(\widetilde{u}^{\varepsilon}_{+})=0$). Consequently, a convenient approximation of the adjoint solution should involve a subspace in $(V^{\varepsilon}_{H})^{\perp}$.

The above argument can be somewhat quantified. Let $u^{\varepsilon}_{+}\in V^{\varepsilon}_{+}$ be the approximation (in the refined space) of the primal solution. We assume that the saturation assumption $|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{+})|\leq\beta|Q(u^{\varepsilon})-Q(\widehat{u}^{\varepsilon}_{H})|$ is satisfied for some $0\leq\beta<1$. Using the fact that $u^{\varepsilon}_{+}-\widehat{u}^{\varepsilon}_{H}\in V^{\varepsilon}_{+}$ and next the fact that $\widetilde{u}^{\varepsilon}_{+}\in V^{\varepsilon}_{+}$, we obtain

[TABLE]

Using the triangle inequality, we get

[TABLE]

which quantifies the approximation (18). The indicator defined in (18) may be used to drive adaptive algorithms, but it does not provide for guaranteed error bounds on $Q$. Alternative techniques are proposed in the next sections.

4.2 Guaranteed error bounds using CRE

In the FEM context, classical approaches to bound the error on a quantity of interest $Q$ consist in solving the primal and adjoint problems with the same approximation space, before using a posteriori bounds available for the error in the energy norm [3]. These approaches can be easily extended to MsFEM with general, conforming or nonconforming, formulations.

Denoting by $\widetilde{u}^{\varepsilon}_{H}$ the approximate adjoint solution obtained using the space $V^{\varepsilon}_{H}$, and potentially post-processing it to define $\widehat{\widetilde{u}}^{\varepsilon}_{H}\in V$ in case $V^{\varepsilon}_{H}\not\subset V$ (see Section 3.1), the equality (17) yields that

[TABLE]

The Cauchy-Schwarz inequality then yields

[TABLE]

Remark 9

In the conforming formulation of MsFEM (i.e. when $V^{\varepsilon}_{H}\subset V$), we can take $\widehat{u}^{\varepsilon}_{H}=u^{\varepsilon}_{H}$ and $\widehat{\widetilde{u}}^{\varepsilon}_{H}=\widetilde{u}^{\varepsilon}_{H}$ so that the Galerkin orthogonality $\displaystyle R^{\varepsilon}\left(\widehat{\widetilde{u}}^{\varepsilon}_{H}\right)=0$ holds. The result (20) can be recast as

[TABLE]

This is a classical result in the monoscale conforming FEM context.

Consequently, guaranteed error bounds on $Q$ can be obtained from any guaranteed a posteriori global error estimate on ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u^{\varepsilon}-\widehat{u}^{\varepsilon}_{H}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{u}^{\varepsilon}-\widehat{\widetilde{u}}^{\varepsilon}_{H}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$. Using the CRE concept introduced in Section 3, (20) yields the a posteriori estimate

[TABLE]

where $\widehat{\boldsymbol{q}}^{\varepsilon}_{H}\in\mathcal{S}$, resp. $\widehat{\widetilde{\boldsymbol{q}}}^{\varepsilon}_{H}\in\widetilde{\mathcal{S}}$, is constructed as a post-processing of $\boldsymbol{q}^{\varepsilon}_{H}$, resp. $\widetilde{\boldsymbol{q}}^{\varepsilon}_{H}$, as detailed in Section 3.2. However, the estimate (21), based on the Cauchy-Schwarz inequality over the whole domain $\Omega$, does not exploit the locality of $Q$. It thus leads to crude upper bounds when MsFEM errors for primal and adjoint problems have disjoint supports. A better estimate would consist in using the Cauchy-Schwarz inequality at the element level, under the form

[TABLE]

However, computable upper bounds on the error in the energy norm restricted to each element $K$ are hardly available, even though recent works such as [43] have proposed first ideas in that direction.

An alternative and more accurate estimate than (21) can be constructed by taking advantage of the property (10) of the CRE concept. Indeed, we infer from (19) that

[TABLE]

with $\widehat{\boldsymbol{q}}^{\varepsilon\ast}_{H}=(\widehat{\boldsymbol{q}}^{\varepsilon}_{H}+\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}\widehat{u}^{\varepsilon}_{H})/2$ and $C^{\varepsilon}_{H}$ a computable term defined by

[TABLE]

We next get from (22) that

[TABLE]

with

[TABLE]

Using the Cauchy-Schwarz inequality and next (10), we obtain that

[TABLE]

Note that, thanks to the factor $1/2$ on the right-hand side, the estimate (23) yields a more accurate estimation of $Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})$ than (21).

The fully computable quantity $Q(u^{\varepsilon}_{H})+\Delta Q+\overline{C}^{\varepsilon}_{H}$ can be interpreted as a corrected approximation of $Q(u^{\varepsilon})$. Furthermore, the estimate (23) can also be written as

[TABLE]

hence

[TABLE]

In contrast to (21), the quantity $\eta^{Q}$ partially takes into account error cancellations.

It is fruitful to notice that the bounds (21) and (23), constructed from the CRE concept, hold for independent (MsFEM) discretizations of the primal and adjoint problems as they are merely based on equilibrium properties. In the following section, we show how a convenient local enrichment of the adjoint approximate solution can thus be advantageously used to get sharp error bounds from (23).

4.3 Non-intrusive local enrichment of the adjoint solution

We extend to the MsFEM context an enrichment technique initially introduced in [30, 17, 39] for goal-oriented error estimation in FEM models. The general idea is to improve the accuracy of the approximate adjoint solution $\widetilde{u}^{\varepsilon}_{H}$ (or $\widehat{\widetilde{u}}^{\varepsilon}_{H}$) by adding a priori information on the problem, i.e. given functions whose purpose is to reproduce localized high-gradient components of $\widetilde{u}^{\varepsilon}$ (which are difficult to numerically capture) in the vicinity of the adjoint loading region. These functions, referred to as handbook functions and denoted $\widetilde{u}^{\varepsilon}_{hand}$ (with associated flux $\widetilde{\boldsymbol{q}}^{\varepsilon}_{hand}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}\widetilde{u}^{\varepsilon}_{hand}$), correspond to generalized Green’s functions. In analogy to (15), they satisfy

[TABLE]

where $Q_{hand}$ is an auxiliary quantity of interest (the eventual quantity of interest $Q$ in (15) may be equal to $Q_{hand}$, or a linear combination of various $Q_{hand}$, or various translations of $Q_{hand}$, …). The domain $\Omega_{hand}\subset\Omega$ is chosen to be sufficiently large, so that the solution to the above problem with $\mathcal{H}^{1}_{0}(\Omega_{hand})$ replaced by $V$ would indeed be close to zero on $\partial\Omega_{hand}$. In practice, $\widetilde{u}^{\varepsilon}_{hand}$ is numerically pre-computed, using an accurate fine scale discretization. Examples of such enrichment functions, for a 2D multiscale medium described by $\mathbb{A}^{\varepsilon}=A^{\varepsilon}\,\mathbb{I}_{2}$ with $A^{\varepsilon}(x_{1},x_{2})=3+\cos(2\pi x_{1}/\varepsilon)+\cos(2\pi x_{2}/\varepsilon)$ and $\varepsilon=1/20$ (see Fig. 2), are given below.

A first example corresponds to a linear functional $Q$ given by (16) with $\widetilde{\boldsymbol{q}}_{Q}=0$, $\widetilde{g}_{Q}=0$ and a localized body force $\widetilde{f}_{Q}=0$ except in the square domain $\omega_{Q}$, centered at the center of $\Omega$ and of size $4\varepsilon\times 4\varepsilon$, where $\widetilde{f}_{Q}=1$. Such a quantity of interest corresponds to the average of $u^{\varepsilon}$ in $\omega_{Q}$. The corresponding handbook function is shown in Fig. 3.

Remark 10

Here and in all what follows, the quantity of interest is given by (16) for some $\widetilde{f}_{Q}$ and $\widetilde{\boldsymbol{q}}_{Q}$, with $\widetilde{g}_{Q}=0$. We take for $\Omega_{hand}$ a domain 10 times as large (in each direction) as the support of $\widetilde{f}_{Q}$ and $\widetilde{\boldsymbol{q}}_{Q}$. For the example discussed above, we thus choose for $\Omega_{hand}$ a square domain centered at the center of $\Omega$ and of size $40\varepsilon\times 40\varepsilon$.

A second example corresponds to a linear functional $Q$ given by (16) with $\widetilde{f}_{Q}=0$, $\widetilde{g}_{Q}=0$ and a localized pre-flux $\widetilde{\boldsymbol{q}}_{Q}=(1,0)^{T}$ in the same domain $\omega_{Q}$ as above, and $\widetilde{\boldsymbol{q}}_{Q}=0$ elsewhere in $\Omega$. Such a quantity of interest corresponds to the average of $\boldsymbol{e}_{1}\cdot\nabla u^{\varepsilon}$ in $\omega_{Q}$. The corresponding handbook function is shown in Fig. 4.

Handbook functions are then locally and seamlessly inserted in the overall procedure, as we explain now. Let $\Omega_{Q}$ be the area of interest over which the quantity $Q$ is defined, i.e. the region which supports the extraction functions $\widetilde{f}_{Q}$ and $\widetilde{\boldsymbol{q}}_{Q}$ (recall (16) and Remark 10). We next introduce some domain $\Omega_{1}^{PUM}\subset\Omega_{hand}$, which is the union of some coarse elements of $\mathcal{T}_{H}$, and which contains $\Omega_{Q}$ (see Fig. 5). Its boundary is denoted $\partial\Omega_{1}^{PUM}$. Let $\mathcal{I}^{PUM}$ be the subset of vertices in $\mathcal{T}_{H}$ which are located in $\overline{\Omega_{1}^{PUM}}$, so that

[TABLE]

We next introduce

[TABLE]

and see that $\displaystyle\Omega_{1}^{PUM}\cup\Omega_{2}^{PUM}=\text{Supp}\left(\sum_{i\in\mathcal{I}^{PUM}}\phi^{0}_{i}\right)\subset\Omega$.

Remark 11

In the numerical examples of Section 6.1, $\Omega_{Q}$ is a square of size $4\varepsilon\times 4\varepsilon$, and we take for $\Omega_{1}^{PUM}$ the coarse element containing $\Omega_{Q}$. In the numerical example of Section 6.2, $\Omega_{Q}$ is a square of size $0.1\times 0.1$, and the initial coarse meshsize is $H=0.1$. We then take for $\Omega_{1}^{PUM}$ the union of the coarse element equal to $\Omega_{Q}$ and of the 8 coarse elements surrounding it.

The enrichment approach consists in approximating the adjoint solution $\widetilde{u}^{\varepsilon}$ by

[TABLE]

where $\widetilde{u}^{\varepsilon}_{res,H}$ is an approximate residual solution which is computed using the MsFEM approximation space $V^{\varepsilon}_{H}$. We emphasize that

no additional degree of freedom is associated with the enrichment part, in contrast to other similar enrichment methods such as XFEM or GFEM (where the approximation would read $\displaystyle\sum_{i\in\mathcal{I}^{PUM}}\alpha_{i}\,\phi^{0}_{i}\,\widetilde{u}^{\varepsilon}_{hand}+\widetilde{u}^{\varepsilon}_{res,H}$ for some additional unknown coefficients $\alpha_{i}$ and with $\widetilde{u}^{\varepsilon}_{res,H}\in V^{\varepsilon}_{H}$). This is due to the fact that $\widetilde{u}^{\varepsilon}_{hand}$ aims at representing the local high-gradients due to the singular adjoint loading, which can be fully assessed a priori;

• 2.

the strategy uses the same MsFEM discretization space $V^{\varepsilon}_{H}$ to approximate both the primal and the adjoint problems, even though a higher level of approximation is obtained for the adjoint solution due to the enrichment term. Consequently, all the numerical technicalities (mesh connectivities, stiffness matrix, …) implemented to solve the primal problem can be directly reused to solve the adjoint problem. Only the loadings need to be changed. For this reason, the strategy is called non-intrusive. It fully makes sense in the MsFEM framework where the discretization space is obtained from costly computations performed only once in the offline stage.

Remark 12

The above enrichment methodology can also be used to address goal-oriented error estimation on pointwise quantities of interest (when these are defined), as shown in [17, 39]. In this case, handbook functions correspond to classical Green’s functions. They can be analytically computed for homogeneous and isotropic media.

Remark 13

We have chosen the notations $\Omega_{1}^{PUM}$, $\Omega_{2}^{PUM}$ and $\mathcal{I}^{PUM}$ in analogy to those introduced in the XFEM and the GFEM approaches. However, here, no partition of unity property is used.

In view of (15) and (26), and assuming that $Q_{hand}\equiv Q$ in (25), the problem to be solved numerically consists in finding $\widetilde{u}^{\varepsilon}_{res,H}\in V^{\varepsilon}_{H}$ such that

[TABLE]

with $\boldsymbol{n}_{12}$ is the outgoing normal vector to $\Omega_{1}^{PUM}$. In the last line, we have used the fact (assuming e.g. that $\widetilde{g}_{Q}=0$ and $\widetilde{\boldsymbol{q}}_{Q}=0$ in (16)) that $-\boldsymbol{\nabla}\cdot(\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}\widetilde{u}^{\varepsilon}_{hand})=\widetilde{f}_{Q}$ in $\Omega_{hand}\supset\Omega_{1}^{PUM}$ (see (25)), and that $\displaystyle\int_{\Omega_{1}^{PUM}}\widetilde{f}_{Q}\,v=\int_{\Omega}\widetilde{f}_{Q}\,v$ because $\widetilde{f}_{Q}$ is supported in $\Omega_{Q}\subset\Omega_{1}^{PUM}$.

Post-processing the numerical flux $\widetilde{\boldsymbol{q}}^{\varepsilon}_{res,H}=\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}\widetilde{u}^{\varepsilon}_{res,H}$ as explained in Section 3.2, we compute a flux field $\widehat{\widetilde{\boldsymbol{q}}}^{\varepsilon}_{res,H}\in\widetilde{\mathcal{S}}_{res}$, with

[TABLE]

We next define

[TABLE]

We also introduce

[TABLE]

where $\widehat{\widetilde{u}}^{\varepsilon}_{res,H}\in V$ is equal to $\widetilde{u}^{\varepsilon}_{res,H}$ when $V_{H}^{\varepsilon}\subset V$, and is obtained after a post-processing of $\widetilde{u}^{\varepsilon}_{res,H}$ otherwise (see Section 3.1).

Using the estimate (23), we eventually obtain

[TABLE]

with

[TABLE]

The fully computable quantity $Q(u^{\varepsilon}_{H})+\Delta Q+\overline{C}^{\varepsilon}_{H+}$ is a corrected approximation of $Q(u^{\varepsilon})$. As in (24), the quantity

[TABLE]

is a computable upper bound of $|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})|$. This bound is thought to be sharp as the enrichment of the adjoint solution is supposed to make the term $\displaystyle E_{CRE}\left(\widehat{\widetilde{u}}^{\varepsilon}_{res,H},\widehat{\widetilde{\boldsymbol{q}}}^{\varepsilon}_{res,H}\right)$ small. When $\eta^{Q}$ is higher than a given error tolerance $\eta_{tol}^{Q}$, an adaptive strategy can be introduced. This is the subject of the next section.

5 Goal-oriented adaptive strategy

The objective here is to design an adaptive algorithm that optimally modifies the parameters of the MsFEM discretization in order to reduce the error on $Q$ and reach a given error tolerance $\eta_{tol}^{Q}$ while keeping the computational cost as small as possible. This is performed by introducing local error indicators associated to each error source and defined from (hierarchical) auxiliary reference problems and CRE properties.

5.1 Recast of MsFEM

We first introduce a hierarchy of MsFEM-type problems, which is identical to that introduced in [19] (we also refer to [31], where a similar hierarchy is introduced). We briefly present them below, and refer to [19, Section 5.1] for more details.

We first introduce a projection operator $I_{H}:V\to V^{0}_{H}\subset V$ chosen such that its image is $V^{0}_{H}$, where $V^{0}_{H}$ is spanned by the Lagrange basis functions associated with each node of $\mathcal{T}_{H}$ (see Section 2.1). A possible choice for $I_{H}$ is the Clément interpolant. We next define the fine-scale space

[TABLE]

and observe that $V=V^{0}_{H}\oplus V^{f}$.

For any function $\varphi_{H}\in V^{0}_{H}$, we define the corrector $\Theta^{\varepsilon}(\varphi_{H})\in V^{f}$ as the solution to the problem

[TABLE]

The associated operator $\Theta^{\varepsilon}:V^{0}_{H}\to V^{f}$ is an affine operator:

[TABLE]

We then define a reconstruction operator $\Lambda^{\varepsilon}:V^{0}_{H}\to V$ as $\Lambda^{\varepsilon}=I+\Theta^{\varepsilon}$, as well as the multiscale functional space \displaystyle\overline{V}^{\varepsilon}_{H}=\{\Lambda^{\varepsilon}(\varphi_{H}),\ \varphi_{H}\in V^{0}_{H}\}=\Theta^{\varepsilon}_{0}+\text{Span}\big{\{}\phi^{0}_{i}+\Theta_{\ell}^{\varepsilon}(\phi^{0}_{i})\big{\}}, where $\{\phi^{0}_{i}\}$ is the set of Lagrange basis functions of $V^{0}_{H}$.

Consider the following MsFEM-like problem:

[TABLE]

which is equivalent to the problem

[TABLE]

It is easy to show that (29) and (30) are well-posed, and that their solution $\overline{u}^{\varepsilon}_{H}$ satisfies the variational formulation (2) of the reference problem. We hence have $\overline{u}^{\varepsilon}_{H}=u^{\varepsilon}$. Note in particular that the fine scale problems (28) are solved exactly, over the whole domain $\Omega$, and considering the reference loading $F$ on the right-hand side.

Starting from the exact formulations (29) or (30), we now construct a fully discrete MsFEM formulation, on the MsFEM space $V^{\varepsilon}_{H}$, by making the following simplifying changes (denoted S#):

S1: replace the right-hand side term with [math] in the fine-scale problem (28), which makes the operator $\Theta^{\varepsilon}$ linear (i.e. $\Theta^{\varepsilon}_{0}=0$) and independent of the loadings. This leads to the low-dimensional vector space \text{Span}\{\Lambda^{\varepsilon}(\phi^{0}_{i})\}=\text{Span}\big{\{}\phi^{0}_{i}+\Theta_{\ell}^{\varepsilon}(\phi^{0}_{i})\big{\}};

• 2.

S2: restrict to $\Omega_{i}$ the correction $\Theta^{\varepsilon}(\phi^{0}_{i})$, where $\Omega_{i}$ is the support of $\phi^{0}_{i}$. This leads to the vector space

[TABLE]

• 3.

S3: localize the fine-scale problem (28) by replacing $\Omega$ by smaller domains. In practice, these are subdomains $S_{K}$ that include the element $K$, on the boundary $\partial S_{K}$ of which we prescribe homogeneous Dirichlet boundary conditions;

• 4.

S4: replace $V^{f}$ by a discrete $P_{1}$-Lagrange FE space $V^{f}_{h}$ associated to a fine grid $\mathcal{T}_{h}$ obtained from a regular refinement of $\mathcal{T}_{H}$. The fine mesh size $h$ should be chosen so that oscillations of the data can be accurately captured.

Applying the simplifications S1-S4 on (29) or (30), we obtain a fully discrete MsFEM formulation. For each $K\in\mathcal{T}_{H}$ and each basis function $\phi^{0}_{i}\in V^{0}_{H}$, the local correction $\Theta^{\varepsilon}_{h}(\phi^{0}_{i|K})\in V^{f}_{h}(S_{K})$ is defined as the solution to

[TABLE]

where $\pi(\phi^{0}_{i|K})$ is the affine function on $S_{K}$ which is equal to $\phi^{0}_{i|K}$ on $K$, and $V^{f}_{h}(S_{K})$ is a discrete subspace of $V^{f}(S_{K})=\mathcal{H}^{1}_{0}(S_{K})$. The local correctors $\Theta^{\varepsilon}_{h}(\phi^{0}_{i|K})$ are then connected together in the global corrector $\displaystyle\Theta^{\varepsilon}_{h}(\phi^{0}_{i}):=\sum_{K\in\mathcal{T}_{H}}\Theta^{\varepsilon}_{h}(\phi^{0}_{i|K})_{|K}$. The MsFEM functional space eventually reads $V^{\varepsilon}_{H}=\text{Span}\{\phi^{0}_{i}+\Theta^{\varepsilon}_{h}(\phi^{0}_{i})\}$.

Remark 14

Applying the simplifications S1-S4 on (29) with the choice $S_{K}=K$ yields the MsFEM approach without oversampling given by (4)–(5).

5.2 Splitting of error sources and error indicators

The MsFEM formulation given in Section 5.1 shows that 3 parameters can be tuned in an adaptive procedure to improve the accuracy of $Q(u^{\varepsilon}_{H})$:

• $\bullet$

the size $H_{K}$ of the elements of the coarse mesh $\mathcal{T}_{H}$. Simplification S1 amounts to ignoring the influence of $f$ and $g$ in micro-scale equations. This assumption is valid provided $H_{K}$ is comparable with the characteristic length of spatial variations of the external loading and macroscale solution. As detailed in [19, Remark 16], the smaller $H_{K}$ is, the better Simplification S1 is justified. The validity of Simplification S2 is also related to $H_{K}$. If the set $\{H_{K}\}$ is not chosen correctly, the coarse-scale discretization error can dominate the overall approximation error;

• $\bullet$

the size $d_{K}$ (minimum patch radius) of the computational domains $S_{K}$ used to solve the micro-scale equations obtained after Simplification S3. If the set $\{d_{K}\}$ is not chosen correctly, the artificial boundary conditions set on $\partial S_{K}$ may strongly dominate the overall approximation error;

• $\bullet$

the local size $h_{K}$ of the fine mesh used to solve the micro-scale equations (in each $S_{K}$) after applying Simplification S4. If the set $\{h_{K}\}$ is not chosen correctly, a fine-scale discretization error can dominate the overall approximation error.

We wish to define a robust adaptive algorithm that is able: (i) to detect regions where the MsFEM approximation is not sufficiently accurate for the quantity of interest $Q(u^{\varepsilon})$; (ii) to adapt locally and optimally the relevant MsFEM parameters among the three ones mentioned above.

The item (i) can be directly addressed with a spatial decomposition, over the macro mesh $\mathcal{T}_{H}$, of the a posteriori error estimate $\eta_{Q}$ defined by (27):

[TABLE]

with

[TABLE]

where $\theta_{\rm max}$ is the maximizer in (27). The contribution $\eta^{Q}_{K}$, which is fully computable, is considered as a local error indicator in the coarse element $K$.

The above item (ii) requires the definition of specific error indicators associated to each individual error source. In order to set up such error indicators, we follow our approach described in details in [19, Section 5.2] and introduce two intermediate reference problems:

the first one, denoted PR1 in the following, is obtained from the initial reference problem (2) (recast as (29) or (30) and denoted PR0) by applying Simplifications S1 and S2. The functional space associated with PR1 is the space $\overline{V}^{\varepsilon}_{H,1}$ defined by (31), and the solution to PR1 is denoted $\overline{u}_{H,1}^{\varepsilon}$. Note that, in the problem PR1, fine-scale computations are still defined over the whole domain $\Omega$;

• 2.

the second one, denoted PR2 in the following, is obtained from the initial reference problem PR0 by applying Simplifications S1, S2 and S3. The functional space associated with PR2 is denoted $\overline{V}^{\varepsilon}_{H,2}$, and the solution to PR2 is denoted $\overline{u}_{H,2}^{\varepsilon}$. In the problem PR2, the fine-scale computations are defined over the subdomains $S_{K}$ and are performed exactly (no discretization of $V^{f}(S_{K})=\mathcal{H}^{1}_{0}(S_{K})$ is yet introduced).

Note that $\overline{u}_{H,1}^{\varepsilon}$ and $\overline{u}_{H,2}^{\varepsilon}$ may be nonconforming. We next write

[TABLE]

where

• $\bullet$

$e^{Q}_{macro}=|Q(u^{\varepsilon}-\overline{u}_{H,1}^{\varepsilon})|$ is the part of the error on $Q$ due to the coarse-scale discretization;

• $\bullet$

$e^{Q}_{over}=|Q(\overline{u}_{H,1}^{\varepsilon}-\overline{u}_{H,2}^{\varepsilon})|$ is the part of the error on $Q$ due to the localization on $S_{K}$ of the problems (28);

• $\bullet$

$e^{Q}_{micro}=|Q(\overline{u}_{H,2}^{\varepsilon}-u^{\varepsilon}_{H})|$ is the part of the error on $Q$ due to the fine-scale discretization.

The strategy to assess $e^{Q}_{macro}$, $e^{Q}_{over}$ and $e^{Q}_{micro}$ uses error indicators constructed from the CRE concept, with admissibility in the weaker sense of the associated reference problems (both for the primal and the adjoint problems). This is illustrated in Fig. 6. We have used the same ideas in [19, Section 5.2] to define specific error indicators (on the whole field $u^{\varepsilon}$) associated to each individual error source. Error indicators $\eta^{Q}_{macro}$, $\eta^{Q}_{over}$ and $\eta^{Q}_{micro}$ associated to each error source can then be computed. We emphasize that the expensive computations (involving the resolution of highly oscillatory problems) can all be performed in the offline stage of our approach. These indicators are used to drive a goal-oriented greedy adaptive algorithm detailed in the next section.

5.3 Adaptive algorithm

From the error estimate $\eta^{Q}$ defined by (27) and the error indicators ($\eta^{Q}_{macro}$, $\eta^{Q}_{over}$, $\eta^{Q}_{micro}$) defined in Section 5.2, we design a goal-oriented greedy adaptive algorithm that yields an adaptive MsFEM approach. The overall algorithm is presented in Fig. 7 and is similar to the algorithm designed in [19, Section 5.3] for the control of the MsFEM error measured in the energy norm.

The algorithm is based on a splitting (32) of the goal-oriented error estimator $\eta^{Q}$ into local contributions $\{\eta^{Q}_{macro,K},\eta^{Q}_{over,K},\eta^{Q}_{micro,K}\}$ in each element $K$ of the mesh $\mathcal{T}_{H}$. These are used to mark coarse elements of the mesh $\mathcal{T}_{H}$ for further enrichment, and thus to drive the adaptive process by selecting optimal MsFEM parameters $(H_{K},d_{K},h_{K})$ in each coarse region $K$ to compute both primal and adjoint approximate solutions.

Each iteration in the algorithm is made of three groups of computations: (i) microscale (possibly costly) computations that can be performed in the offline phase of the iteration; (ii) macroscale computations that can be performed in this offline phase; (iii) macroscale computations that need to be performed in the online phase of the iteration, but which are inexpensive.

The adaptive procedure is associated with the following considerations:

• $\bullet$

the algorithm is initialized with the coarsest MsFEM configuration that can be employed (regular coarse mesh, no oversampling, fine-scale problems solved with a single macro element). In this configuration, the diffusion tensor that locally appears to compute the MsFEM shape functions is the average of $\mathbb{A}^{\varepsilon}$ on $K$, as discussed in [19, Remark 20];

• $\bullet$

when modifying the parameters $h_{K}$ from their initial value $h_{K}=H_{K}$ (in Step 5.5), two values are considered in the adaptive process: $h_{K}=\varepsilon/5$ or $h_{K}=\varepsilon/20$ (finest mesh size at the microscale);

• $\bullet$

when modifying the parameters $d_{K}$ from their initial value $d_{K}=H_{K}$ i.e. $S_{K}=K$ (in Step 5.6), the oversampling size is determined by adding layers of progressive thickness $\varepsilon$, $2\varepsilon$, $3\varepsilon$, …, around the element $K$. Such an adaptation is not performed when $h_{K}=H_{K}$;

• $\bullet$

refining the mesh $\mathcal{T}_{H}$, i.e. modifying the parameters $H_{K}$ (in Step 5.7), is performed using a quadtree (or octree in 3D) method with nested elements. This requires to handle hanging nodes;

• $\bullet$

the adapted parameters $h_{K}$, $d_{K}$ and $H_{K}$ are defined from classical adaptive strategies (see [42]), based on the convergence rates given by the a priori estimates (6) or (7), and with an error target $\eta^{Q}_{macro/over/micro}=\eta^{Q}_{tol}/3$;

• $\bullet$

the technique proposed in [5] and based on composition rules (see Section 2.2.4) is beneficially used throughout the adaptive process. It enables independent computations without coming back to the fine scale offline computations (no additional costly computations). The technique is used in its $p$-refinement version in Steps 4.3 and 4.4, while it is used in its $H$-refinement version in Step 4.2 (when the macro mesh has to be refined). In short, this technique allows us to refine the coarse mesh without the need to solve new fine-scale problems of the type (8).

The adaptive algorithm allows us to adjust automatically (and hopefully optimally) the calculation parameters during the simulation. We expect that it provides MsFEM parameters which ensure a correct balance between computational cost and accuracy on the output of interest. In particular, it should lead to fine-scale solutions only in the parts of $\Omega$ which have an impact on the value of the quantity of interest. As will be shown below, and as expected, it also leads (for the same tolerance on the error) to less expensive computations than the adaptive algorithm introduced in [19] to drive the MsFEM computations according to the error measured in the energy norm.

6 Numerical results

We consider two test-cases. The first one corresponds to a structure with a local defect, for which we consider two possible loadings (see Section 6.1). The second one is a flow problem in fractured porous media (see Section 6.2), which has already been considered in [20]. For all these test-cases, the conclusions of our numerical investigations are the same:

It is worth working with a goal-oriented adaptive strategy rather than an adaptive strategy based on the global error: for the same error tolerance, the computational cost is much lower, in the sense that the obtained mesh (satisfying the prescribed error tolerance) is much coarser.

• 2.

As the mesh is iteratively adapted, the error on the quantity of interest monotically decreases until reaching the prescribed tolerance.

• 3.

The error estimate $\eta^{Q}$ is always an upper bound on the exact error, in view of (27). We see below that it only slightly overestimates the error (roughly by 10% on the considered numerical examples), showing thus its accuracy.

6.1 Structure with a local defect

To start with, we consider a highly oscillatory problem with defects, as introduced and analyzed in details in [10, 11, 12] (see also the recent works [13, 14]). We assume that the oscillating coefficient is a local perturbation (which models the defect) of a function with simple geometric properties, for instance a periodic function. We thus consider problem (1) with a (scalar, to simplify) coefficient $A^{\varepsilon}(\boldsymbol{x})=a(\boldsymbol{x}/\varepsilon)$ with $a=a_{per}+b$ where $b\in L^{r}(\mathbb{R}^{d})$ is the local perturbation. The homogenization of such a problem is addressed in [10] for the Hilbertian case $r=2$ and in [11, 12] for the cases $1<r<+\infty$ (see also [13, 14]). In general, the defect does not modify the macroscopic behavior (i.e. the homogenized solution) but it affects the solution locally, and at the small scale.

Inspired by these works, we consider here the following example. We set $\Omega=(-1,1)^{2}$, choose $\varepsilon=1/N$ for some integer $N$, and consider problem (1) with $\mathbb{A}^{\varepsilon}(\boldsymbol{x})=A^{\varepsilon}(\boldsymbol{x})\mathbb{I}$, where $A^{\varepsilon}(\boldsymbol{x})=a_{per}(\boldsymbol{x}/\varepsilon)+b(\boldsymbol{x}/\varepsilon)$ with

[TABLE]

In the following, we choose $\varepsilon=1/20$. The evolution of $A^{\varepsilon}$ is shown on Fig. 8.

6.1.1 Sinusoidal loading

We first consider the loading $f(x_{1},x_{2})=\sin(\pi x_{1})\cos(\pi x_{2})$. We take an initial mesh $\mathcal{T}_{H}$ made of $9\times 9$ macro elements, with $h_{K}=\varepsilon/3$ for any $K$ and no oversampling. The approximate MsFEM solution $u^{\varepsilon}_{H}$ is represented below and compared to the exact solution $u^{\varepsilon}$, which is obtained in practice through an expensive computation over a $500\times 500$ fine mesh. The solutions $u^{\varepsilon}$ and $u^{\varepsilon}_{H}$ are shown on Fig. 9, the gradients $\boldsymbol{\nabla}u^{\varepsilon}$ and $\boldsymbol{\nabla}u^{\varepsilon}_{H}$ are shown on Fig. 10, and the fluxes $\boldsymbol{q}^{\varepsilon}$ and $\boldsymbol{q}^{\varepsilon}_{H}$ are shown on Fig. 11. We also show the distribution of the energy of the exact solution and of the global error on Fig. 12.

We also show on Fig. 13 the distribution of the error on $u^{\varepsilon}$, $\boldsymbol{\nabla}u^{\varepsilon}\cdot\boldsymbol{e}_{1}$ and $\boldsymbol{q}^{\varepsilon}\cdot\boldsymbol{e}_{1}$.

We choose the quantity of interest $\displaystyle Q(u^{\varepsilon})=\frac{1}{|\omega|}\int_{\omega}\boldsymbol{q}^{\varepsilon}\cdot\boldsymbol{e}_{1}=\frac{1}{|\omega|}\int_{\omega}\boldsymbol{e}_{1}\cdot\mathbb{A}^{\varepsilon}\boldsymbol{\nabla}u^{\varepsilon}$ with $|\omega|=4\varepsilon\times 4\varepsilon$. We consider two cases:

Case 1: the subregion $\omega$ is centered on the defect ($Q=Q_{1}$);

• 2.

Case 2: the subregion $\omega$ is far from the defect ($Q=Q_{2}$).

Note that $u^{\varepsilon}$ is close to zero in the region around the defect (see Fig. 9), so considering as quantity of interest the average of $u^{\varepsilon}$ would not be interesting. In contrast, $\boldsymbol{q}^{\varepsilon}\cdot\boldsymbol{e}_{1}$ is not close to zero in that region (see Fig. 11). The adjoint solutions $\widetilde{u}^{\varepsilon}$ and fluxes $\widetilde{\boldsymbol{q}}^{\varepsilon}$ corresponding to the quantities of interest $Q_{1}$ and $Q_{2}$ are shown on Figs. 14 and 15. We observe, as expected, that the adjoint solutions are essentially supported in the domain of interest $\omega$.

The error threshold is fixed to $1\%$. We show on Fig. 16 (resp. Fig. 17) the final MsFEM discretization obtained after using the adaptive algorithm presented in Section 5.3, in order to respect the error threshold for the quantity of interest $Q_{1}$ (resp. $Q_{2}$). This discretization is to be compared with that obtained when controlling the global error (in energy norm) using the tools presented in [19], which is shown on Fig. 18. Obviously, the mesh on Figs. 16 and 17 is much coarser than the mesh on Fig. 18. This shows, as expected, that it is worth working with a goal-oriented adaptive strategy, rather than an adaptive strategy based on the global error: for the same error, the computational cost is much lower. Comparing Figs. 16 and 17, we also observe that the mesh is only refined close to the region of interest $\omega$. In the second example, when this region is far from the defect, it is not needed to adapt the mesh close to the defect.

As regards the accuracy of the error estimate (with respect to the actual error), we observe that $\eta^{Q}/|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})|=1.11$ (resp. 1.13) at the first iteration of the adaptive algorithm, and $\eta^{Q}/|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})|=1.07$ (resp. 1.06) at the final iteration of the adaptive algorithm, for the quantities of interest $Q_{1}$ and $Q_{2}$, respectively. As pointed out above, the error estimate thus only very slightly overestimates the actual error.

We eventually show on Fig. 19 the convergence of the error estimate on $Q$ with respect to the iterations of the adaptive process. As the iterations proceed, the error (which is initially close to 25%, due purposely to the coarse mesh that is used) monotically decreases until reaching the prescribed tolerance of 1% in about 10 iterations. We also show on Fig. 20 the evolution of the error indicators along the adaptive process. We observe that, for this case, the fine mesh is first refined before adapting the coarse mesh and finally extending the oversampling zone.

6.1.2 Exponential loading

With the aim to consider a quantity of interest that depends on $u^{\varepsilon}$ (and not on the flux $\boldsymbol{q}^{\varepsilon}$ as in the previous example), we consider here a different loading: $f(x_{1},x_{2})=\exp(-(x_{1}^{2}+x_{2}^{2}))$. We again take an initial mesh $\mathcal{T}_{H}$ made of $9\times 9$ macro elements, with $h_{K}=\varepsilon/3$ for any $K$ and without any oversampling. The reference solution is computed on a $500\times 500$ fine mesh. The MsFEM solution $u^{\varepsilon}_{H}$ is compared to the reference solution $u^{\varepsilon}$ on Fig. 21, their gradients $\boldsymbol{\nabla}u^{\varepsilon}_{H}$ and $\boldsymbol{\nabla}u^{\varepsilon}$ on Fig. 22, and fluxes $\boldsymbol{q}^{\varepsilon}_{H}$ and $\boldsymbol{q}^{\varepsilon}$ on Fig. 23. We also represent the distribution of the energy of the exact solution and of the global error on Fig. 24. In contrast to the previous test-case, $u^{\varepsilon}$ is not close to 0 in the zone where the defect is located (see Fig. 21). This is thus a relevant test-case for considering as quantity of interest the average of $u^{\varepsilon}$ in that zone.

We also plot on Fig. 25 the distribution of the error on $u^{\varepsilon}$, $\boldsymbol{\nabla}u^{\varepsilon}\cdot\boldsymbol{e}_{1}$ and $\boldsymbol{q}^{\varepsilon}\cdot\boldsymbol{e}_{1}$.

As pointed out above, we choose the quantity of interest $\displaystyle Q(u^{\varepsilon})=\frac{1}{|\omega|}\int_{\omega}u^{\varepsilon}$ with $|\omega|=4\varepsilon\times 4\varepsilon$. We consider two cases:

Case 1: the subregion $\omega$ is centered on the defect ($Q=Q_{1}$);

• 2.

Case 2: the subregion $\omega$ is far from the defect ($Q=Q_{2}$).

The corresponding adjoint solutions $\widetilde{u}^{\varepsilon}$ and fluxes $\widetilde{\boldsymbol{q}}^{\varepsilon}$ are shown on Figs. 26 and 27. Again, as expected, the adjoint solutions are essentially supported in the domain of interest $\omega$.

The error threshold is fixed to $1\%$. We show on Fig. 28 (resp. Fig. 29) the final MsFEM discretization obtained after using the adaptive algorithm presented in Section 5.3, in order to respect the error threshold for the quantity of interest $Q_{1}$ (resp. $Q_{2}$). This discretization can be compared with that obtained when controlling the global error (in energy norm) using the tools presented in [19], which is shown on Fig. 30. The conclusions are the same as in the case of sinusoidal loading (see Section 6.1.1): (i) it is worth working with a goal-oriented adaptive strategy rather than with an adaptive strategy based on the global error and (ii) the mesh needs to be refined only in the region close to the region of interest $\omega$.

As regards the accuracy of the error estimate, we again observe that it is very good with only a slight overestimation of the error. Computations yield $\eta^{Q}/|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})|=1.14$ (resp. 1.15) at the first iteration of the adaptive algorithm, and $\eta^{Q}/|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})|=1.07$ (resp. 1.04) at the final iteration of the adaptive algorithm, for the quantities of interest $Q_{1}$ and $Q_{2}$, respectively.

We show on Fig. 31 how the error estimate on $Q$ decreases as the adaptive iterations of the mesh proceed. For instance, after the first adaptation, the error has been divided by a factor 2, going from 80% to 50% for $Q_{1}$ (resp. 40% to 20% for $Q_{2}$).

6.2 Flow problem

We now consider a flow problem in a fractured porous medium, with a source and a sink. This problem has already been considered in [20], with different multiscale techniques. The conductivity within the fractures can be several orders of magnitude higher than the conductivity within the matrix. In flow applications, one needs to obtain a good approximation of the pressure in a neighborhood of the wells.

We consider the domain $\Omega=(0,1)^{2}$ and the inflow-outflow source term $f=\chi_{K_{1}}-\chi_{K_{2}}$ , where $K_{1}=[0.1,0.2]\times[0.8,0.9]$ and $K_{2}=[0.8,0.9]\times[0.1,0.2]$. The pressure $u^{\varepsilon}$ is solution to

[TABLE]

and we consider the quantity of interest $\displaystyle Q(u^{\varepsilon})=\frac{1}{|K_{2}|}\int_{K_{2}}u^{\varepsilon}$, that is the average value of the solution (i.e. the pressure) on the outflow region $K_{2}$ (where the well is located).

The scalar-valued high-contrast function $A^{\varepsilon}(\boldsymbol{x})$ is shown on Fig. 32. We observe a high-conductivity channel crossing the domain $\Omega$, and separating the inflow and outflow regions. In the blue region, we have $A^{\varepsilon}=1$, and in the red region we set $A^{\varepsilon}=10^{4}$. The contrast between the conductivities in the blue and red regions is thus very large.

We take an initial mesh $\mathcal{T}_{H}$ made of $10\times 10$ macro elements, with $h_{K}=H/10$ for any $K$ and no oversampling. We show on Figs 33, 34 and 35 the exact solution $u^{\varepsilon}$ and its MsFEM approximation $u^{\varepsilon}_{H}$, their gradients $\boldsymbol{\nabla}u^{\varepsilon}$ and $\boldsymbol{\nabla}u^{\varepsilon}_{H}$, and the fluxes $\boldsymbol{q}^{\varepsilon}$ and $\boldsymbol{q}^{\varepsilon}_{H}$. We can see that the solution $u^{\varepsilon}$ is essentially supported in the source region $K_{1}$ and the sink region $K_{2}$. The initial approximation $u^{\varepsilon}_{H}$ is not very accurate, especially in the top half of $\Omega$. As can be seen on Fig. 40, the initial error on the quantity of interest is almost 60%. The discretization is purposely too coarse.

The error between the exact solution and its MsFEM approximation is shown on Fig. 36.

The adjoint solution $\widetilde{u}^{\varepsilon}$ and flux $\widetilde{\boldsymbol{q}}^{\varepsilon}$ corresponding to the quantity of interest $Q$ are shown on Fig. 37.

The error threshold is fixed to $1\%$. We show on Fig. 38 the final MsFEM discretization that we obtain in order to respect the error threshold for the quantity of interest. As in the previous test cases, this discretization is much coarser than the discretization (shown on Fig. 39) that is obtained when controlling the global error (in energy norm) using tools presented in [19].

The accuracy of the error estimate is again very good in this example. We get $\eta^{Q}/|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})|=1.14$ at the first iteration of the adaptive algorithm, and $\eta^{Q}/|Q(u^{\varepsilon})-Q(u^{\varepsilon}_{H})|=1.08$ at the final iteration of the adaptive algorithm.

The error estimate on the quantity of interest is shown on Fig. 40, as the iterations of the adaptive process proceed. We again observe that this error monotically decreases.

7 Conclusions and perspectives

We have introduced efficient numerical tools for the a posteriori error estimation of specific quantities of interest in the context of MsFEM computations. These tools are based on the Constitutive Relation Error concept, that we have previously used in [19] to estimate the error in the energy norm, and that we have here successfully extended to the goal-oriented context.

Our approach leads to fully computable and guaranteed upper bounds on the error, and enables to perform adaptive multiscale computations by selecting the relevant MsFEM parameters to reach a prescribed error tolerance. Our procedure is consistent with the offline/online paradigm of MsFEM, hence yielding affordable additional computational costs.

The performances of the approach have been investigated on several multiscale problems. Fine scale features of the solution are adaptively recovered during the iterations, if and where needed. In particular, the fine-scale features of the solution in some parts of the physical domain may not affect the quantity of interest, and our procedure then prevents the discretization approach to capture them, thus achieving a significant computational gain at no loss in accuracy. Our method thus outperforms the standard a posteriori error estimation strategy (which assesses and controls the error in the energy norm): for an equal tolerance on the error, a coarser discretization is obtained, leading to a smaller computational cost.

Acknowledgements

We thank C. Le Bris for stimulating discussions on the topics addressed in this article. LC thanks Inria for enabling his two-year leave (2014–2016) in the MATHERIALS project team. The work of FL is partially supported by ONR under grant N00014-15-1-2777 and by EOARD under grant FA9550-17-1-0294.

Références

• [1]

Abdulle A., Nonnenmacher A.

A posteriori error analysis of the heterogeneous multiscale method for homogenization problems.

C.R. Acad. Sci. Paris, Série I 2009; 347:1081–1086.

• [2]

Abdulle A., Nonnenmacher A.

A posteriori error estimates in quantities of interest for the finite element heterogeneous multiscale method.

Numerical Methods for Partial Differential Equations 2013; 29(5):1629–1656.

• [3]

Ainsworth M., Oden J.T.

A posteriori error estimation in finite element analysis.

John Wiley & Sons 2000.

• [4]

Akbari Rahimabadi A., Kerfriden P., Bordas S.

Scale selection in nonlinear fracture mechanics of heterogeneous materials.

Philosophical Magazine 2015; 95(28-30):3328–3347.

• [5]

Allaire G., Brizzi R.

A multiscale finite element method for numerical homogenization.

SIAM Multiscale Modeling & Simulation 2006; 4(3):790–812.

• [6]

Becker R., Rannacher R.

A feed-back approach to error control in finite element methods: Basic analysis and examples.

East-West Journal of Numerical Mathematics 1996; 4:237–264.

• [7]

Becker R., Rannacher R.

An optimal control approach to shape a posteriori error estimation in finite element methods.

Acta Numerica, Cambridge University Press, 2001; 10:1–120.

• [8]

Ben Dhia H., Rateau G.

Mathematical analysis of the mixed Arlequin method.

C.R. Acad. Sci. Paris 2001; 332(7):649–654.

• [9]

Bensoussan A., Lions J.-L., Papanicolaou G.

Asymptotic Analysis for Periodic Structures.

Studies in Mathematics and its Applications, vol. 5, North-Holland, Amsterdam, New York 1978.

• [10]

Blanc X., Le Bris C., Lions P.-L.

A possible homogenization approach for the numerical simulation of periodic microstructures with defects.

Milan Journal of Mathematics 2012; 80(2):351–367.

• [11]

Blanc X., Le Bris C., Lions P.-L.

Local profiles and elliptic problems at different scales with defects.

C.R. Acad. Sci. Paris, Série I 2015; 353:203–208.

• [12]

Blanc X., Le Bris C., Lions P.-L.

Local profiles for elliptic problems at different scales: defects in, and interfaces between periodic structures.

Communications in Partial Differential Equations 2015; 40(12):2173–2236.

• [13]

Blanc X., Le Bris C., Lions P.-L.,

On correctors for linear elliptic homogenization in the presence of local defects.

Communications in Partial Differential Equations 2018; 43(6):965–997.

• [14]

Blanc X., Le Bris C., Lions P.-L.,

On correctors for linear elliptic homogenization in the presence of local defects: the case of advection-diffusion.

Journal de Mathématiques Pures et Appliquées 2019; 124:106–122.

• [15]

Brezzi F., Lions J.-L., Pironneau O.

Analysis of a chimera method.

C.R. Acad. Sci. Paris 2001; 332:655–660.

• [16]

Cao T., Kelly D.W.

Pointwise and local error estimates for the quantities of interest in two-dimensional elasticity.

Computers and Mathematics with Applications 2003; 46(1):69–79.

• [17]

Chamoin L., Ladevèze P.

A non-intrusive method for the calculation of strict and efficient bounds of calculated outputs of interest in linear viscoelasticity problems.

Computer Methods in Applied Mechanics and Engineering 2008; 197(9-12):994–1014.

• [18]

Chamoin L., Legoll F.

A pedagogical view on a posteriori error estimation in finite element analysis.

in preparation.

• [19]

Chamoin L., Legoll F.

A posteriori error estimation and adaptive strategy for the control of MsFEM computations.

Computer Methods in Applied Mechanics and Engineering 2018; 336:1–38.

• [20]

Chung E.T., Leung W.T., Pollock S.

Goal-oriented adaptivity for GMsFEM.

Journal of Computational and Applied Mathematics 2016; 296:625–637.

• [21]

Cirak F., Ramm E.

A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem.

Computer Methods in Applied Mechanics and Engineering 1998; 156:351–362.

• [22]

E W., Engquist B., Huang Z.

Heterogeneous multiscale method: a general methodology for multiscale modeling.

Physical Review B 2003; 67(9):092101.

• [23]

Efendiev Y., Galvis J., Hou T.

Generalized multiscale finite element methods.

Journal of Computational Physics 2013; 251:116–135.

• [24]

Efendiev Y., Hou T.

Multiscale Finite Element Methods: Theory and Applications.

Springer, New York 2009.

• [25]

Efendiev Y., Hou T., Wu X.H.

Convergence of a nonconforming multiscale finite element method.

SIAM Journal of Numerical Analysis 2000; 37(3):888–910.

• [26]

Feyel F.

A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua.

Computer Methods in Applied Mechanics and Engineering 2003; 192:3233–3244.

• [27]

Fraeijs de Veubeke B., Hugge M.A.

Dual analysis for heat conduction problems by finite elements.

International Journal for Numerical Methods in Engineering 1972; 5(1):65–82.

• [28]

Gendre L., Allix O., Gosselet P., Comte F.

Non-intrusive and exact global/local techniques for structural problems with local plasticity.

Computational Mechanics 2009; 44(2):233–245.

• [29]

Giles M.B., Suli E.

Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality.

Acta Numerica 2002:145–236.

• [30]

Grätsch T., Hartmann F.

Finite element recovery techniques for local quantities of linear problems using fundamental solutions.

Computational Mechanics 2003; 33:15–21.

• [31]

Henning P., Ohlberger M., Schweizer B.

An adaptive multiscale finite element method.

SIAM Multiscale Modeling & Simulation 2014; 12(3):1078–1107.

• [32]

Hesthaven J., Zhang S., Zhu X.

High-order multiscale finite element method for elliptic problems.

SIAM Multiscale Modeling & Simulation 2014; 12(2):650–666.

• [33]

Hou T., Wu X.H.

A multiscale finite element method for elliptic problems in composite materials and porous media.

Journal of Computational Physics 1997; 134:169–189.

• [34]

Hou T., Wu X.H., Cai Z.

Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients.

Mathematics of Computation 1999; 68(227):913–943.

• [35]

Hou T., Wu X.H., Zhang Y.

Removing the cell resonance error in the multiscale finite element method via a Petrov-Galerkin formulation.

Communications in Mathematical Sciences 2004; 2:185–205.

• [36]

Hughes T.J.R., Feijoo G.R., Mazzei L., Quincy J.B.

The variational multiscale method – A paradigm for computational mechanics.

Computer Methods in Applied Mechanics and Engineering 1998; 166(1-2):3–24.

• [37]

Jhurani C., Demkowicz L.

Multiscale modeling using goal-oriented adaptivity and numerical homogenization. Part 1: mathematical formulation and numerical results.

Computer Methods in Applied Mechanics and Engineering 2012; 213-216:399–417.

• [38]

Ladevèze P.

Strict upper error bounds for calculated outputs of interest in computational structural mechanics.

Computational Mechanics 2008; 42(2):271–286.

• [39]

Ladevèze P., Chamoin L.

Calculation of strict error bounds for finite element approximations of nonlinear pointwise quantities of interest.

International Journal for Numerical Methods in Engineering 2010; 84:1638–1664.

• [40]

Ladevèze P., Leguillon D.

Error estimate procedure in the finite element method and applications.

SIAM Journal of Numerical Analysis 1983; 20(3):485–509.

• [41]

Ladevèze P., Maunder E.A.W.

A general method for recovering equilibrating element tractions.

Computer Methods in Applied Mechanics and Engineering 1996; 137:111–151.

• [42]

Ladevèze P., Pelle J.-P.

Mastering Calculations in Linear and Nonlinear Mechanics.

Springer, 2004.

• [43]

Ladevèze P., Pled F., Chamoin L.

New bounding techniques for goal-oriented error estimation applied to linear problems.

International Journal for Numerical Methods in Engineering 2013; 93(13):1345–1380.

• [44]

Larson M.G., Malqvist A.

Adaptive variational multiscale methods based on a posteriori error estimation: duality techniques for elliptic problems.

In B. Engquist, P. Lötstedt, and O. Runborg, editors, Multiscale Methods in Science and Engineering, volume 44 of Lecture Notes in Computational Science and Engineering, pages 181–193. Springer, 2005.

• [45]

Larson M.G., Malqvist A.

Adaptive variational multiscale methods based on a posteriori error estimation: energy norm estimates for elliptic problems.

Computer Methods in Applied Mechanics and Engineering 2007; 196(21-24):2313–2324.

• [46]

Le Bris C., Legoll F., Lozinski A.

An MsFEM type approach for perforated domains.

SIAM Multiscale Modeling & Simulation 2014; 12(3):1046–1077.

• [47]

Le Bris C., Legoll F., Thomines F.

Multiscale finite element approach for “weakly” random problems and related issues.

Mathematical Modelling and Numerical Analysis 2014; 48(3):815–858.

• [48]

Lions J.-L., Pironneau O.

Domain decomposition methods for CAD.

C.R. Acad. Sci. Paris, Série I 1999; 328(1):73–80.

• [49]

Lozinski A., Pironneau O.

Numerical zoom for localized multiscales.

Numerical Methods for Partial Differential Equations 2011; 27:197–207.

• [50]

Maier M., Rannacher R.

Duality-based adaptivity in finite element discretization of heterogeneous multiscale problems.

Journal of Numerical Mathematics 2016; 24(3):167–187.

• [51]

Malqvist A., Peterseim D.

Localization of elliptic multiscale problems.

Mathematics of Computation 2014; 83(290):2583–2603.

• [52]

Nolen J., Papanicolaou G., Pironneau O.

A framework for adaptive multiscale methods for elliptic problems.

SIAM Multiscale Modeling & Simulation 2008; 7(1):171–196.

• [53]

Oden J.T., Prudhomme S.

Goal-oriented error estimation and adaptivity for the finite element method.

Computers and Mathematics with Applications 2001; 41:735–756.

• [54]

Oden J.T., Prudhomme S.

Estimation of modeling error in Computational Mechanics.

Journal of Computational Physics 2002; 182:496–515.

• [55]

Oden J.T., Reddy J.N.

On dual complementary variational principles in mathematical physics.

International Journal of Engineering Science 1974; 12:1–29.

• [56]

Oden J.T., Vemaganti K.

Estimation of local modeling error and goal-oriented modeling of heterogeneous materials; Part 1: Error estimates and adaptive algorithms.

Journal of Computational Physics 2000; 164:22-47.

• [57]

Oden J.T., Vemaganti K., Moës N.

Hierarchical modeling of heterogeneous solids.

Computer Methods in Applied Mechanics and Engineering 1999; 172(1-4):3–25.

• [58]

Panetier J., Ladevèze P., Chamoin L.

Strict and effective bounds in goal-oriented error estimation applied to fracture mechanics problems solved with XFEM.

International Journal for Numerical Methods in Engineering 2010; 81:671–700.

• [59]

Paraschivoiu M., Peraire J., Patera A.T.

A posteriori finite element bounds for linear functional outputs of elliptic partial differential equations.

Computer Methods in Applied Mechanics and Engineering 1997; 150:289–312.

• [60]

Pled F., Chamoin L., Ladevèze P.

On the techniques for constructing admissible stress fields in model verification: performances on engineering examples.

International Journal for Numerical Methods in Engineering 2011; 88(5):409–441.

• [61]

Prudhomme S., Chamoin L., Ben Dhia H., Bauman P.T.

An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations.

Computer Methods in Applied Mechanics and Engineering 2009; 198(21-26):1887–1901.

• [62]

Prudhomme S., Oden J.T.

On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors.

Computer Methods in Applied Mechanics and Engineering 1999; 176:313–331.

• [63]

Rannacher R., Suttmeier F.T.

A feedback approach to error control in finite element methods: application to linear elasticity.

Computational Mechanics 1997; 19:434–446.

• [64]

Strouboulis T., Copps K., Babuska I.

The generalized finite element method.

Computer Methods in Applied Mechanics and Engineering 2001; 190:4081–4193.

• [65]

Strouboulis T., Zhang L., Wang D., Babuska I.

A posteriori error estimation for generalized finite element methods.

Computer Methods in Applied Mechanics and Engineering 2006; 195:852–879.

• [66]

Vemaganti K., Oden J.T.

Estimation of local modeling error and goal-oriented modeling of heterogeneous materials; Part 2: A computational environment for adaptive modeling of heterogeneous elastic solids.

Computer Methods in Applied Mechanics and Engineering 2001; 190:6089-6124.

• [67]

Verfürth R.

A review of a posteriori error estimates and adaptive mesh-refinement techniques.

Wiley-Teubner 1996.

• [68]

Waeytens J., Chamoin L., Ladevèze P.

Guaranteed error bounds on pointwise quantities of interest for transient viscodynamics problems.

Computational Mechanics 2012; 49:291–307.