Numerical homogenization of elliptic PDEs with similar coefficients

TL;DR

This paper introduces a parallelizable Petrov-Galerkin localized orthogonal decomposition algorithm for efficiently solving sequences of elliptic PDEs with similar, rapidly varying coefficients, applicable in time-dependent and stochastic contexts.

Contribution

The paper develops an adaptive PG-LOD method that selectively recomputes local correctors, improving efficiency for sequences of similar elliptic PDEs.

Findings

The method effectively handles sequences with similar coefficients.

Adaptive recomputation enhances computational efficiency.

Application demonstrated on 3D time-dependent Darcy flow.

Abstract

We consider a sequence of elliptic partial differential equations (PDEs) with different but similar rapidly varying coefficients. Such sequences appear, for example, in splitting schemes for time-dependent problems (with one coefficient per time step) and in sample based stochastic integration of outputs from an elliptic PDE (with one coefficient per sample member). We propose a parallelizable algorithm based on Petrov-Galerkin localized orthogonal decomposition (PG-LOD) that adaptively (using computable and theoretically derived error indicators) recomputes the local corrector problems only where it improves accuracy. The method is illustrated in detail by an example of a time-dependent two-pase Darcy flow problem in three dimensions.

Figures

Equations

$$-\mathrm{div}{A}^n\nabla {\bar{u}}^n$$

$${\bar{u}}^n$$

$$\mathbf{n}\cdot A^n\nabla {\bar{u}}^n$$

$$(A\nabla u,\nabla v)=(f,v)-(A\nabla g,\nabla v).$$

$$H^{-1}\Vert v-{\mathcal{I}}_{H}v{\Vert }_{L^2(T)}+\Vert \nabla (v-{\mathcal{I}}_{H}v){\Vert }_{L^2(T)}\le C_{\mathcal{I}}\Vert \nabla v{\Vert }_{L^2(U(T))}.$$

$$U(T)=\bigcup \{T^{\prime }\in {\mathcal{T}}_H:\overline{T}\cap \overline{T^{\prime }}\ne \varnothing \}.$$

$$(A\nabla (u_H+u^{\mathrm{f}}),\nabla v_H)$$

$$(A\nabla u^{\mathrm{f}},\nabla v^{\mathrm{f}})$$

$$(A\nabla \mathcal{Q}v,\nabla v^{\mathrm{f}})$$

$$(A\nabla \mathcal{R}f,\nabla v^{\mathrm{f}})$$

$$(A\nabla u^{\mathrm{ms}},\nabla v_H)=(f,v_H)-(A\nabla g,\nabla v_H)-(A\nabla \mathcal{R}f,\nabla v_H)+(A\nabla \mathcal{Q}g,\nabla v_H).$$

$$U_{k+1}(T)=\bigcup \{T^{\prime }\in {\mathcal{T}}_H:\overline{U_k(T)}\cap \overline{T^{\prime }}\ne \varnothing \}.$$

$$V^{\mathrm{f}}(U_k(T))=\{v\in V^{\mathrm{f}}:v{|}_{\Omega \setminus U_k(T)}=0\},$$

$$(A\nabla {\mathcal{Q}}_{k,T}v,\nabla v^{\mathrm{f}})$$

$$(A\nabla {\mathcal{R}}_{k,T}f,\nabla v^{\mathrm{f}})$$

$$(A\nabla u_k^{\mathrm{ms}},\nabla v_H)$$

$$(A\nabla {\mathcal{R}}_{k}f,\nabla v_H)+(A\nabla {\mathcal{Q}}_{k}g,\nabla v_H),$$

$$({\tilde{A}}_T\nabla {\tilde{\mathcal{Q}}}_{k,T}v,\nabla v^{\mathrm{f}})$$

$$({\tilde{A}}_T\nabla {\tilde{\mathcal{R}}}_{k,T}f,\nabla v^{\mathrm{f}})$$

$$(A\nabla {\tilde{u}}_k^{\mathrm{ms}},\nabla v_H)=(f,v_H)-(A\nabla g,\nabla v_H)-(A\nabla {\tilde{\mathcal{R}}}_{k}f,\nabla v_H)+(A\nabla {\tilde{\mathcal{Q}}}_{k}g,\nabla v_H)$$

$$\tilde{a}(u,v):=\sum _{T\in {\mathcal{T}}_H}{\tilde{a}}_T(u,v):=\sum _{T\in {\mathcal{T}}_H}({\tilde{A}}_T({\chi }_T\nabla -\nabla {\tilde{\mathcal{Q}}}_{k,T}){\mathcal{I}}_{H}u,\nabla v)$$

$$\tilde{L}(v)$$

$$:=\sum _{T\in {\mathcal{T}}_H}(f,v{)}_T-(A\nabla g,\nabla v{)}_T-({\tilde{A}}_T\nabla {\tilde{\mathcal{R}}}_{k,T}f,\nabla v)+({\tilde{A}}_T\nabla {\tilde{\mathcal{Q}}}_{k,T}g,\nabla v).$$

$$\tilde{a}({\hat{u}}_k^{\mathrm{ms}},v_H)=\tilde{L}(v_H).$$

$$|{\mathcal{Q}}_{k,T}v-{\tilde{\mathcal{Q}}}_{k,T}v{|}_A$$

$$|{\mathcal{R}}_{k,T}f-{\tilde{\mathcal{R}}}_{k,T}f{|}_A$$

$$|{\mathcal{Q}}_{k,T}g-{\tilde{\mathcal{Q}}}_{k,T}g{|}_A$$

$$e_{u,T}$$

$$e_{f,T}$$

$$e_{g,T}$$

$$e_u=\underset{T\in {\mathcal{T}}_H}{\max }{e}_{u,T},e_f=\underset{T\in {\mathcal{T}}_H}{\max }{e}_{f,T},\text{ and }{e}_g=\underset{T\in {\mathcal{T}}_H}{\max }{e}_{g,T}.$$

$$|z{|}_{A,U_k(T)}^2$$

$$=(({\tilde{A}}_T-A)\nabla {\tilde{\mathcal{Q}}}_{k,T}v,\nabla z{)}_{U_k(T)}-(({\tilde{A}}_T-A)\nabla v,\nabla z{)}_T$$

$$\le \Vert ({\tilde{A}}_T-A)A^{-1/2}({\chi }_T\nabla v-\nabla {\tilde{\mathcal{Q}}}_{k,T}v){\Vert }_{L^2(U_k(T))}\cdot |z{|}_{A,U_k(T)}.$$

$$\Vert ({\tilde{A}}_T-A)A^{-1/2}({\chi }_T\nabla w-\nabla {\tilde{\mathcal{Q}}}_{k,T}w){\Vert }_{L^2(U_k(T))}$$

$$\le \Vert ({\tilde{A}}_T-A)A^{-1}{\Vert }_{L^{\infty }(T)}+$$

$$\Vert ({\tilde{A}}_T-A)A^{-1/2}{\tilde{A}}_T^{-1/2}{\Vert }_{L^{\infty }(U_k(T))}\Vert A^{-1/2}{\tilde{A}}_T^{1/2}{\Vert }_{L^{\infty }(T)}.$$

$$|z{|}_{A,U_k(T)}^2$$

Full text

1 Introduction

We consider a sequence of elliptic partial differential equations (PDEs) with different, but in some sense similar, rapidly varying coefficients. In some applications, the difference between consecutive coefficients in the sequence is localized, for example for certain Darcy flow applications and in the simulation of random defects in composite materials. This paper studies an opportunity to exploit that the differences are localized to save computational work in the context of the localized orthogonal decomposition method (LOD, [20]).

The accuracy of Galerkin projection onto standard finite element spaces generally suffers from variations in the coefficient that are not resolved by the finite element mesh. The work [3] studies an elliptic equation in 1D with a rapidly varying coefficient and notes that coefficient variations within the element lead to inaccurate solutions for the standard finite element method. Replacing the coefficient with its elementwise harmonic average leads to an accurate method. This result, however, does not easily generalize to higher dimensions. For periodic and semi-periodic coefficients varying on an asymptotically fine scale, a homogenized coefficient can be computed and used for coarse scale computations also in higher dimensions [5]. The early multiscale method [15] is based on homogenization theory and works under assumptions on scale separation and periodicity. Many recent contributions within the field of numerical homogenization can be used without assumptions on periodicity and in higher dimensions, see e.g. [4, 16, 18, 20, 23]. In this work, we consider the LOD technique [20] in the Petrov–Galerkin formulation (PG-LOD) studied in detail in [9].

The fundamental idea of the LOD method is that a low-dimensional function space (multiscale space) with good approximation properties is constructed by computing localized fine-scale correctors to the basis functions of a standard low-dimensional coarse finite element space based on a coarse mesh. Each localized corrector problem is posed only within a patch of a certain radius around its coarse basis function and thus depends only on the diffusion coefficient in that patch. The PG-LOD method has several good properties from a computational perspective. The main advantage is that the PG-LOD corrector problems can be computed completely in parallel with the only communication being a final reduction to form a low-dimensional global stiffness matrix. Further, the fine-scale coefficient only needs to be accessible and stored in memory for one localized corrector problem at a time. Additionally, the method is robust in the sense that both the localized corrector problems and the global low-dimensional problems are typically small enough to be solved with a direct solver.

Once computed, the correctors can be reused for problems with the same or similar diffusion coefficient. We study the case when the diffusion coefficient varies in a sequence of problems. In such situations, there is an opportunity to reuse previously computed localized correctors if the coefficients do not vary too much between consecutive problems. Since the computational cost is proportional to the number of localized corrector problems that have to be recomputed, it is most advantageous if the perturbations of the coefficient are localized. Two practical examples are two-phase flow where the coefficient depends on the saturation of the two fluids, or when the coefficient is a deviation from a base coefficient as in the case with defects in composite materials.

In this work we derive computable error indicators for the error introduced by refraining from recomputing a corrector after a perturbation in the coefficient. The method we propose computes all localized correctors and global stiffness matrix contributions for the first coefficient in the sequence of elliptic PDEs. For the subsequent coefficients, we use the error indicators to adaptively recompute only the correctors that need to be recomputed in order to get a sufficiently accurate solution. The coefficients that have not been recomputed we call lagging coefficients. The method is completely parallelizable over the elements of the coarse mesh. A particularly interesting setting is when only quantities on the coarse mesh are required from the solution, for example upscaled Darcy fluxes in a Darcy flow problem, or the coarse interpolation of the full solution. Any computed fine scale quantities can then be forgotten between the iterations in the sequence and the memory requirement becomes very low.

The paper is divided into five sections: Problem formulation in Section 2, method description in Section 3, error analysis in Section 4, implementation in Section 5, and numerical experiments in Section 6. Both the method description and the error analysis are divided into four steps, with increasing level of approximation in each step: (i) reformulation by variational multiscale method (VMS), (ii) localization by LOD, (iii) approximation of localized correctors by lagging coefficient, and (iv) approximation of global stiffness matrix contribution by lagging coefficient. The main results are the method (12) in Section 3.4, the error bound in Theorem 4 and Algorithm 1.

2 Problem formulation

Let $\Omega$ be a polygonal domain in $\mathbb{R}^{d}$ (with $d=1,2$ or $3$) with the boundary partitioned into disjoint subsets $\Gamma_{D}$ (for Dirichlet boundary conditions) and $\Gamma_{N}$ (for Neumann boundary conditions). Suppose we have a sequence of elliptic equations: for $n=1,2,\ldots$, solve for $\bar{u}^{n}$, such that

[TABLE]

where $f\in L^{2}(\Omega)$, $g\in H^{1/2}(\Gamma_{D})$, ${\bf n}$ is the outward normal of the boundary, and $A^{n}\in L^{\infty}(\Omega)$ is a coefficient varying significantly over small distances. To keep the presentation short, we limit ourselves to the case where $f$ and $g$ are independent of $n$, however, the analysis in this paper can be generalized to $n$-dependent $f$ and $g$. We will refer to the sequence index or rank as time step throughout the paper, although it does not need to correspond to a step from a time-disceratization. For instance, in Section 5 we briefly discuss an application for simulation of weakly random defects in composite materials, where the sequence index corresponds to a Monte Carlo sample member index.

In the remainder of this section and Sections 3–4, we consider a fixed step $n$ and drop this index for all quantities. We call the coefficient $A=A^{n}$ at the current time step the true coefficient. Ideally, only the true coefficient $A$ would be used in the solution at the current time step. However, in order to lower the computational cost, computations from previous time steps will be reused. This means coefficients from previous time steps (lagging coefficients) will enter the analysis through the definition of the localized correctors and in the assembly of the global stiffness matrix. These lagging coefficients will be denoted by $\tilde{A}$. We also want to emphasize that the error indicators derived here are applicable also to situations where the coefficient deviates from a base coefficient, for example within the application of simulations of weakly random defects in composite materials.

We will work with a weak formulation of the above problem. Let $V=\{v\in H^{1}(\Omega)\,:\,v|_{\Gamma_{D}}=0\}$. In case $\Gamma_{D}$ is empty, we instead consider only solutions and test functions in the quotient space $V=H^{1}(\Omega)/\mathbb{R}$. Let $(\cdot,\cdot)$ denote the $L^{2}$-scalar product over $\Omega$, and $(u,v)_{\omega}=\int_{\omega}uv$. Further, we define $\|v\|^{2}_{L^{2}(\omega)}=\int_{\omega}v^{2}$, $\|v\|_{L^{2}}=\|v\|_{L^{2}(\Omega)}$, and the bilinear form $a(u,v)=(A\nabla u,\nabla v)$. We let $\bar{u}=u+g$, where $u\in V$ and $g\in H^{1}(\Omega)$ is an extension of the boundary condition $g$ to the full domain and seek to find $u\in V$, such that for all $v\in V$,

[TABLE]

Assuming there exist constants $0<\alpha$ and $\beta<\infty$, so that $\alpha\leq A\leq\beta$ a.e., $(A\nabla\cdot,\nabla\cdot)$ is bounded and coercive on $V$ and existence of a unique solution is guaranteed by the Lax–Milgram theorem. We further define the energy norm $|\cdot|_{A}=(A\nabla\cdot,\nabla\cdot)^{1/2}$ on $V$, and the semi-norm $|\cdot|_{A,\omega}=(A\nabla\cdot,\nabla\cdot)_{\omega}^{1/2}$.

3 Method description

In this section, we describe the proposed numerical method in a series of steps, each of which introduces another level of approximation for the problem (2) above.

3.1 Variational multiscale method

The first step is to reformulate the problem using the variational multiscale method [16, 17]. This formulation forms the basis for the LOD approximation and makes it possible to reduce the dimensionality of the problem once the corrector problems have been solved.

Let $\mathcal{T}_{H}$ be a regular and quasi-uniform family of conforming subdivisions of $\Omega$ into elements of maximum diameter $H$, and $V_{H}\subset V$ be a family of conforming first order finite element spaces on this mesh, e.g. $\mathcal{P}_{1}$ or $\mathcal{Q}_{1}$ depending on the shape of the elements. The choice of a linear projective quasi-interpolation operator $\mathcal{I}_{H}:V\to V_{H}$ defines the fine space as its kernel $V^{\rm f}=\ker\mathcal{I}_{H}=\{v\in V\,:\,\mathcal{I}_{H}v=0\}$. We assume there exists a constant $C$ independent of $H$ so that for all $v\in V$ and $T\in\mathcal{T}_{H}$, it holds

[TABLE]

Here $U(T)$ is the union of all neighboring elements to $T$, i.e.

[TABLE]

Since we assume $\mathcal{I}_{H}$ is projective (this is not strictly necessary, see e.g. [9]), we have the decomposition $V=V_{H}\oplus V^{\rm f}$ and can decompose the solution $u=u_{H}+u^{\rm f}$ and test function $v=v_{H}+v^{\rm f}$ and test (2) with the two spaces separately:

[TABLE]

We note that $u^{\rm f}$ is linear in $f$ and $u_{H}$, and we define the linear correction operators $\mathcal{Q}:H^{1}(\Omega)\to V^{\rm f}$ and $\mathcal{R}:L^{2}(\Omega)\to V^{\rm f}$, so that $u^{\rm f}=-\mathcal{Q}u_{H}+\mathcal{R}f-\mathcal{Q}g$, i.e., find $\mathcal{Q}v\in V^{\rm f}$ and $\mathcal{R}f\in V^{\rm f}$, such that for all $v^{\rm f}\in V^{\rm f}$,

[TABLE]

These equations have unique solutions, since $(A\nabla\cdot,\nabla\cdot)$ is still bounded and coercive on a subspace $V^{\rm f}\subset V$.

We introduce a new space, the multiscale space, $V^{\rm ms}=V_{H}-\mathcal{Q}V_{H}=\{v_{H}-\mathcal{Q}v_{H}\,:\,v_{H}\in V_{H}\}$, and note that we have the orthogonality relation $V^{\rm ms}\perp_{a}V^{\rm f}$. The solutions $\mathcal{Q}u_{H}$, $\mathcal{R}f$, and $\mathcal{Q}g$ can be plugged into (4a) and we get the following low-dimensional Petrov–Galerkin problem, find $u^{\rm ms}\in V^{\rm ms}$, such that for all $v_{H}\in V_{H}$,

[TABLE]

The full solution is then $u=u^{\rm ms}+\mathcal{R}f-\mathcal{Q}g$.

Remark 1** (Right hand side correction $\mathcal{R}f$).**

It is possible to obtain an approximate solution even if neglecting the right hand side correction term, i.e. letting $\mathcal{R}=0$ above. See for example [14, 20].

3.2 Localized orthogonal decomposition

The second step is to localize the corrector computations by means of localized orthogonal decomposition (LOD). The basic idea is to solve the corrector problems (5) only on localized patches instead of on the full domain to reduce the computational cost.

For the localization, we define element patches for $T\in\mathcal{T}_{H}$, $U_{k}(T)\subset\Omega$, where $0\leq k\in\mathbb{N}$. With trivial case $U_{0}(T)=T$, $U_{k}(T)$ (a $k$-layer element patch around $T$) is defined by the recursive relation

[TABLE]

See Figure 1 for an illustration of element patches.

We further define localized fine spaces

[TABLE]

consisting of fine functions which are zero outside element patches. Throughout the paper, localized quantities are subscripted with the patch size $k$.

Instead of solving (5), we compute the operators $\mathcal{Q}_{k}=\sum_{T\in\mathcal{T}_{H}}\mathcal{Q}_{k,T}$ and $\mathcal{R}_{k}=\sum_{T\in\mathcal{T}_{H}}\mathcal{R}_{k,T}$, with $\mathcal{Q}_{k,T}$ and $\mathcal{R}_{k,T}$ defined by

[TABLE]

for all $v^{\rm f}\in V^{\rm f}(U_{k}(T))$ and all $T\in\mathcal{T}_{H}$.

We define the localized multiscale space $V_{k}^{\rm ms}=V_{H}-\mathcal{Q}_{k}V_{H}$. Our localized multiscale problem reads find ${u}_{k}^{\rm ms}\in V_{k}^{\rm ms}$, such that for all $v_{H}\in V_{H}$,

[TABLE]

and the full solution for the second approximation is $u_{k}={u}_{k}^{\rm ms}+\mathcal{R}_{k}f-\mathcal{Q}_{k}g$.

3.3 Lagging multiscale space

In the third approximation we compute the localized element correctors using a lagging coefficient ${\tilde{A}}_{T}$ rather than the true $A$. This makes it possible to reuse correctors that have been computed at earlier time steps, so that localized correctors only for a small number of elements $T$ need to be recomputed.

We define the lagging localized corrector operators ${\tilde{\mathcal{Q}}}_{k}=\sum_{T}{\tilde{\mathcal{Q}}}_{k,T}$ and ${\tilde{\mathcal{R}}}_{k}=\sum_{T}{\tilde{\mathcal{R}}}_{k,T}$. The element corrector operators ${\tilde{\mathcal{Q}}}_{k,T}v,{\tilde{\mathcal{R}}}_{k,T}f\in V^{\rm f}(U_{k}(T))$ are defined such that for all $v^{\rm f}\in V^{\rm f}(U_{k}(T))$,

[TABLE]

Note that lagging coefficients ${\tilde{A}}_{T}$ are not necessarily the same for all $T$.

Example 1** (Relation between lagging coefficient and time steps).**

As an example, for the current time step $A=A^{n}$, for element $T^{\prime}$ the coefficient can be one time step old, i.e. ${\tilde{A}}_{T^{\prime}}=A^{n-1}$ and for $T^{\prime\prime}$ three time steps old, i.e. ${\tilde{A}}_{T^{\prime\prime}}=A^{n-3}$. That is, different lagging localized element correctors may be defined in terms of coefficients from different time steps in history.

In analogy with previous multiscale spaces, we define a lagging multiscale space ${\tilde{V}}^{\rm ms}_{k}=V_{H}-\mathcal{\tilde{Q}}_{k}V_{H}$ and the problem is then to find ${\tilde{u}}^{\rm ms}_{k}\in{\tilde{V}}^{\rm ms}_{k}$, such that for all $v_{H}\in V_{H}$,

[TABLE]

and the full solution for the third approximation is ${\tilde{u}}_{k}={\tilde{u}}^{\rm ms}_{k}+\mathcal{\tilde{R}}_{k}f-\mathcal{\tilde{Q}}_{k}g$.

3.4 Lagging global stiffness matrix contribution

The fourth approximation involves not only using a lagging multiscale space, but also a lagging coefficient in the assembly of the global stiffness matrix and right hand side. The rationale behind this is that computing the integrals in the stiffness matrix and right hand side for (9) requires that all precomputed element correctors are stored. To circumvent this, we propose the following approximation. First we define a lagging bilinear form $\tilde{a}$ (and its elementwise contributor $\tilde{a}_{T}$), based on the same lagging coefficients ${\tilde{A}}_{T}$ as was used for the multiscale space in the previous section,

[TABLE]

where $\chi_{T}$ is the indicator function for subset $T\subset\Omega$. We also define a lagging linear functional $\tilde{L}$ (and its elementwise contributor $\tilde{L}_{T}$),

[TABLE]

Then the problem is posed as to find ${\hat{u}}^{{\rm ms}}_{k}\in{\tilde{V}}^{\rm ms}_{k}$, such that for all $v_{H}\in V_{H}$,

[TABLE]

The full solution for the final approximation step is then ${\hat{u}}_{k}={\hat{u}}^{{\rm ms}}_{k}+\mathcal{\tilde{R}}_{k}f-\mathcal{\tilde{Q}}_{k}g$.

We note that (12) coincides with (9) when ${\tilde{A}}_{T}=A$ for all $T$. Also, we note that the coefficients for the linear system can be computed immediately after ${\mathcal{\tilde{Q}}}_{k,T}$ and ${\mathcal{\tilde{R}}}_{k,T}$ have been computed. This means no correctors need to exist simultaneously.

This method is independent of the true coefficient $A$, if not ${\tilde{A}}_{T}=A$ for any $T$. In order to construct a numerical method with control of the error from this approximation, we use error indicators on the element correctors to determine whether they need to be recomputed or not. Next, we define three computable error indicators, $e_{u}$, $e_{f}$ and $e_{g}$, for the error introduced by using a lagging coefficient.

3.5 Error indicators

As can be seen in later Sections 4.3 and 4.4, the differences $|\mathcal{Q}_{k,T}v-\mathcal{\tilde{Q}}_{k,T}v|_{A}$ for $v\in V_{H}$, $|\mathcal{R}_{k,T}f-\mathcal{\tilde{R}}_{k,T}f|_{A}$, and $|\mathcal{Q}_{k,T}g-\mathcal{\tilde{Q}}_{k,T}g|_{A}$ constitute the sources to the error in the approximation from using lagging coefficients. In this section, we define three elementwise error indicators ($e_{u,T}$, $e_{f,T}$, and $e_{g,T}$) and relate them to the above differences in Lemma 2.

Lemma 2** (Error indicators: definitions and bounds).**

The following bounds hold,

[TABLE]

where

[TABLE]

We additionally define

[TABLE]

Proof.

For any $v\in V_{H}$, let $z=\mathcal{Q}_{k,T}v-\mathcal{\tilde{Q}}_{k,T}v$, then using (7) and (8), we get

[TABLE]

Then, clearly $e_{u,T}$ (if it exists) constitute the asserted bound. The following inequality gives a bound for the norm being maximized in the definition of $e_{u,T}$ (assuming that $|w|_{A,T}=1$),

[TABLE]

The maximum is thus attained and exists by the extreme value theorem.

Similarly, for $z=\mathcal{R}_{k,T}f-\mathcal{\tilde{R}}_{k,T}f$, we have

[TABLE]

which motivates the definition of $e_{f,T}$ and the asserted bound. The result for $e_{g,T}$ holds analogously. ∎

Regarding the computation of these error indicators, both $e_{f,T}$ and $e_{g,T}$ are straight-forward to compute, being a ratio of two computable norms. The error indicator $e_{u,T}$ is also easy to compute. It is the square root of a Rayleigh quotient for a generalized eigenvalue problem (where the restriction $|w|_{A,T}=1$ removes the singularity of the denominator matrix):

[TABLE]

with the matrices

[TABLE]

for all $i,j=1,\ldots,m-1$ where $m$ is the number of basis functions in $T$ (i.e. one of them removed). The squared maximum $e_{u,T}^{2}$ corresponds to the maximum eigenvalue $\max_{\ell}\mu_{\ell}$. We emphasize that the matrices $\bf B$ and $\bf C$ are very small: the same size as the number of degrees of freedom in the coarse element $T$ (minus one for removing the constant), e.g., $2\times 2$ for 2D simplicial meshes or $7\times 7$ for 3D hexahedral meshes.

3.5.1 Coarse error indicators

In order to compute the error indicators $e_{u,T}$, $e_{f,T}$, and $e_{g,T}$ we need access to the true coefficient $A$ and lagging correctors $\mathcal{\tilde{Q}}_{k,T}\phi_{i}$, $\mathcal{\tilde{R}}_{k,T}f$, and $\mathcal{\tilde{Q}}_{k,T}g$ at the same time. This implies all lagging correctors need to be saved in order to compute the error indicators. Since the correctors in practice are defined on patches of a fine mesh, and the patch overlap can be substantial, the memory requirements for saving them might be large. In this section, we construct an additional bound that makes it possible to discard the lagging correctors after they have been computed.

We construct the following bound starting from the definition of $e_{u,T}$ in Lemma 2,

[TABLE]

where $\delta_{T}=({\tilde{A}}_{T}-A)A^{-1/2}{\tilde{A}}_{T}^{-1/2}$, and we used that

[TABLE]

in the last inequality. We further define $E_{u}=\max_{T\in\mathcal{T}_{H}}E_{u,T}$.

The maximum in (13) corresponds to a maximum eigenvalue of a low-dimensional generalized eigenvalue problem, as was the case for $e_{u,T}$ in Section 3.5. More specifically, it is the square root of the maximum eigenvalue

[TABLE]

of ${\bf B}{\bf x}_{\ell}=\mu_{\ell}{\bf C}{\bf x}_{\ell}$ with the matrices

[TABLE]

for $i,j=1,\ldots,m-1$, where $m$ is the number of basis functions in $T$. We note that the quantity ${\tilde{\mu}}_{T,T^{\prime}}$ can be computed directly after corrector $\mathcal{\tilde{Q}}_{k,T}$ has been computed for the basis functions in element $T$. Now, $\mathcal{\tilde{Q}}_{k,T}$ does not need to be saved for computing $E_{u,T}$ later, and it can be discarded. In particular, the memory required for storing ${\tilde{\mu}}_{T,T^{\prime}}$ (which, however, is needed to compute $E_{u,T}$) scales like $\mathcal{O}(k^{d}H^{-d})$.

Still, ${\tilde{A}}_{T}$ needs to be available to compute $\|A^{-1/2}{\tilde{A}_{T}}^{1/2}\|^{2}_{L^{\infty}(T)}$ and $\delta_{T}$. This might not be a problem in applications where there is a low-dimensional description of the coefficient, for example if the coefficient is defined by a set of geometric shapes which can be described by location, size, shape and so on. In the section for numerical experiments, we will study an example of upscaled two-phase Darcy flow, where we illustrate a way to avoid saving $\tilde{A}_{T}$.

The error indicator $E_{u}$ can replace $e_{u}$ in all results and algorithms in this work. Similar coarse error indicators can be derived for $e_{f}$ and $e_{g}$.

4 Error analysis

In this section we study the approximation error of the three approximations $u_{k}$, ${\tilde{u}}_{k}$ and ${\hat{u}}_{k}$, and the inf-sup stability for the systems yielding the solutions $u^{\rm ms}$, ${u}_{k}^{\rm ms}$, ${\tilde{u}}^{\rm ms}_{k}$ and ${\hat{u}}^{{\rm ms}}_{k}$. Finally, in Theorem 4 in Section 4.4, we present a bound on the error $u-{\hat{u}}_{k}$ of the full approximation.

We use $C$ to denote a constant that is independent of the regularity of $u$, patch size $k$ and coarse mesh size $H$. It can, however, depend on the contrast $\frac{\beta}{\alpha}$. The value of the constant is not tracked between steps in inequalities. By the notation $a\lesssim b$, we mean $a\leq Cb$.

4.1 Variational multiscale method

Since the varitional multiscale formulation (6) is only a reformulation of the original problem, without any approximations, there is no error. However, the well-posedness of the formulation is still of interest.

4.1.1 Stability

Uniqueness of a solution to (6) is guaranteed by an inf-sup condition for $a$ on $V^{\rm ms}$ and $V_{H}$,

[TABLE]

The existence inf-sup condition holds analogously. We let $\gamma$ denote the inf-sup stability constant and note that it is depends on the contrast for general $\mathcal{I}_{H}$. See [12, 24] for corrector localization results independent of the contrast.

4.2 Localized orthogonal decomposition

For the error analysis of LOD we recite previous exponential decay results (first presented in [20]) of the localized corrector operators by means of the following lemmas. For example, the proof in [14, Lemma 3.6] is almost directly applicable here.

Lemma 3** (Localization error).**

Let $k>0$ be a fixed integer and let $p_{T}\in V^{\rm f}$ be the solution of

[TABLE]

for all $v^{\rm f}\in V^{\rm f}$, where $F_{T}\in V^{*}$ such that $F_{T}(v^{\rm f})=0$ for all $v^{\rm f}\in V^{\rm f}(\Omega\setminus T)$. Furthermore, we let $p_{k,T}\in V^{\rm f}(U_{k}(T))$ be the solution of

[TABLE]

for all $v^{\rm f}\in V^{\rm f}(U_{k}(T))$. Then there exists a constant $0<\theta<1$ that depends on the contrast but not on $H$ or the variations of $A$, such that

[TABLE]

This lemma can be applied for the localization error of both $\mathcal{Q}-\mathcal{Q}_{k}$ and $\mathcal{R}-\mathcal{R}_{k}$. In analogy with the definition of $\mathcal{Q}_{k}$ as a sum of $\mathcal{Q}_{k,T}$, we can define $\mathcal{Q}=\sum_{T\in\mathcal{T}_{H}}\mathcal{Q}_{T}$ with $\mathcal{Q}_{T}=\mathcal{Q}_{\infty,T}$. Then for any $v\in V$, we can identify $\mathcal{Q}_{T}v$ with $p_{T}$ and $\mathcal{Q}_{k,T}v$ with $p_{k,T}$ in the lemma above (and similarly for $\mathcal{R}$).

4.2.1 Stability

Using Lemma 3, we get the following result for $\mathcal{Q}v-\mathcal{Q}_{k}v$, with $v\in H^{1}(\Omega)$,

[TABLE]

If in addtion $v\in V_{H}$, we can use the stability of $\mathcal{I}_{H}$ and continue to get

[TABLE]

Using the result above, we can derive an inf-sup constant for $a$ and the pair of spaces $V_{k}^{\rm ms}$ and $V_{H}$,

[TABLE]

For sufficiently large $k$, there is a uniform bound $\gamma_{0}\leq\gamma_{k}$. See [9] for more details on stability of this approximation.

4.2.2 Error

For arbitrary $u_{I}\in V_{k}^{\rm ms}$, using the equations (6) and (7), we have for all $v\in V_{H}$,

[TABLE]

The inf-sup condition for uniqueness above yields the following approximation result, for arbitrary $u_{I}\in V_{k}^{\rm ms}$,

[TABLE]

In analogy with (15) we get the following result for $\mathcal{R}f-\mathcal{R}_{k}f$,

[TABLE]

Recall that $u=\mathcal{I}_{H}u^{\rm ms}-\mathcal{Q}\mathcal{I}_{H}u^{\rm ms}+\mathcal{R}f-\mathcal{Q}g$ and $u_{k}=\mathcal{I}_{H}{u}_{k}^{\rm ms}-\mathcal{Q}_{k}\mathcal{I}_{H}{u}_{k}^{\rm ms}+\mathcal{R}_{k}f-\mathcal{Q}_{k}g$. Now, if we choose $u_{I}=\mathcal{I}_{H}u^{\rm ms}-\mathcal{Q}_{k}\mathcal{I}_{H}u^{\rm ms}\in V_{k}^{\rm ms}$, then $u^{\rm ms}-u_{I}=-(\mathcal{Q}-\mathcal{Q}_{k})\mathcal{I}_{H}u^{\rm ms}$ and using the approximation result we get

[TABLE]

Then, using $\mathcal{I}_{H}u^{\rm ms}=\mathcal{I}_{H}u$, interpolation stability (3) and stability of the continuous problem, we have for the full error

[TABLE]

This result was first shown in [20] and is noteworthy, since the error of the approximation decays exponentially with increasing $k$, independently of the regularity of the solution $u$.

4.3 Lagging multiscale space

For this step, we use a lagging multiscale space ${\tilde{V}}^{\rm ms}_{k}$ and need to establish an inf-sup stability constant for $a$ with respect to ${\tilde{V}}^{\rm ms}_{k}$ and $V_{H}$. We will use the results from Lemma 2 both for deriving stability and the approximation error. The following full corrector error can be derived using Lemma 2,

[TABLE]

The bounds $|\mathcal{R}_{k}f-\mathcal{\tilde{R}}_{k}f|^{2}_{A}\lesssim k^{d}e^{2}_{f}\|f\|^{2}_{L^{2}}$ and $|\mathcal{Q}_{k}g-\mathcal{\tilde{Q}}_{k}g|^{2}_{A}\lesssim k^{d}e^{2}_{g}|g|^{2}_{A}$ hold similarly.

We note that if ${\tilde{A}}_{T}=A$, then $e_{u,T}=e_{f,T}=e_{g,T}=0$. Obviously, updating a lagging coefficient for an element corrector leads to no error for this element corrector.

4.3.1 Stability

We can now derive an inf-sup constant for $a$ on ${\tilde{V}}^{\rm ms}_{k}$ and $V_{H}$, using similar techniques as in (16),

[TABLE]

We note that $k$ enters the constant, but that it can be compensated by a small $e_{u}$. Since $e_{u,T}$ is computable, a rule to recompute all element correctors $T$ with $e_{u,T}\geq{\rm TOL}(k)$ for some small enough ${\rm TOL}(k)=\mathcal{O}(k^{-d/2})$, will (after recomputation) make ${\tilde{A}}_{T}=A$ and $e_{u,T}=0$. This makes $e_{u}<{\rm TOL}(k)$. Following this adaptive rule makes it possible to find a lower bound $\tilde{\gamma}_{0}\leq\tilde{\gamma}_{k}$ for sufficiently large $k$ and sufficiently small ${\rm TOL}$.

4.3.2 Error

Again, we get an approximation result from the inf-sup stability. In complete analogy with (LABEL:eq:through_interpolation) and (18), we get

[TABLE]

4.4 Lagging global stiffness matrix contribution

In the fourth approximation (12), the coefficients for the integration of the global stiffness matrix and (parts of) the right hand side are also lagging.

4.4.1 Stability

We derive an inf-sup constant for $\tilde{a}$ (see (10)) with respect to ${\tilde{V}}^{\rm ms}_{k}$ and $V_{H}$,

[TABLE]

Again, $k$ enters, but can be compensated by a small $e_{u}$ according to the discussion in Section 4.3.1. Thus, there is a bound $\hat{\gamma}_{0}\leq\hat{\gamma}_{k}$ for all sufficiently large $k$.

4.4.2 Error

To study the error $|{\tilde{u}}_{k}-{\hat{u}}_{k}|_{A}$, we first note that $|{\tilde{u}}_{k}-{\hat{u}}_{k}|_{A}=|{\tilde{u}}^{\rm ms}_{k}-{\hat{u}}^{{\rm ms}}_{k}|_{A}$, since the right hand side and boundary condition corrections are the same in both cases. We form the following difference from (9) and (12),

[TABLE]

Add and subtract $a({\hat{u}}^{{\rm ms}}_{k},v_{H})$ and use Lemma 2 to get

[TABLE]

Inf-sup stability for $a$ and $\tilde{a}$ finally gives, for any $v_{H}\in V_{H}$,

[TABLE]

We conclude this section by presenting the main theoretical result of this paper. It gives a bound of the full error of ${\hat{u}}_{k}$ (in energy norm) in terms of the patch size $k$ and the error indicators $e_{u}$, $e_{f}$, and $e_{g}$ defined in Lemma 2. This theorem forms the basis for the implementation of a method that updates the multiscale space adaptively while iterating through the sequence of coefficients.

Theorem 4** (Error bound for multiscale method with lagging coefficient).**

Assume $k$ is sufficiently large, so that $\gamma_{k}\geq\gamma_{0}$ holds. Let ${\rm TOL}=ck^{-d/2}$ and (by recomputation of element correctors) $\max(e_{u},e_{g},e_{f})\leq{\rm TOL}$. Choose $c$ sufficiently small so that $\tilde{\gamma}_{k}\geq\tilde{\gamma}_{0}$, and $\hat{\gamma}_{k}\geq\hat{\gamma}_{0}$. Further, let $u$ solve (2) and ${\hat{u}}^{{\rm ms}}_{k}$ solve (12). Let ${\hat{u}}_{k}={\hat{u}}^{{\rm ms}}_{k}+\mathcal{\tilde{R}}_{k}f-\mathcal{\tilde{Q}}_{k}g$. Then

[TABLE]

where the hidden constant depends on the contrast but is independent of mesh size $H$, patch size $k$ and regularity of the solution $u$.

Proof.

The estimate of the full error $|u-{\hat{u}}_{k}|_{A}$ is obtained by combining (18), (LABEL:eq:error_falsespace), and (22), and using the triangle inequality,

[TABLE]

and finally using the assumed bounds of $e_{u}$, $e_{f}$ and $e_{g}$. ∎

Remark 5** (Selecting parameters $H$, $k$ and ${\rm TOL}$).**

The coarse mesh size parameter $H$ is typically chosen based the desired accuracy of the computation. The localization parameter $k$ is chosen to be proportional to $|\log(H)|$ guaranteeing a perturbation of the approximation of the order $H|\log(H)|^{d/2}$. Finally, ${\rm TOL}$ is chosen proportional to $H$. The resulting error bound in energy norm then reads $\lesssim|\log(H)|^{d/2}H(\|f\|_{L^{2}}+|g|_{A})$.

5 Implementation

In this section, we present an algorithm for computing approximate solutions to a sequence of problems as described by (1). In a practical implementation we can not let $V$ be an infinite dimensional space. We will assume that there is a finite element space $V_{h}$ based on a mesh that resolves the coefficient, which if used to solve (2), yields an approximate solution $u_{h}$ with satisfactory small error $|u-u_{h}|_{A}$. The analysis in the previous sections holds also if replacing $V$ with $V_{h}$, however, the error estimates will then of course be bounding $|u_{h}-{\hat{u}}_{k}|_{A}$ instead of $|u-{\hat{u}}_{k}|_{A}$. In the end of the section we discuss the memory requirements of the algorithm.

The key idea is that, as time $n$ progresses, we do not update the full multiscale space, but only the parts where it is necessary for a sufficiently small error. If $A^{n}$ only changes slightly between two consecutive $n$, it is possible that many of the element correctors ($\mathcal{\tilde{Q}}_{k,T}v_{H}$, $\mathcal{\tilde{R}}_{k,T}f$ and $\mathcal{\tilde{Q}}_{k,T}g$) based on lagging coefficients do not need to be recomputed. We use the error indicators $e_{u}$, $e_{f}$ and $e_{g}$ to determine for which elements to recompute correctors. This results in an algorithm that is completely parallelizable over $T$, except for the solution of the low-dimensional (posed in $V_{H}$) global system. Even the assembly of the global stiffness matrix ${\bf K}=(K_{ij})_{ij}$ and right hand side ${\bf b}=(b_{i})_{i}$ can be done in parallel, as it becomes a reduction over $T$.

The algorithm is presented in Algorithm 1. We denote by $\phi_{i}\in V_{H}$, $i=1,2,\ldots$ the finite element basis functions spanning $V_{H}$.

Note that the if-statement in this algorithm together with properly chosen $k$ and ${\rm TOL}$, ensures that the conditions for Theorem 4 are fulfilled. The numerical experiment in Section 6.1 investigates the relations between the error and ${\rm TOL}$ and the fraction of recomputed element correctors.

The memory required to perform the main algorithm grows with $k$ in the following manner. Suppose $\mathcal{O}(h^{-d})$ is the number of elements in the fine discretizations, as is the case for quasi-uniform meshes. To compute $e_{u}$, $e_{f}$ and $e_{g}$, we need to keep $\tilde{A}_{T}$, $\mathcal{\tilde{Q}}_{k,T}\phi_{i}$, $\mathcal{\tilde{R}}_{k,T}f$, and $\mathcal{\tilde{Q}}_{k,T}g$ between the iterations in Algorithm 1 (see Lemma 2). Since the patches $U_{k}(T)$ overlap by $\mathcal{O}(k^{d})$ coarse elements, the amount of memory required between the iterations scales like $\mathcal{O}(k^{d}h^{-d})$. In high dimensions for very fine meshes, the amount of memory needed for storage can become a limitation. Depending on the application, it is possible to reduce the memory requirements. Below we give two examples of such applications.

Example 2** (Defects in composite materials).**

For simulations on weakly random materials [2], we consider the coefficient of a reference material and a material with random defects shown in Figure 2. If each ball in this material has a certain low probability to be missing (a localized point defect), the proposed method can be used to solve the model problem with the defect material on the right (true coefficient) using correctors precomputed on the reference material on the left (lagging coefficient). In sample based methods for stochastic integration (e.g. Monte Carlo), the proposed method for determining what correctors to recompute can reduce the computational cost for the full simulation.

The lagging coefficient $\tilde{A}_{T}$ in this example is the single reference coefficient and thus the same for all $n$ and all $T$. Because of this, no additional memory is required to store lagging coefficients in this case. If we additionally use the (less efficient) coarse error indicators $E_{u,T}$, $E_{f,T}$, and $E_{g,T}$ presented in Section 3.5.1, our memory requirement scales with $\mathcal{O}(k^{d}H^{-d}+h^{-d})$ between the iterations in the algorithm.

Example 3** (Two-phase Darcy flow).**

In a discretization of a two-phase Darcy flow system of equations (pressure and saturation equation) for an injection scenario, the permeability coefficient $A=A(s^{n})$ varies over time $n$ indirectly through the dependence on the saturation $s^{n}$. Typically, the change in saturation between time steps is localized to the front of the plume of the injected fluid. Thus, most corrector problems can be expected to be reused between iterations. An approach to reducing the memory requirements for the solution of this problem is revisited in detail in Section 6.2.

6 Numerical experiments

In all numerical experiments, we use $\mathcal{Q}_{1}$ Lagrange finite elements in 2D or 3D on rectangular or rectangular cuboid elements. The degrees of freedom are the values of the polynomial in the corners of the element.

We define the interpolation operator $\mathcal{I}_{H}$ to be used throughout the experiments. Let $P_{1}$ denote the polynomials of no partial degree greater than $1$, i.e. $\partial^{2}p/\partial x^{2}=0$ for all independent variables $x$ if $p\in P_{1}$. We define the broken finite element space,

[TABLE]

We denote by $\Pi_{H}$ the $L^{2}$-projection onto $S_{H,\rm{b}}$ and by $E_{H}:S_{H,\rm{b}}\to V_{H}$ the boundary condition conforming node averaging operator (Oswald interpolation operator), for all nodes in $\mathcal{T}_{H}$,

[TABLE]

where $T_{x}=\{T\in\mathcal{T}_{H}\,:\,x\in\overline{T}\}$ and $\operatorname{card}$ is the cardinality. Then we define $\mathcal{I}_{H}=E_{H}\circ\Pi_{H}$. This operator satisfies (3), see e.g. [8].

6.1 Experiments studying the effects of $k$ and ${\rm TOL}$

We let $\Omega=[0,1]^{2}$, $\Gamma_{D}=\{x\in\partial\Omega:x_{1}=0\text{ or }x_{1}=1\}$, $\Gamma_{N}=\partial\Omega\setminus\Gamma_{D}$, $f=0$, $g=1-x_{1}$, and $A_{\rm b}$ as shown in Figure 4(a). $A$ was constructed by taking a uniform grid with $512\times 512$ cells, and assigning each grid cell a value $10^{c}$, where $c$ was drawn from a uniform distribution between $[-2,0]$, for each cell independently. Then, the values in cells whose midpoint $x_{\rm m}=(x_{{\rm m},1},x_{{\rm m},2})$ satisfied $15/32\leq x_{{\rm m},1}\leq 1/2$ were set to $10^{-2}$. Finally, the values in cells whose midpoint $x_{\rm m}$ satisfied $1/4\leq x_{{\rm m},2}\leq 5/16$ were set to $1$.

The space $V$ is discretized on a $\mathcal{Q}_{1}$ finite element space on a uniform grid of size $512\times 512$, see Figure 3. The $V_{H}$ is chosen as a $\mathcal{Q}_{1}$ finite element space on the coarse mesh shown in the same figure.

6.1.1 Error decay with $k$

First, we let $A=A_{\rm b}$ and solve for $u_{k}$ for $k=1,2,3,$ and $4$. We solve for $u$ on the fine mesh and use that as reference solution. The exponential convergence in terms of $k$ can be observed in Figure 5.

6.1.2 Error decay with ${\rm TOL}$

Now we fix $k=3$ and define a sequence of coefficients $A^{n}$ for $n=0,\ldots,127$,

[TABLE]

This describes a perturbation of a factor up to $3$ over the full domain, sweeping from the left to the right. We emphasize that the difference $A_{n+1}-A_{n}$ is nonzero everywhere, which means a strategy to determine which correctors to recompute is necessary. We use Algorithm 1 to compute the approximate solution ${\hat{u}}_{k}$ for every time step $n$. A reference solution $u$ is also computed. We do this for four values of ${\rm TOL}=0.5,0.1,0.05$, and $0.01$.

The (relative) error in energy norm versus the time step $n$ is plotted to the left in Figure 6. The right plot in the same figure shows the fraction of all element correctors $\mathcal{\tilde{Q}}_{k,T}$ that were recomputed in each time step. We note that the error decreases with decreasing ${\rm TOL}$ as expected and that the fraction of recomputed element correctors increase with decreasing ${\rm TOL}$. Without an adaptive strategy, all element correctors would have to be recomputed in every time step. See Figure 7 for two maps over the recomputed element correctors in time step $n=31$ for two different values of ${\rm TOL}$.

6.2 Low-memory Darcy flow upscaling algorithm

In order to continue with two additional numerical experiments (in Section 6.3), we describe an algorithm for pressure solution upscaling for Darcy flows that reduces the space complexity to $\mathcal{O}(k^{d}H^{-d}+h^{-d})$ (from $\mathcal{O}(k^{d}h^{-d})$). This is done by solving the saturation equation on the coarse mesh and the pressure equation with saturation dependent diffusion coefficient on a fine mesh using the proposed adaptive multiscale method. This is possible in a situation where the diffusion coefficient cannot be averaged on the coarse mesh, but the saturation solution can. The low space complexity and the possibility to parallelize the corrector computations enable the solution of large-scale problems of this kind.

6.2.1 A Two-phase Darcy flow model problem

We consider the immiscible non-capillary two-phase Darcy flow problem using the fractional flow formulation [13, 21]. This leads to a system of a coupled pressure and saturation equation

[TABLE]

where space-time functions: $u$, $s$ and $f$ are pressure, saturation for the wetting phase, and sources/sinks, respectively; space function $K$ is intrinsic permeability; and nonlinear scalar functions: $\lambda$ and $\lambda_{w}$ are total mobility and wetting phase mobility, respectively. A common technique used for solving this system is sequential splitting, where the pressure equation and saturation equation are solved separately within a time step $n$. This means that, as we iterate in time, we need to solve a sequence of pressure equations with coefficient $A^{n}(x)=\lambda(s(x,t_{n-1}))K(x)$. Since the wetting saturation $s$ changes only significantly along the plume front between time steps, we are in the setting where consecutive differences in the coefficient is localized.

The permeability $K$ varies on a fine scale, requring these variations to be resolved by a fine mesh with mesh size $h$ in order to obtain an accurate pressure solution. We consider the case when the saturation equation needs only be solved on a coarser mesh with mesh size $H>h$ to obtain a sufficiently accurate saturation solution. We let the fine mesh $\mathcal{T}_{h}$ be a refinement of the coarse mesh $\mathcal{T}_{H}$. The pressure and saturation equations are solved sequentially: Given initial data for the saturation, the pressure equation is solved. An approximation of the coarse element face flux is computed and used to solve for the next saturation using a zeroth order upwind discontinuous Galerkin method with explicit Euler forward time-stepping.

We use the same discretization scheme as in [22]. Let $P_{0}(\mathcal{T}_{H})$ be the space of piecewise constants on the elements of $\mathcal{T}_{H}$. Let $\mathcal{F}_{H}$ denote set of faces of $\mathcal{T}_{H}$. Each face $F\in\mathcal{F}_{H}$ has a normal direction ${\bf n}_{F}$ (outward pointing for boundary faces). We define the jump operator over face $F$ as $[\![v]\!]=(v|_{T_{1}})|_{F}{\bf n}_{1}\cdot{\bf n}_{F}+(v|_{T_{2}})|_{F}{\bf n}_{2}\cdot{\bf n}_{F}$, where ${\bf n}_{1}$ and ${\bf n}_{2}$ are the outward pointing face normals of the two elements $T_{1}$ and $T_{2}$ adjacent to $F$. Let $\langle\cdot,\cdot\rangle_{\omega}$ denote the $L^{2}$-scalar product when $\omega$ is $d-1$-dimensional. We let the flow be completely driven by boundary conditions, i.e. $f=0$. We use the following discretization for the saturation equation. Given $s_{H}^{n-1}\in P_{0}(\mathcal{T}_{H})$, ${\sigma}^{n}\in L^{1}(\mathcal{F}_{H})$, find $s_{H}^{n}\in P_{0}(\mathcal{T}_{H})$ such that for all $r_{H}\in P_{0}(\mathcal{T}_{H})$,

[TABLE]

Here ${\sigma}^{n}$ is an upscaled total flux quantity approximating (over face $F$) ${\sigma^{n}}|_{F}\approx-{\bf n}_{F}\cdot\lambda(s)K\nabla(u+g)$; the sets $\mathcal{F}_{H,I}$, $\mathcal{F}_{H,\rm{out}}$, and $\mathcal{F}_{H,\rm{in}}$ contain interior faces, Dirichlet boundary faces with outgoing and ingoing flux, respectively; $s_{B}$ is the saturation boundary condition; and $s_{H,\rm{upw}}^{n}$ is the upwind saturation

[TABLE]

where $T_{1}$ and $T_{2}$ are adjacent to $F$ and ${\bf n}_{F}$ points from $T_{1}$ to $T_{2}$; and the function $\psi(s)=\lambda_{w}(s)/\lambda(s)$ is the so called fractional flow function. The discretization of the pressure equation is: find $u_{h}^{n}\in V_{h}$, so that for all $v_{h}\in V_{h}$,

[TABLE]

where $A^{n}=\lambda(s_{H}^{n-1})K$. As suggested in [22], we define the non-conservative pre-flux ${\bar{\sigma}}^{n}$ be a harmonic average of the (discontinuous) element face flux $-{\bf n}_{F}\cdot A^{n}\nabla(u_{h}+g)$ at the faces. Then we use the post-processing technique presented in the same paper (with non-weighted minimization) to post-process ${\bar{\sigma}}^{n}$ and obtain the conservative flux ${\sigma}^{n}$ used in the saturation equation.

We make two observations on the information exchange between the two equations when using this discretization:

In the pressure equation, we are only interested in the coarse scale saturation $s_{H}^{n-1}$ from the saturation equation. 2. 2.

In the saturation equation, we are only interested in the upscaled flux $\int_{F}\sigma^{n}$ from the pressure equation.

6.2.2 Coarse error indicators

The first observation allows us to compute $E_{u}$, $E_{f}$ and $E_{g}$ from Section 3.5.1, without saving $\tilde{A}_{T}$. Suppose that for element $T$, the lagging coefficient ${\tilde{A}}_{T}=A^{m}$ is from time step $m<n$, and $A=A^{n}$. We note that $\delta_{T}$ is a coarse quantity, since fine-scale $K$ cancels,

[TABLE]

Thus, to compute $\delta_{T}$ in (13), only ${\tilde{\lambda}}_{T,T^{\prime}}:=\lambda(s^{m-1}_{H})|_{T^{\prime}}$ for $T^{\prime}\subset U_{k}(T)$ need to be saved from previous time steps. The memory required to store ${\tilde{\lambda}}_{T,T^{\prime}}$ behaves like $\mathcal{O}(k^{d}H^{-d})$. Also, ${\tilde{\lambda}}_{T,T^{\prime}}$ can be used to compute $A^{-1/2}{\tilde{A}_{T}}^{1/2}$ in (13).

To summarize, no lagging fine scale information needs to be stored. Only the coarse quantities ${\tilde{\mu}}_{T,T^{\prime}}$ and ${\tilde{\lambda}}_{T,T^{\prime}}$ need to be saved between iterations.

6.2.3 Coarse face flux

The second observation is that we only need the coarse element face flux $\int_{F}\sigma^{n}$ for the saturation equation. Since this quantity is defined on the coarse mesh, we precompute

[TABLE]

for all $T,T^{\prime}\in\mathcal{T}_{H}$, for all faces $F\subset\overline{T^{\prime}}$ and all basis functions $\phi_{i}$ with support in $T$. Here, we used the harmonic average $\{\{v\}\}|_{F}=2\frac{(v|_{T_{1}})|_{F}(v|_{T_{2}})|_{F}}{(v|_{T_{1}})|_{F}+(v|_{T_{2}})|_{F}}$, where $T_{1}$ and $T_{2}$ are the two elements adjacent to $F$, if $F$ is an interior face, and $\{\{v\}\}|_{F}=2(v|_{T})|_{F}$, where $T$ is adjacent to $F$, if $F$ is a boundary face. The memory required for storing ${\tilde{\sigma}}^{n}_{u,T,T^{\prime},F,i}$ and ${\tilde{\sigma}}^{n}_{fg,T,T^{\prime},F}$ scales with $\mathcal{O}(k^{d}H^{-d})$, since the two quantities are zero for all pairs $(T^{\prime},T)$ except when $T^{\prime}\subset U_{k}(T)$.

If we now let the coarse component of the multiscale solution be expressed as $\mathcal{I}_{H}{\hat{u}}_{k}=\sum_{i}\alpha_{i}\phi_{i}$, then we can compute the upscaled non-conservative face flux by

[TABLE]

The final conservative face flux $\sigma^{n}|_{F}$ is then computed using the post-processing technique developed in [22].

We conclude this section by listing the upscaling algorithm (Algorithm 2) for this two-phase Darcy flow problem. In this algorithm, the memory requirements are $\mathcal{O}(k^{d}H^{-d}+h^{-d})$ (where $h^{-d}$ is for the coefficient $A$ which can be distributed on different computational nodes). This allows for very refined fine meshes. Also, the coarse element loop is still completely parallel, and this algorithm serves as a good candidate for a scalable memory efficient upscaling algorithm.

6.3 Darcy flow upscaling numerical experiments

In the following two experiments, we investigate the properties of the upscaling algorithm presented in the previous section. We pick the following mobility functions $\lambda_{w}(s)=s^{3}$, $\lambda_{n}(s)=(1-s)^{3}$, and $\lambda(s)=\lambda_{w}(s)+\lambda_{n}(s)$.

6.3.1 2D random field data

We let $\Omega=[0,1]^{2}$, $\Gamma_{D}=\{x\in\partial\Omega:x_{1}=0\text{ or }x_{1}=1\}$, $\Gamma_{N}=\partial\Omega\setminus\Gamma_{D}$, $f=0$, $g=1-x_{1}$. We use a $512\times 512$ rectangular grid as fine mesh for $V_{h}$ and a $64\times 64$ grid as coarse mesh for $V_{H}$. The permeability $K$ is realized as a piecewise constant function on the fine mesh from a lognormal distribution with exponential spatial correlation and standard deviation $3$, i.e. for a fine element midpoint $x_{\rm m}$,

[TABLE]

where $\kappa(x_{\rm m})\sim\mathcal{N}(0,1)$ and covariance between points are

[TABLE]

with correlation length $d=0.05$ and where $\|\cdot\|_{2}$ denotes Euclidian norm. The initial saturation is set to $s^{0}=0$, and boundary conditions are set to $s_{B}=1$ on the left boundary, which is the only boundary with ingoing flux. The number of time steps and their size are set to $N=2000$ and $\Delta t=N^{-1}$, respectively.

The upscaling algorithm was run with this setup with $k=1,2,3$ and ${\rm TOL}=0.4$, 0.2, 0.1, 0.05, 0.025, and 0.0125. A reference solution $s^{n}_{H,\rm ref}$, where the pressure equation was solved on the fine mesh using the $Q_{1}$ standard finite element method in every iteration was computed. To illustrate the need for upscaling, we also computed the pressure equation on the coarse mesh using the standard $\mathcal{P}_{1}$ finite elements. See Figure 8 for plots over the error in the saturaton solution at the final time step and average fraction of recomputed correctors. In the error plot we can see that both parameters $k$ and ${\rm TOL}$ affect the error in the chosen regimes. We note from the recomputation plot that there is no dependency between the fraction of recomputed element correctors and patch size $k$. Figure 9 shows an example of the saturation solution and the number of times correctors have been computed.

6.3.2 3D random field data

We let $\Omega=[0,1]^{3}$, $\Gamma_{D}=\{x\in\partial\Omega:x_{1}=0\text{ or }x_{1}=1\}$, $\Gamma_{N}=\partial\Omega\setminus\Gamma_{D}$, $f=0$, $g=1-x_{1}$. We use a $128\times 128\times 128$ rectangular grid as fine mesh for $V_{h}$ and a $16\times 16\times 16$ grid as coarse mesh for $V_{H}$. We use a sample $\omega_{i,j,k,\ell}$ of independent uniformly distributed random numbers between [math] and $1$. The permeability $K$ is a piecewise constant function on the fine uniform mesh and is defined by

[TABLE]

where $x_{\rm m}=(x_{{\rm m},1},x_{{\rm m},2},x_{{\rm m},3})$ are the fine element midpoints, and $\lceil\cdot\rceil$ denotes the ceiling function. See Figure 10 for the particular realization used. The boundary conditions are set to $s_{B}=0$ on the boundary with ingoing flux, and the initial saturation is set to a piecewise constant function $s^{0}_{H}$ on the coarse mesh with the following values in the coarse element midpoints $x_{\rm m}$,

[TABLE]

The number of time steps are set to $N=200$ and the time step to $\Delta t=1$.

The upscaling algorithm was run for the three parameter combinations I: $k=1$, ${\rm TOL}=0.1$, II: $k=2$, ${\rm TOL}=0.1$, and III: $k=1$, ${\rm TOL}=0.01$. Figure 10 gives an illustration of the solution at $n=0$ and $n=200$. One of the images shows the recomputed elements as blue boxes and we can see that many elements are not recomputed. In this case, we were not able to compute a reference solution using the available computational resources, but we can estimate the sensitivity of the solutions with respect to the parameters $k$ and ${\rm TOL}$. Let $s^{200}_{H,{\rm I}}$, $s^{200}_{H,{\rm II}}$, and $s^{200}_{H,{\rm III}}$ denote the saturation solutions at time step $n=200$ for the parameter combinations I, II, and III, respectively. We get

[TABLE]

These numbers suggest that the error due to localization (controlled by the parameter $k$) dominates in this case.

7 Conclusion

Elliptic equations with similar rapidly varying coefficients occur for instance in time-dependent problems for two-phase Darcy flow and in stochastic simulations on defect composite materials. We consider a sequence of elliptic equations, each with different coefficients $A^{n}$, $n=1,2,\ldots$. We define a method that computes and updates an LOD multiscale space as we iterate through the coefficients. This is done by the computation of localized element correctors that depend on the coefficient in the vicinity of the element. These computations can be performed completely in parallell. We derive error indicators $e_{u,T}$, $e_{f,T}$, and $e_{g,T}$ that indicate whether or not to update the corrector at element $T$ while iterating the sequence of coefficients. By selecting a small enough tolerance ${\rm TOL}$ for the error indicators, the multiscale space will keep its approximation properties through the sequence of coefficients. It is shown analytically and numerically that the error indicators bound the error in energy norm of the solution. We present a memory efficient upscaling algorithm for a particular application of two-phase Darcy flows.