Numerical upscaling of perturbed diffusion problems

TL;DR

This paper introduces a numerical method for efficiently solving perturbed elliptic PDEs with rapidly varying coefficients by reusing local computations, leveraging PG-LOD for parallelization and low memory use.

Contribution

The paper develops a PG-LOD based approach for perturbed diffusion problems, enabling efficient reuse of computations for local defects and global mappings.

Findings

Method achieves efficient solutions for perturbed problems.

Parallelization reduces computational overhead.

Numerical examples demonstrate accuracy and efficiency.

Abstract

In this paper we study elliptic partial differential equations with rapidly varying diffusion coefficient that can be represented as a perturbation of a reference coefficient. We develop a numerical method for efficiently solving multiple perturbed problems by reusing local computations performed with the reference coefficient. The proposed method is based on the Petrov--Galerkin Localized Orthogonal Decomposition (PG-LOD) which allows for straightforward parallelization with low communcation overhead and memory consumption. We focus on two types of perturbations: local defects which we treat by recomputation of multiscale shape functions and global mappings of a reference coefficient for which we apply the domain mapping method. We analyze the proposed method for these problem classes and present several numerical examples.

Equations

$$-\nabla \cdot A\nabla \bar{u}$$

$$\bar{u}$$

$$A(x)\approx A_{\text{ref}}(x),f(x)\approx f_{\text{ref}}(x).$$

$$0<\alpha$$

$$\infty >\beta$$

$$V:=V_h\subset H_0^1(\Omega )=\{v\in H^1(\Omega )|\text{tr}(v)=0\},$$

$$a(u,v)=F(v)-a(g,v)$$

$$a(u,v):={\int }_{\Omega }\left(A\nabla u\right)\cdot \nabla v,F(v):={\int }_{\Omega }fv.$$

$${\int }_{\Omega }{A}_{\text{ref}}\nabla (u_{\text{ref}}+g)\cdot \nabla vdx={\int }_{\Omega }{f}_{\text{ref}}v.$$

$$\begin{array}{rl}{|||u-u_{\text{ref}}|||}^2& \le (A\nabla u-A_{\text{ref}}\nabla u_{\text{ref}},\nabla (u-u_{\text{ref}}))+(A_{\text{ref}}\nabla u_{\text{ref}}-A\nabla u_{\text{ref}},\nabla (u-u_{\text{ref}}))\\ & =(f-f_{\text{ref}},u-u_{\text{ref}})+((A_{\text{ref}}-A)\nabla (u_{\text{ref}}+g),\nabla (u-u_{\text{ref}}))\\ & \le \left(\frac{C_{\text{P}}}{{\alpha }^{1/2}}\Vert f-f_{\text{ref}}{\Vert }_{L^2(\Omega )}+\frac{C_{\text{P}}}{{\alpha }^{3/2}}\Vert A_{\text{ref}}-A{\Vert }_{L^{\infty }(\Omega )}\Vert f_{\text{ref}}{\Vert }_{L^2(\Omega )}\right)|||u-u_{\text{ref}}|||,\end{array}$$

$$A_{\text{ref}}={\chi }_{{\Omega }_1}+\alpha {\chi }_{{\Omega }_{\alpha }},$$

$${\mathcal{J}}_{ij}(x)=\left[\frac{\partial {\psi }_i}{\partial x_j}\right],$$

$$\begin{array}{rl}-{\nabla }_y\cdot A_y{\nabla }_y{u}_y& =f_y,\text{ in }\Omega ,\\ u_y& =g_y,\text{ on }\Gamma ,\end{array}$$

$${\nabla }_{y}v(x):={\left[\frac{\partial v}{\partial y_i}(x)\right]}_i={\left[\sum _j\frac{\partial v}{\partial x_j}(x)\frac{\partial x_j}{\partial y_i}(x)\right]}_i={\mathcal{J}}^{-T}(x){\nabla }_{x}v(x).$$

$$\begin{array}{rl}a(v,w)& ={\int }_{\Omega }\left((A_{\text{ref}}-D)\circ {\psi }^{-1}\right){\nabla }_y\left(v\circ {\psi }^{-1}\right)\cdot {\nabla }_y\left(w\circ {\psi }^{-1}\right)\text{d}y\\ & ={\int }_{\Omega }\det (\mathcal{J}){\mathcal{J}}^{-1}(A_{\text{ref}}-D){\mathcal{J}}^{-T}{\nabla }_{x}v\cdot {\nabla }_{x}w\text{d}x\end{array}$$

$$F(w)={\int }_{\Omega }{f}_y\left(w\circ {\psi }^{-1}\right)\text{d}y={\int }_{\Omega }\det (\mathcal{J})\left(f_y\circ \psi \right)w\text{d}x.$$

$$A=\det (\mathcal{J}){\mathcal{J}}^{-1}(A_{\text{ref}}-D){\mathcal{J}}^{-T},f=\det (\mathcal{J})\left(f_y\circ \psi \right),g{|}_{\Gamma }=g_y\circ \psi ,$$

$$V_H:=V\cap {\mathcal{P}}_1({\mathcal{T}}_H),$$

$$U_0(\omega )$$

$$U_{k+1}(\omega )$$

$${\mathcal{I}}_{H}v:=\sum _{x\in \mathcal{N}}{\zeta }_x{\lambda }_x,$$

$${\int }_T(P_{T}v)w_H={\int }_{T}v{w}_H,$$

$${\zeta }_x=\mathrm{card}({\mathcal{T}}_x{)}^{-1}\sum _{T^{\prime }\in {\mathcal{T}}_x}(P_{T^{\prime }}v)(x).$$

$$H_T^{-1}\Vert v-{\mathcal{I}}_{H}v{\Vert }_{L^2(T)}+\Vert \nabla {\mathcal{I}}_{H}v{\Vert }_{L^2(T)}\le C_{{\mathcal{I}}_H}\Vert \nabla v{\Vert }_{L^2(U_1(T))},$$

$$V^{\text{f}}=\ker ({\mathcal{I}}_H)=\left\{v\in V|{\mathcal{I}}_H(v)=0\right\}$$

$$a(\mathcal{Q}v,v^{\text{f}})=a(v,v^{\text{f}}),$$

$$a(\mathcal{R}f,v^{\text{f}})={\int }_{\Omega }f{v}^{\text{f}}.$$

$$a(u_H+u^{\text{f}},v_H)$$

$$a(u^{\text{f}},v^{\text{f}})$$

$$a(u^{\text{ms}},v_H)=F(v_H)-a(\mathcal{R}f,v_H)-a(g-\mathcal{Q}g,v_H),$$

$$\left\{{\lambda }_x-\mathcal{Q}{\lambda }_x|x\in \mathcal{N}\right\}.$$

$$V^{\text{f}}(U_{k}(\omega)):=\left\{v\in V^{\text{f}}\,\middle|\,v\big{|}_{\Omega\setminus U_{k}(\omega)}=0\right\}.$$

Full text

Abstract. In this paper we study elliptic partial differential equations with rapidly varying diffusion coefficient that can be represented as a perturbation of a reference coefficient. We develop a numerical method for efficiently solving multiple perturbed problems by reusing local computations performed with the reference coefficient. The proposed method is based on the Petrov–Galerkin Localized Orthogonal Decomposition (PG-LOD) which allows for straightforward parallelization with low communcation overhead and memory consumption. We focus on two types of perturbations: local defects which we treat by recomputation of multiscale shape functions and global mappings of a reference coefficient for which we apply the domain mapping method. We analyze the proposed method for these problem classes and present several numerical examples.

1 Introduction

Manufactured heterogeneous materials, such as composites with tailored properties, are crucial tools in engineering. The challenge of performing accurate computer simulations involving such materials have driven the development of multiscale methods over decades [16, 25, 15, 20, 22]. Multiscale methods have turned out to be successful in computing coarse-scale representations of the solutions to such problems. However, when the heterogeneous data is perturbed it is not obvious how multiscale methods can be adapted. Understanding the effect of perturbations is important since manufactured materials, in general, will not be perfect. Manufacturing tolerances and faults lead to perturbations in the material distribution. There are also other problems, such as time stepping with time dependent diffusion coefficient, optimization of material distribution, and non-linear diffusion problems, which call for iterative procedures where the data in the current iterate can be seen as a perturbation of the data in the previous.

In this paper, we study elliptic problems with diffusion coefficients that are perturbations of a single reference diffusion coefficient. We consider the following Dirichlet type problem, which we will refer to as the strong form of the (inhomogeneous) perturbed problem: find $\bar{u}$ such that

[TABLE]

on a bounded polygonal/polyhedral domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, with boundary $\Gamma$. We assume the right hand side $f\in L^{2}(\Omega)$, diffusion coefficient $A\in L^{\infty}(\Omega,\mathbb{R}^{d\times d})$ is symmetric positive definite and rapidly varying, and the trace of the function $g\in H^{1}(\Omega)$ defines the Dirichlet boundary conditions. We further consider $A$ and $f$ to be a perturbations of a reference diffusion coefficient and a reference right hand side, respectively,

[TABLE]

The aim of the paper is to reuse computations made with the reference quantities $A_{\text{ref}}$ and $f_{\text{ref}}$ in a reliable way when solving a set of perturbed problems. The perturbations we study are local defects and domain mapppings. Figure 1.1 illustrates the reference coefficient (left), random defects (center), and domain mapping to the physical diffusion coefficient (right).

Perturbations of the diffusion coefficient in elliptic problems have been studied extensively. This work was inspired by a series of papers by Le Bris and coworkers on weakly random homogenization [18, 19], where they consider weakly random coefficient problems, similar to the once illustrated in Figure 1.1, using the multiscale finite element method [15] for the spatial discretization. This allows the authors to consider both rapidly varying and perturbed diffusion coefficients. There are several works on partial differential equations posed on random domains using domain mapping including [8, 3]. In the context of perturbed coefficients, which is considered in this paper, domain mapping is instead used to transform a reference coefficient to the perturbed coefficient. A perturbation in position of the material distribution is transformed back to a difference in value of the reference diffusion coefficient, which is advantageous from a numerical perspective.

Multiscale methods have been a vibrant area of research for decades [15, 25, 22, 1]. The main idea is to solve local fine scale problems to compute an improved basis which is used to solve a global coarse-scale problem. These techniques are often parallel by construction. One method that has proven to give accurate results also for non-periodic diffusion coefficient is the Localized Orthogonal Decomposition (LOD) method [20]. LOD is based on an orthogonal split of the solution space, using the scalar product induced by the weak form of the problem (1.1). In recent years it has been improved and reformulated. In [5] a Petrov–Galerkin version of the method was presented and analyzed. PG-LOD has the advantage that the assembly of the modified stiffness matrix is much faster than for the original method. Concerning the implementation of the LOD, a detailed algebraic overview has been given in [6]. In the recent work [12] a sequence of problems with similar coefficients are considered, with applications in time dependent diffusion problems. This approach is also useful for studying perturbations of the type presented in Figure 1.1, see [17].

In this paper we apply the PG-LOD methodology introduced in [12] to solve elliptic problems with perturbed diffusion coefficient. The PG-LOD method allows for local recomputation of basis functions to handle perturbations in the data from defects and domain mappings. We derive error indicators to decide where recomputations are necessary. In this way we can efficiently simulate a vast number of perturbations by mainly solving (upscaled) coarse-scale problems.

The paper is organized as follows. In Section 2 we formulate the problem and the types of perturbations that are considered in this paper. In Section 3 we present the proposed numerical method based on PG-LOD. In Section 4 we derive error bounds and in Section 5 we discuss implementation, memory consumption and parallelization. Finally in Section 6 we present numerical examples.

2 Problem formulation

We assume that the coefficient matrix $A\in L^{\infty}(\Omega,\mathbb{R}^{d\times d})$ is symmetric and uniformly elliptic such that

[TABLE]

We further let $f\in L^{2}(\Omega)$ and $g\in H^{1}(\Omega)$ with a trace on $\Gamma$ which defines the Dirichlet boundary values. We introduce a function space $V$ where we seek a solution of equation (1.1) posed on variational form. In this paper we primarily consider a conforming finite element space

[TABLE]

defined on a computational mesh $\mathcal{T}_{h}$ that is assumed to be fine enough to resolve the variations in the diffusion coefficients well. However, we may as well choose $V=H^{1}_{0}(\Omega)$ and the analysis presented in the paper will still go through. To simplify the notation we stick to the notation $u\in V=V_{h}$. On weak form we get: find $u\in V$ such that

[TABLE]

for all $v\in V$, where

[TABLE]

The Lax–Milgram Lemma guarantees existence and uniqueness of a solution $u\in V$. The full solution including the boundary data is given by $u+g$. As mentioned above, (2.3) will be referred to as the perturbed problem.

To motivate why we cannot simply replace the perturbed coefficient with the reference coefficient, we formulate an artificial problem based on the reference coefficient and right hand side: find $u_{\text{ref}}\in V$, such that for all $v\in V$,

[TABLE]

The error between $u_{\text{ref}}$ and $u$ can be bounded in the energy norm ${|\kern-1.29167pt|\kern-1.29167pt|\cdot|\kern-1.29167pt|\kern-1.29167pt|}:=a(\cdot,\cdot)^{1/2}=\lVert A^{1/2}\nabla\cdot\rVert_{L^{2}(\Omega)}$ in the following way,

[TABLE]

where $C_{\text{P}}$ is the Poincaré constant for $\Omega$. This error bound suggests that even local perturbations in the structure of the coefficient or the right hand side, e.g. by a defect or shift, may lead to very poor accuracy. This occurs for example if we consider a problem with a highly conductive thin channel in the diffusion coefficient and a right hand side $f$ which has support inside the channel. If the channel is moved slightly so that the support of the right hand side is now outside the channel the solution will be have very differently. The error with respect to perturbations in $f$ is less severe since it is measured in the $L^{2}$-norm.

For this reason it is not clear how to reuse computations for the standard finite element method in a reliable way. In this paper we will treat this difficulty using a multiscale approach where solutions to localized subproblems, based on the reference coefficient, can be reused when solving for the perturbed problem.

2.1 Perturbations

To simplify the presentation we consider perturbations of coefficients that only takes the two values $1$ and $0<\alpha<1$. We emphasize that this is not necessary for the proposed method to work or for theory to hold. However, it highlights the application to composite materials which has inspired this work. Let $\Omega_{1},\Omega_{\alpha}\subseteq\Omega$ be two disjoint subdomains of $\Omega$ with $\Omega_{1}\cup\Omega_{\alpha}=\Omega$. Let $A_{\text{ref}}$ be defined by

[TABLE]

where $\chi$ is an indicator function.

We formalize the two types of perturbations that we consider (see Figure 1.1) by introducing a defect perturbation $D$ and a domain mapping $\psi$. A perturbation from a defect can be expressed by D=(1-\alpha){\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}_{\omega} where $\omega\subseteq\Omega_{1}$ and that $A_{\text{ref}}-D$ can be considered as the perturbed coefficient. For shift perturbations, we assume that the domain mapping perturbation can be described as a variable transformation with a perturbation function $\psi:\Omega\rightarrow\Omega$ which maps the reference coefficient (expressed in $x$-coordinates) to a mapped coefficient (expressed in $y$-coordinates). We assume that $\psi$ maps the boundary to itself (i.e. $\Gamma=\{\psi(x)\,:\,x\in\Gamma\}$) and that it is a one-to-one mapping in $\Omega$. Figure 2.1 provides an example of a variable transformation. We denote the corresponding Jacobi matrix

[TABLE]

for $1\leq i,j\leq d$, and assume it to be bounded with bounded inverse for a.e. $x\in\Omega$. The two perturbation types can either be combined, or be considered individually by letting $D\equiv 0$ or $\psi=\operatorname{Id}$.

Using $A_{\text{ref}}$, $D$ and $\psi$, we formulate the mapped problem in the $y$-variable with $y=\psi(x)$ by

[TABLE]

where the coefficient is defined by $A_{y}=(A_{\text{ref}}-D)\circ\psi^{-1}$, and the derivatives have been distorted accordingly. The $y$-variable corresponds to the physical spatial variable in a typical situation. Depending on the physics being modeled, either the mapped right hand side $f_{y}\in L^{2}(\Omega)$ or the perturbed $f$ (below) can be considered given. It makes no difference for the development of the numerical method, but may affect the choice of $f_{\text{ref}}$. The Dirichlet boundary value $g_{y}\in H^{1/2}(\Gamma)$ is defined only on the boundary $\Gamma$, which is mapped to itself. The solution in the perturbed domain is denoted $u_{y}(y)$.

Next, we use the mapped problem to define the perturbed problem in equation (2.3). The gradient operator $\nabla_{y}$ in the mapped domain can be expressed in $\nabla_{x}=\nabla$ by

[TABLE]

Based on the elliptic operator in equation (2.7) we define the perturbed bilinear form for the mapped problem as

[TABLE]

and the corresponding linear functional

[TABLE]

We see that this now fits the formulation of the perturbed problem (2.3) with

[TABLE]

and mapped solution $u_{y}=(u+g)\circ\psi^{-1}$. For the problem to be well posed we assume $\mathcal{J}$ and $\mathcal{J}^{-1}$ to be bounded almost everywhere and $A$ to be symmetric positive definite. We note that the perturbed coefficient $A$ can be computed from the reference coefficient by means of the Jacobian matrix and that the Dirichlet boundary value function $g\in H^{1}(\Omega)$ can be chosen arbitrarily in the interior of $\Omega$. We emphasize that the domain mapping transforms a shift defect to a change-in-value perturbation. For many coefficients this is advantageous as seen in equation (2.6), where now the $L^{\infty}$-norm can be expressed entirely in terms of how much $\mathcal{J}$ differs from identity. The domain mapping covers continuous (possibly global) perturbations while defects cover discontinuous (often local) perturbations. With respect to Figure 1.1 we see that the middle picture corresponds to $\psi=\operatorname{Id}$ and the right picture to $D\equiv 0$.

Remark 2.1** (Discretization of $\psi$)**

In the numerical examples in Section 6 we let $V=V_{h}$ be the space of quadrilateral finite elements and $\psi$ to be a linear combination of the corresponding bilinear shape functions. This leads to an isoparametric finite element formulation. Isoparametric finite elements are for instance used to guarantee accuracy when solving problems on curved domains. Here the aim is instead to map the perturbed diffusion coefficient to a reference. The theoretical justification for using bilinear domain mappings $\psi$ follows directly from the theory of isoparametric finite elements, see [4, 2].

3 Adaptive PG-LOD method

In this section, we develop the proposed method that reuses reference computations peformed using the reference coefficient $A_{\text{ref}}$ and right hand side $f_{\text{ref}}$ in an adaptive manner to solve the perturbed problem (2.3) expressed in $A$, $u$, $f$, and $g$. The perturbed data can stem from a domain mapping, defects or both. We consider the Petrov–Galerkin version of the LOD method as presented in [5] and briefly derive it in Sections 3.1–3.3. Section 3.4 defines the computable error indicators and Section 3.5 presents the method that reuses reference computations while balancing the error. Recently, a similar approach was applied to problems with time dependent diffusion, see [12].

3.1 Preliminaries

Let ${\mathcal{T}}_{H}$ be a coarse, shape regular, conforming mesh family of the domain $\Omega$. We denote the maximum diameter of an element in ${\mathcal{T}}_{H}$ with $H$ and $\mathcal{N}$ the set of all corresponding interior nodes of the mesh ${\mathcal{T}}_{H}$. We let

[TABLE]

where $\mathcal{P}_{1}({\mathcal{T}}_{H})$ denotes the space of ${\mathcal{T}}_{H}$-piecewise affine functions that are continuous on the domain $\Omega$. Since our full space $V$ is a conforming finite element space we assume that the meshes and spaces are nested $V_{H}\subset V$. However, it is possible, with a minor modification of the proposed method, to violate this condition and still get convergence, see [21].

We introduce the concept of element patches as they will be used in the definition of the interpolation and for the localization of the PG-LOD method. For arbitrary $\omega\subseteq\Omega$ and $0\leq k\in\mathbb{N}$, we define coarse grid patches $U_{k}(\omega)\subset\Omega$ by

[TABLE]

If $\omega=\left\{x\right\}$, for a node $x\in\mathcal{N}$, we call $U_{k}(x)$ a $k$-layer nodal patch. For $\omega=T$, where $T\in{\mathcal{T}}_{H}$, we call $U_{k}(T)$ a $k$-layer element patch, see Figure 3.1

3.2 Multiscale decomposition

Fine scale features that occur in the solutions are not captured in the space $V_{H}$. We characterize the fine scale parts of $V$ as the kernel of an linear surjective (quasi-)interpolation operator $\mathcal{I}_{H}:V\to V_{H}$ that maps a function $v\in V$ to a function $v_{H}\in V_{H}$ in the coarse FE space. Let ${\mathcal{T}}_{x}=\{T\in{\mathcal{T}}_{H}\,:\,x\in\overline{T}\}$ be the set of elements neigboring $x$. We let the interpolation operator $\mathcal{I}_{H}:V\to V_{H}$ be an $L^{2}$-projection to the broken finite element space in composition with averaging in the nodes. In other words,

[TABLE]

where coefficients $\zeta_{x}$ are determined as follows. Let the operator $P_{T}$ be the element $L^{2}$-projection P_{T}v\in V_{H}\big{|}_{T} such that

[TABLE]

for all w_{H}\in V_{H}\big{|}_{T}, then

[TABLE]

This interpolation operator is linear, continuous and its restriction to $V_{H}$ is an isomorphism. Furthermore it fulfills the stability result

[TABLE]

for every $v\in V$ and $T\in{\mathcal{T}}_{H}$, with a generic constant $C_{\mathcal{I}_{H}}>0$, see [23]. We refer to [11, 24] for other possible choices of interpolation operators.

We let the kernel of $\mathcal{I}_{H}$,

[TABLE]

define the fine scales of the space $V$. Since $\mathcal{I}_{H}$ is a projection, this allows for the split $V=V_{H}\oplus V^{\text{f}}$. We define a correction operator ${\mathcal{Q}}v\in V^{\text{f}}$, for a given $v\in V$, to be the solution of

[TABLE]

for all $v^{\text{f}}\in V^{\text{f}}$, and define the multiscale space $V^{\text{ms}}:=V_{H}-{\mathcal{Q}}V_{H}$.

For any $v^{\text{f}}\in V^{\text{f}}$ and $v^{\text{ms}}\in V^{\text{ms}}$, we observe $a(v^{\text{ms}},v^{\text{f}})=0$. This leads to the orthogonal decomposition with respect to the $a$-scalar product, $V=V^{\text{ms}}\oplus_{a}V^{\text{f}}$. Right hand correction can be used to improve accuracy in LOD based methods, see e.g. [13, 12]. We define this correction by ${\mathcal{R}}f\in V^{\text{f}}$ such that, for all $v^{\text{f}}\in V^{\text{f}}$,

[TABLE]

We now derive the Petrov–Galerkin LOD method for the perturbed problem, see also [5, 7]. We use $V=V_{H}\oplus V^{\text{f}}$ to decompose (2.3) into two equations

[TABLE]

for all $v_{H}\in V_{H}$ and $v^{\text{f}}\in V^{\text{f}}$.

With the definitions of ${\mathcal{Q}}$ and ${\mathcal{R}}$ we obtain from (3.3) that $u^{\text{f}}=-{\mathcal{Q}}u_{H}+{\mathcal{R}}f-{\mathcal{Q}}g$. Plugging this into (3.2) gives

[TABLE]

for all $v_{H}\in V_{H}$. Hence, solving (3.4) gives the exact solution $u=u^{\text{ms}}+{\mathcal{R}}f-{\mathcal{Q}}g$.

In order to use the multiscale space in a practical implementation, we need a computable basis. Since $V^{\text{ms}}$ and $V_{H}$ have equal dimensions, it suffices to apply the fine scale corrector ${\mathcal{Q}}$ on every single basis function $\lambda_{x}$ of $V_{H}$ to obtain a corrected basis, i.e.

[TABLE]

We note that a global fine scale computation for each node is necessary to compute all ${\mathcal{Q}}\lambda_{x}$ which is computationally expensive. However, Figure 3.2 suggests that the computation can be localized to a small area around the support of the original basis function. It was shown in [20] that the corrected basis functions decay exponentially, and that localized computations are possible.

Remark 3.1** (Mixed boundary conditions)**

It is also possible to have mixed Neumann and Dirichlet boundary conditions in the PG-LOD method, see [13]. With mixed boundary conditions we also need to compute fine scale corrections for the inhomogeneous Neumann boundary data (as we do for the Dirichlet data) to get optimal convergence rate.

3.3 Localized multiscale method

The fine scale space $V^{\text{f}}$ can be restricted to patches $U_{k}(\omega)$ with the intuitive definition, for $\omega\subseteq\Omega$ and $k\in\mathbb{N}$,

[TABLE]

These local fine scale patches enable the truncation of the corrector. We define localized element correction operators ${\mathcal{Q}}_{k}$ and ${\mathcal{R}}_{k}$ by

[TABLE]

where ${\mathcal{R}}_{k,T},{\mathcal{Q}}_{k,T}\,:\,V\to V^{\text{f}}(U_{k}(T))$ solves

[TABLE]

for all $v^{\text{f}}\in V^{\text{f}}(U_{k}(T))$. We construct a localized multiscale space $V_{k}^{\text{ms}}$ using $V_{H}$ and the local correctors

[TABLE]

This space is spanned by $\left\{\lambda_{x}-{\mathcal{Q}}_{k}\lambda_{x}\right\}_{x\in\mathcal{N}}$.

We formulate the localized version of the Petrov–Galerkin LOD in (3.4): find $u^{\text{ms}}_{k}\in V_{k}^{\text{ms}}$ such that for all $v\in V_{H}$,

[TABLE]

where the full approximation of $u$ is

[TABLE]

The main reason for using a Petrov–Galerkin formulation is that it avoids the expensive computation of products between corrected basis functions without losing convergence order [5].

The exponential decay of the correctors yield the following error bounds (see [20]) of the localized correctors in terms of $k$:

[TABLE]

where the notation $a\lesssim b$ means $a\leq Cb$ with a constant $C$ independent of $H$, $k$, and TOL (which is used in the coming sections). The well-posedness of (3.5) was studied previously in [5, 12] and appears to be conditioned on sufficiently large $k$ in general. (This condition will be revisited in the proof of Theorem 4.1). Furthermore, from [12, Section 4.2] we obtain an error bound for $u-u_{k}$ which reads

[TABLE]

3.4 Error indicators

Since the aim of this paper is to reuse local corrector computations performed with the reference coefficients $A_{\text{ref}}$ and $f_{\text{ref}}$, we need a notation for a modified ${\mathcal{Q}}_{k,T}$ and ${\mathcal{R}}_{k,T}$, computed using $A_{\text{ref}}$ instead of $A$ and $f_{\text{ref}}$ instead of $f$. We let ${\mathcal{Q}}^{\text{ref}}_{k,T}\,:\,V\to V^{\text{f}}(U_{k}(T))$ and ${\mathcal{R}}^{\text{ref}}_{k,T}\,:\,L^{2}(\Omega)\to V^{\text{f}}(U_{k}(T))$ be the solutions to

[TABLE]

for all $v^{\text{f}}\in V^{\text{f}}(U_{k}(T))$. In order to decide for which $T\in{\mathcal{T}}_{H}$ we need to recompute the correctors (${\mathcal{Q}}_{k,T}$ and ${\mathcal{R}}_{k,T}$) and for which we can still use the reference correctors (${\mathcal{Q}}^{\text{ref}}_{k,T}$ and ${\mathcal{R}}^{\text{ref}}_{k,T}$) we need computable error indicators. For the indicators to be useful and efficient they have to be computable and independent of ${\mathcal{Q}}_{k,T}$ and ${\mathcal{R}}_{k,T}$.

For every $T\in{\mathcal{T}}_{H}$ we define error indicators only depending on the reference corrector and coefficients $A$ and $A_{\text{ref}}$. The definition of the error indicators is motivated by Lemma 4.2 and by Theorem 4.1 in Section 4 below. See [12] for a similar construction.

Definition 3.2** (Error indicators)**

For each $T\in{\mathcal{T}}_{H}$, we define

[TABLE]

where {\mathchoice{\raisebox{0.0pt}{\displaystyle\chi}}{\raisebox{0.0pt}{\textstyle\chi}}{\raisebox{0.0pt}{\scriptstyle\chi}}{\raisebox{0.0pt}{\scriptscriptstyle\chi}}}_{T} denotes the indicator function for an element $T\in{\mathcal{T}}_{H}$ and

[TABLE]

We define the square root of a symmetric positive definite matrix $A^{1/2}$ as the unique principal square root, also positive definite. The error indicators are defined with the goal to reduce memory consumption. With the definitions above, all quantities depending on the reference coefficient and right hand side can be computed in advance. Further implementation details will be discussed in Section 5.

The numerical method we propose in the next subsection exploits the possibility to replace e.g. ${\mathcal{R}}_{k,T}f$ with the precomputed ${\mathcal{R}}^{\text{ref}}_{k,T}f_{\text{ref}}$, and to reuse integrals for the stiffness matrix based on $A_{\text{ref}}$ instead of using $A$. The error indicators will be used to identify which local problems that need to be recomputed and which local problems that can be computed using the reference data.

Remark 3.3** (Perturbed boundary data)**

We do not cover the case when $g$ is perturbed. It can be perturbed explicitly from a reference boundary condition (in which case a reference $g$ similar to $f_{\text{ref}}$ needs to be introduced), or it can be perturbed by a domain mapping $g|_{\Gamma}=g_{y}\circ\psi$. We remark that if $\psi$ is the identity mapping for all points on the boundary $\Gamma$, then $g$ (including its values in the interior of the domain) can be picked independent of $\psi$.

Remark 3.4** (Pure domain mapping)**

In the domain mapping setting, assuming $A_{\text{ref}}$ to be scalar and $D\equiv 0$, we have that

[TABLE]

i.e. it is independent of $A_{\text{ref}}$. If the Jacobi matrix can be written as a $\epsilon$-perturbation of the identity the size of $\delta$ will be proportional to $\epsilon$.

3.5 Adaptive method

We are now ready to present the full method with adaptively updated correctors. The main idea is to compute the perturbed correctors only for a subset of all elements, and reuse the reference correctors for all other elements. This means we effectively solve the problem in a mixed multiscale space using a bilinear form that is defined as a combination of the two coefficients.

Definition 3.5** (PG-LOD with adaptively updated correctors)**

The proposed method follows five steps:

Provided reference data $(A_{\text{ref}},f_{\text{ref}})$: Compute (for all $T\in{\mathcal{T}}_{H}$) reference correctors ${\mathcal{Q}}^{\text{ref}}_{k,T}\lambda_{x}$ (for all basis functions $\lambda_{x}$), ${\mathcal{R}}^{\text{ref}}_{k,T}f_{\text{ref}}$, and ${\mathcal{Q}}^{\text{ref}}_{k,T}g$, based on the reference coefficient $A_{\text{ref}}$ and reference right hand side $f_{\text{ref}}$. 2. 2.

Provided perturbed data $(A,f)$: Compute (for all $T\in{\mathcal{T}}_{H}$) error indicators $E_{{\mathcal{Q}}V_{H},T}$, $E_{f,T}$, $E_{{\mathcal{R}}f,T}$, and $E_{{\mathcal{Q}}g,T}$ and mark the elements $T$ for which all of the following inequalities hold true,

[TABLE]

Denote the set of marked elements by ${\mathcal{T}}^{\text{ref}}_{H}\subset{\mathcal{T}}_{H}$. 3. 3.

Compute (for all $T\in{\mathcal{T}}_{H}\setminus{\mathcal{T}}^{\text{ref}}_{H}$) the mixed correctors $\tilde{{\mathcal{Q}}}_{k,T}\lambda_{x}$, $\tilde{{\mathcal{Q}}}_{k,T}g$, and $\tilde{{\mathcal{R}}}_{k,T}\tilde{f}$, based on the following definitions of the mixed right hand side and correcctors: $\tilde{f}\in L^{2}(\Omega)$ and

[TABLE]

Note that only element correctors in ${\mathcal{T}}_{H}\setminus{\mathcal{T}}^{\text{ref}}_{H}$ need to be recomputed since the reference correctors from Step 1 will be used for the elements in ${\mathcal{T}}^{\text{ref}}_{H}$. We further let $\tilde{{\mathcal{Q}}}_{k}=\sum_{T\in{\mathcal{T}}_{H}}\tilde{{\mathcal{Q}}}_{k,T}$ and $\tilde{{\mathcal{R}}}_{k}=\sum_{T\in{\mathcal{T}}_{H}}\tilde{{\mathcal{R}}}_{k,T}$. 4. 4.

Assemble the adaptively updated LOD stiffness matrix

[TABLE]

using the mixed unsymmetric bilinear form $\tilde{b}$ defined in terms of a element-wise reference $b_{T}^{\text{ref}}$ and a perturbed $b_{T}$:

[TABLE] 5. 5.

Similarly, define the functional $\tilde{c}$ for the right hand side correctors by

[TABLE]

and solve for $\tilde{u}^{H}_{k}\in V_{H}$ in

[TABLE]

for all $v\in V_{H}$, and compute the solution as

[TABLE]

Remark 3.6** (Individual marking)**

The described algorithm can be enhanced in terms of efficiency by separating ${\mathcal{T}}^{\text{ref}}_{H}$ for each corrector type ${\mathcal{Q}}_{k,T}\lambda_{x}$, ${\mathcal{R}}_{k,T}f$, and ${\mathcal{Q}}_{k,T}g$, only updating the correctors with respect to their corresponding error indicator. For example, if $A=A_{\text{ref}}$ but $f\neq f_{\text{ref}}$, obviously the $f$-independent correctors ${\mathcal{Q}}_{k,T}\lambda_{x}$ and ${\mathcal{Q}}_{k,T}g$ need not to be recomputed in Step 3, since only ${\mathcal{R}}_{k,T}f$ can differ from its reference counterpart ${\mathcal{R}}^{\text{ref}}_{k,T}f_{\text{ref}}$. As stated, however, the algorithm would (unnecessarily) recompute the $f$-independent correctors as well. For readability of this paper, we decided to omit individual marking.

4 Error analysis

This section is devoted to the theoretical justification of the proposed method. We present the main theorem of this work. The theorem justifies local recomputation of the correctors based on the value of the error indicators.

Theorem 4.1** (Error bound for the PG-LOD with adaptively updated correctors)**

If

[TABLE]

then there exist $k_{0}>0$ and $\tau_{0}>0$ such that for all $k>k_{0}$ and $\tau<\tau_{0}$ with $\text{TOL}=\tau k^{-d/2}$ the error bound

[TABLE]

is satisfied. Here, $0<\theta<1$ is independent of $H$, $k$, $\tau$ and TOL.

Before proving the theorem, we require the following lemma.

Lemma 4.2** (Error indicators bound the errors in reference correctors and integrals)**

For all $v\in V_{H}$, the following bounds hold

[TABLE]

Proof** (Lemma 4.2):**

For any $v\in V_{H}$, we define $z:={\mathcal{Q}}_{k,T}v-{\mathcal{Q}}^{\text{ref}}_{k,T}v\in V^{\text{f}}(U_{k}(T))$ and observe

[TABLE]

which yields the first inequality of the first part. We proceed to get the second inequality,

[TABLE]

with $\delta=A^{-1/2}(A-A_{\text{ref}})A_{\text{ref}}^{-1/2}$.

The second result follows analogously with $z^{\prime}:={\mathcal{R}}_{k,T}f-{\mathcal{R}}^{\text{ref}}_{k,T}f_{\text{ref}}\in V^{\text{f}}(U_{k}(T))$ and

[TABLE]

Similar arguments yield the third result of the lemma.

Proof** (Theorem 4.1):**

The full error is ${|\kern-1.29167pt|\kern-1.29167pt|u-\tilde{u}_{k}|\kern-1.29167pt|\kern-1.29167pt|}\leq{|\kern-1.29167pt|\kern-1.29167pt|u-u_{k}|\kern-1.29167pt|\kern-1.29167pt|}+{|\kern-1.29167pt|\kern-1.29167pt|u_{k}-\tilde{u}_{k}|\kern-1.29167pt|\kern-1.29167pt|}$. The first term stems from the localization and is bounded by ${|\kern-1.29167pt|\kern-1.29167pt|u-u_{k}|\kern-1.29167pt|\kern-1.29167pt|}\lesssim k^{d/2}\,\theta^{k}\left(\|f\|_{L^{2}(\Omega)}+{|\kern-1.29167pt|\kern-1.29167pt|g|\kern-1.29167pt|\kern-1.29167pt|}\right)$ according to (3.8).

Before proceeding with the second term, we note that we can bound the error in the global correctors $\tilde{{\mathcal{Q}}}_{k}$ and $\tilde{{\mathcal{R}}}_{k}$ in terms of the patch overlap $k^{d/2}$ and TOL by using Lemma 4.2 and the assumption on TOL stated in this theorem: for all $v\in V_{H}$,

[TABLE]

Analogously, we get ${|\kern-1.29167pt|\kern-1.29167pt|({\mathcal{R}}_{k}-\tilde{{\mathcal{R}}}_{k})f|\kern-1.29167pt|\kern-1.29167pt|}+{|\kern-1.29167pt|\kern-1.29167pt|({\mathcal{Q}}_{k}-\tilde{{\mathcal{Q}}}_{k})g|\kern-1.29167pt|\kern-1.29167pt|}\lesssim k^{d/2}\text{TOL}(\|f\|_{L^{2}(\Omega)}+{|\kern-1.29167pt|\kern-1.29167pt|g|\kern-1.29167pt|\kern-1.29167pt|})$. Using (3.7) and (4.2), we additionally bound

[TABLE]

Next, we proceed with the second term, using (3.6) and (3.13) and the bounds above,

[TABLE]

where we added $\pm\tilde{\mathcal{Q}}_{k}u^{H}_{k}$ in the second estimate and use ${|\kern-1.29167pt|\kern-1.29167pt|u^{H}_{k}|\kern-1.29167pt|\kern-1.29167pt|}\mkern-3.0mu\lesssim\mkern-3.0mu\|f\|_{L^{2}(\Omega)}+{|\kern-1.29167pt|\kern-1.29167pt|g|\kern-1.29167pt|\kern-1.29167pt|}$ for the last step.

It remains to bound the energy norm of $u^{H}_{k}-\tilde{u}^{H}_{k}\in V_{H}$. The first step is to establish a coercivity inequality for $\tilde{b}$ on $V_{H}$. To do this, we define the auxiliary bilinear form $b=\sum_{T\in{\mathcal{T}}_{H}}b_{T}$ (with $b_{T}$ defined in Definition 3.5) and bound the consistency error $[b-\tilde{b}](v,w)$ for $v,w\in V_{H}$,

[TABLE]

since

[TABLE]

where we again use Lemma 4.2. Additionally, for $v\in V_{H}$, we note that

[TABLE]

and further that $a(v-{\mathcal{Q}}_{k}v,v)=b(v,v)$. Using (4.5), (3.7), (4.4), and the triangle inequality, we get for $v\in V_{H}$,

[TABLE]

where $C_{1}$, $C_{2}$, and $C_{3}$ are independent of $H$, $k$ and TOL. From this inequality, we note that there exist $k_{0}$ and $\tau_{0}$ such that, for all $k>k_{0}$ and $\tau<\tau_{0}$ (with $\text{TOL}=\tau k^{-d/2}$), we have the coercivity inequality

[TABLE]

for all $v\in V_{H}$, where $\tilde{\gamma}$ depends on $k_{0}$ and $\tau_{0}$ but not on $H$, $k$, $\tau$ or TOL.

We define $c=\sum_{T\in{\mathcal{T}}_{H}}c_{T}$ and observe that the localized problem (3.5) can be expressed in the same manner as the proposed approximation in (3.12) as

[TABLE]

for all $v\in V_{H}$. Next, we use the coercivity inequality with $v=u^{H}_{k}-\tilde{u}^{H}_{k}=:e\in V_{H}$ and the reformulation of (3.5) together with (3.12) to obtain

[TABLE]

where we used (4.4) for $[\tilde{b}-b](u^{H}_{k},e)$ (and analogously derived results for $[\tilde{c}-c](e)$ and $[\tilde{b}-b](g,e)$). Combining the last bound with (3.8) and (4.3) yields the asserted bound.

Remark 4.3** (Right hand side correction)**

It is possible to ignore ${\mathcal{R}}$ in the method. This leads to an additional error term that is proportional to $H$ in (4.1). However, for localized right hand sides or if high accuracy is needed ${\mathcal{R}}$ should be included, see [14] where this term was first analyzed for LOD.

Remark 4.4** ($V=H^{1}_{0}(\Omega)$)**

The assumption that $V=V_{h}$ (a finite element space) means that the error will be with respect to the finite element approximation $u\in V_{h}$. The same analysis goes through if we instead let $V=H^{1}_{0}(\Omega)$ but the corresponding continuous solution is not computable. It is assumed that there is a fine mesh with mesh size $h$ for which the error in the fine scale finite element solution is small enough. An additional a priori error bound then gives the full error with respect to the continuous solution using the triangle inequality.

5 Implementation

This section discusses implementation details specific to the presented method, with emphasis on the computation of the error indicators, parallel computations and memory consumption for large-scale problems. For implementation details on the LOD corrector problems we refer to [6]. However, we would like to emphasize that the localized computations of correctors on patches makes it possible to avoid any global computations in the (typically very large) space $V_{h}$.

5.1 Computing the error indicators

It is important for the method that the local error indicators $E_{{\mathcal{Q}}V_{H},T}$, $E_{f,T}$, $E_{{\mathcal{R}}f,T}$, and $E_{{\mathcal{Q}}g,T}$ can be computed efficiently. Consider the definition of $E_{{\mathcal{Q}}V_{H},T}$ from (3.9):

[TABLE]

We denote the maximum factor by $\tilde{\mu}_{T,T^{\prime}}$ and note that computing the error corrector for an element $T$ amounts to computing $\|A_{\text{ref}}A^{-1}\|_{L^{\infty}(T)}$, and additionally $\lVert\delta\rVert^{2}_{L^{\infty}(T^{\prime})}=\lVert A^{-1/2}(A_{\text{ref}}-A)A_{\text{ref}}^{-1/2}\rVert^{2}_{L^{\infty}(T^{\prime})}$ and $\tilde{\mu}_{T,T^{\prime}}$ for the coarse elements $T^{\prime}$ in the patch $U_{k}(T)$. The values $\tilde{\mu}_{T,T^{\prime}}$ can be precomputed as $\tilde{\mu}_{T,T^{\prime}}=\max_{l}\mu_{l}$, where $\mu_{l}$ is the solution of the eigenvalue problem

[TABLE]

with

[TABLE]

for $i,j=1,\dots,m-1$ and where $m$ denotes the number of basis functions with support in element $T$. For 2D quadrilateral mesh elements this results in a $3\times 3$ system. Given a perturbed coefficient $A$, we then have to compute only the $\lVert\delta\rVert^{2}_{L^{\infty}(T^{\prime})}$ terms and multiply them with the precomputed $\tilde{\mu}_{T,T^{\prime}}$ before summing, and also computing the $\lVert A_{\text{ref}}A^{-1}\rVert_{L^{\infty}(T)}$ factor. An important consequence of this procedure is that the fine scale reference correctors $\mathcal{Q}_{k,T}^{\text{ref}}\lambda_{j}$ can be discarded and only the coarse scale quantity $\tilde{\mu}_{T,T^{\prime}}$ needs to be stored.

The error indicator $E_{{\mathcal{Q}}g,T}$ for ${{\mathcal{Q}}}^{\text{ref}}_{k,T}g$ can be computed in a similar way as $E_{{\mathcal{Q}}V_{H},T}$, the difference being that no eigenvalue problem over functions in $V_{H}$ restricted to $T$ needs to be solved since $g$ is known a-priori.

The right hand side error indicator $E_{f,T}+E_{{\mathcal{R}}f,T}$ is slightly different, since perturbations of both $f$ and $A$ affect the accuracy of ${{\mathcal{R}}}^{\text{ref}}_{k,T}f_{\text{ref}}$. The second term $E_{{\mathcal{R}}f,T}$ can be computed similarly to $E_{{\mathcal{Q}}V_{H},T}$ and $E_{{\mathcal{Q}}g,T}$. The first term $E_{f,T}$ contains an element local Poincaré-type inequality for functions in the fine space (i.e. the null space of $\mathcal{I}_{H}$),

[TABLE]

We have from the a-priori approximability bound of $\mathcal{I}_{H}$ in (3.1) that this quantity scales with $H$. If the constant $C_{\mathcal{I}_{H}}$ is known it can be used here. In many practical situations it is possible to get a sharp estimate by computational means since $V$ is the finite element space $V_{h}$ and the bound we seek is the maximum eigenvalue to a problem posed on a 1-layer element patch $U(T)$ on the fine scale. For all the experiments in Section 6 we discretize the two dimensional unit square with a uniform grid and get three possible cases: $U(T)$ consists 4, 6 or 9 elements (depending on whether it is in the corner, on the edge or in the interior of the domain). By computing (5.3) with $V=V_{h}$ for a decreasing range of small $h$, it was possible to obtain the estimate $\nu_{T}\leq 0.25H$ for all $T$. This value was used in the experiments in Section 6.

5.2 Algorithm, parallelization and memory consumption

Algorithm 1 shows an example of how the computational steps can be carried out. It is based on the assumption that we have to solve for multiple perturbations of a certain reference coefficient and therefore consists of an initial stage when reference quantities are computed. We make the observations that:

• •

The amount of data stored from the initial stage scales like $k^{d}H^{-d}$ since it consists only of overlapping quantities on the coarse mesh.

• •

The initial stage can be computed in parallel over $T$ thanks to the PG-LOD formulation since $b^{\text{ref}}_{T}(\lambda_{j},\lambda_{i})$, $c^{\text{ref}}_{T}(\lambda_{i})$, and $b^{\text{ref}}_{T}(g,\lambda_{i})$ depend only on the correctors for $T$.

Following the initial stage, a loop over the perturbed coefficients follows. For each perturbed coefficient, there is another loop over the elements $T$. It computes updated correctors for the elements for which the error indicators tell it is needed. Finally, the coarse scale linear system is assembled and solved and the full solution can be computed. We note that:

• •

The iterations in the loop over the perturbed coefficients can be executed in parallel.

• •

The iterations in the loop over the elements $T$ can be executed in parallel.

• •

In order to compute the error indicators, the reference coefficient $A_{\text{ref}}$ needs to either be stored explictly (amount $h^{-d}$), be transmitted patch-wise (amount $k^{d}(h/H)^{-d}$) between computers, or be generated from a low-dimensional representation on demand.

• •

There is a reduction over $T$ to assemble the stiffness matrix and the load vector, but only coarse scale data of amount $k^{d}H^{-d}$ is needed for this reduction. This means that fine scale information does not need to be transmitted between computers during reduction.

• •

The coarse scale solution $\tilde{u}^{H}_{k}$ is readily available after solving the coarse scale linear system. If the full solution $\tilde{u}_{k}$ is requested in some area of the domain, the correctors in that area have to be recomputed or stored from a previous corrector computation.

Remark 5.1** (Periodicity)**

Applications such as composite materials often lead to a periodic structure of the underlying reference coefficient. In this case, correctors can be reused. In the full periodic case this leads to only one corrector problem for full patches and comparably few for the boundary patches. This means that memory consumption as well as complexity decreases significantly.

6 Numerical experiments

In this section we present three experiments where the diffusion coefficient is perturbed by defects with either no domain mapping, local domain mapping or global domain mapping. We performed this experiments with gridlod in Python (see [9]). The entire code for all our experiments is available on GitHub in [10]. Our experiments are performed on a $2$D quadrilateral mesh on $\Omega=[0,1]^{2}$. For the fine and the coarse scale discretization $V_{h}$ and $V_{H}$, we use standard $\mathcal{P}_{1}$ finite element spaces on a fine mesh $\mathcal{T}_{h}$ with $256\times 256$ fine elements and a coarse mesh $\mathcal{T}_{H}$ with $32\times 32$ coarse elements. We use $V_{h}$ for the fine reference solution and for the fine scale discretization of each corrector problem. For the mesh size $H=1/32$ we conclude from the error estimate (3.8) that the localization parameter $k=4\geq|\log(H)|$ suffices.

In all our experiments we use a reference coefficient that is piece-wise constant on every fine mesh element. Furthermore $A_{\text{ref}}$ can be expressed as stated in Section 2.1, i.e. it takes two values $1$ and $\alpha$. For the background $\Omega_{\alpha}$ we choose $\alpha=0.1$. In addition, $\Omega_{1}$ is always a non connected subdomain of $\Omega$ which is conforming with respect to $\mathcal{T}_{h}$. For the sake of convenience we neglect an explicit definition of $A_{\text{ref}}$ and $A$ as we visualize them in the figures. We call a perturbation a (local) defect when a (fully connected) subdomain of $\Omega_{1}$ becomes a subdomain of $\Omega_{\alpha}$. In all our experiments $\Omega_{1}$ is represented by black squares and a defect means that a square gets equalized to the white background (compare Figure 6.1). These defects always occur with a probability of $2\%$. The right hand side $f_{y}$ is always defined by $f_{y}(y)=\chi_{[1/8,7/8]^{2}}(y)$. For the sake of simplicity we use zero Dirichlet boundary conditions $g\equiv 0$.

In the experiments we consider the relative error

[TABLE]

where $u_{k}$ is the best PGLOD solution of (3.6) and $\tilde{u}_{k}$ is the solution of (3.13) for a specific the tolerance TOL in the algorithm of Definition 3.5. For $\text{TOL}=\infty$, we clearly have $0\%$ updates of the correctors whereas $\text{TOL}=0$ corresponds to $100\%$ updates. For $100\%$ updates we then end up with the standard Petrov–Galerkin LOD error that is dependent on our data and discretizations. In order to observe the complete behavior of $\mathcal{E}_{\text{rel}}(u_{k},\tilde{u}_{k})$, we compute $\mathcal{E}_{\text{rel}}(u_{k},\tilde{u}_{k})$ for every possible choice of TOL (and thus for every percentage of updates). The relative best PGLOD error $\mathcal{E}_{\text{rel}}(u,u_{k})$ is always around $10^{-3}$ which means that we are comparing to a sufficiently accurate solution.

6.1 Defects

In the first experiment we let $\psi=\operatorname{Id}$ which means we only consider defects. Figure 6.1 displays the coefficient and its perturbation whereas the error indicators $E_{{\mathcal{Q}}V_{H},T}$ and ${E_{{\mathcal{R}}f,T}}$ are plotted for each $T$ in Figure 6.2. Note that ${E_{f,T}}$ and ${E_{{\mathcal{Q}}g,T}}$ are zero for this example. The coarse mesh is visible in the background. We can clearly see that the error indicator $E_{{\mathcal{Q}}V_{H},T}$ detects the defects in the coefficient correctly. Furthermore $E_{{\mathcal{Q}}V_{H},T}$ is exponentially decaying away from each defect. In a view of ${E_{{\mathcal{R}}f,T}}$, we see that its support coincided with the support of $f_{\text{ref}}$ since for $\psi=\operatorname{Id}$ we have $f=f_{\text{ref}}$ which means that whenever $f$ is zero on $T$, both ${\mathcal{R}}_{k,T}$ and ${\mathcal{R}}^{\text{ref}}_{k,T}$ are zero and thus $E_{{\mathcal{R}}f,T}=0$. We also observe that $E_{{\mathcal{Q}}V_{H},T}$ is significantly greater than $E_{{\mathcal{R}}f,T}$. From $\mathcal{E}_{\text{rel}}(u_{k},\tilde{u}_{k})$ in Figure 6.3 we see a high improvement for few updates of the correctors and a sufficiently fast convergence to the best PGLOD solution. From this experiment we conclude that the method can efficiently be used for local defects.

6.2 Local domain mappings

In the second experiment we choose $\psi$ to be a local distortion in the middle of the domain which can be seen in Figure 6.4. With the help of domain mappings the reference coefficient $A_{\text{ref}}$ is subjected to a simple change in value which means that $\Omega_{1}$ does not change its position. This is visualized in the right picture of Figure 6.4. The domain mapping as well as the defects can be clearly seen in $E_{{\mathcal{Q}}V_{H},T}$, whereas the defects stick out compared to the domain mapping. In contrast to the diffusion coefficient, $f$ is a distortion of $f_{\text{ref}}=f_{y}$. This means that the support of the two is not the same which can be seen in ${E_{f,T}}$ in Figure 6.5 . Furthermore $f$ is not affected by the defects in $A$ which explains that ${E_{f,T}}$ only detects the domain mapping and not the defects. In a view of Figure 6.6, we observe a similar effect as in the first experiment. However, due to the domain mapping, it takes comparatively longer to converge to the optimal PGLOD solution.

6.3 Global domain mapping

In our third experiment we address the situation that $\psi$ has support on the whole domain $\Omega$. The perturbation in the physical domain as well as the change in value of the reference coefficient can be seen in Figure 6.7. The yellow background in $E_{{\mathcal{Q}}V_{H},T}$ in Figure 6.8 visualizes the effect of the domain mapping and the defects are clearly notable. In ${E_{f,T}}$ we once again see the change of support in $f$ (left and right black channel) and the defects are visible in ${E_{{\mathcal{R}}f,T}}$.

In Figure 6.9, we observe that local defects are still resolved efficiently whereas the global map causes a relatively low convergence to the optimal PGLOD solution. Thus, this example shows that global domain mappings are difficult to handle. However, for instance for the case that an accuracy of $10^{-2}$ is accurate enough, we still get a reasonable result. Note that in this particular example no use of domain mappings would surely result in $100\%$ recomputation as the complete coefficient changes.

Remark 6.1** (Choosing TOL)**

In all our experiments, $\mathcal{E}_{\text{rel}}(u_{k},\tilde{u}_{k})$ shows a promising behavior. However, we point out that in practice, we clearly do not know $\mathcal{E}_{\text{rel}}(u_{k},\tilde{u}_{k})$ a priori. Thus, we start with a rather high tolerance TOL and consider an update error $\mathcal{E}_{\text{rel}}(\tilde{u}_{k,\text{old}},\tilde{u}_{k})$, where $\tilde{u}_{k,\text{old}}$ denotes on old approximation with respect to the former tolerance. Whenever the number of updates from one tolerance to another is strictly positive and $\mathcal{E}_{\text{rel}}(\tilde{u}_{k,\text{old}},\tilde{u}_{k})$ does still exhibit a significant gain, we should proceed with a smaller TOL and continue this algorithm until $\mathcal{E}_{\text{rel}}(\tilde{u}_{k,\text{old}},\tilde{u}_{k})$ is small enough.