¹

Keywords: numerical homogenization, stochastic homogenization, super-localization, quantitative theory, error estimates, uncertainty
AMS subject classifications: 35R60, 65N12, 65N15, 65N35, 73B27

1 Introduction↩︎

This paper presents a novel numerical stochastic homogenization method for the prototypical random diffusion problem \[-\mathrm{div}(\boldsymbol{A}\nabla \boldsymbol{u}) = f\] subject to homogeneous Dirichlet boundary conditions. Microscopic features of the problem are encapsulated in the random diffusion coefficient \(\boldsymbol{A}\). In this paper, we are particularly interested in coefficients with small correlation lengths.

For deterministic coefficients, numerical homogenization techniques have been studied extensively in the last decades. For recent monographs and reviews on this topic, we refer to [1]–[5]. The random case has not received a similar attention. However, there are several numerical approaches. Let us mention, for example, MsFEM-based numerical stochastic homogenization methods that assume so-called weakly random coefficients; see [6] for an overview. Furthermore, there are many approaches to approximate the effective coefficient in stochastic homogenization, such as [7]–[11]. Closely related to the present work is the numerical stochastic homogenization method proposed and analyzed in [12], [13]. This method is based on the Localized Orthogonal Decomposition (LOD) introduced in [14], [15]; see also [3], [4] for an overview. More precisely, using a reformulation of the LOD for deterministic problems based on a quasi-local discrete integral operator, as discussed in [16], one can derive an effective model of the problem at hand by taking the expectation. This effective model is deterministic, and its solution gives an accurate coarse scale approximation to the expected value of the solution. However, the method is strongly tied to linear finite elements on simplicial meshes (with piecewise constant gradients), which seems to be an artificial limitation and practically unfavorable, especially since structured deterministic and random diffusion coefficients are often based on Cartesian meshes.

In addition to this technical shortcoming, the localization technique underlying the LOD framework has recently been improved, evolving into the Super-Localized Orthogonal Decomposition (SLOD) introduced in [17] (see also [18]–[22]). While the LOD has exponentially decaying basis functions that lead to an exponentially decaying localization error with respect to the diameter of the basis supports relative to the underlying mesh, the localization error of the SLOD actually decays super-exponentially, cf. [17]. This results in smaller local fine-scale problems when computing the basis functions and increased sparsity for a given tolerance.

In this paper, we propose a computationally simple and efficient numerical stochastic homogenization method based on a special collocation-type formulation of the SLOD. This collocation-type formulation leads to a coarse stiffness-type matrix that can be assembled without any communication between the basis functions defined on the patches of the coarse mesh. This allows each patch to be considered separately for sampling, allowing for improved parallelization and a significant speed-up of the method’s assembly process. Furthermore, the favorable localization properties of the SLOD allow for a computationally efficient sampling procedure. In the case of a random diffusion coefficient with a small correlation length and under standard assumptions of quantitative stochastic homogenization, this paper provides an error estimate for the coarse scale approximation of the proposed method, where certain SLOD-specific quantities contribute in an a posteriori manner. The proof of this error estimate is based on the theory of quantitative stochastic homogenization; see, e.g., [23]–[26]. Classical LOD-techniques [3], [14], [15] are used to further evaluate these SLOD-specific quantities, and a worst-case a priori error analysis is conducted for one of them. Several numerical experiments are performed to quantitatively study the effect of the correlation length and other discretization parameters on the accuracy of the approximation.

This manuscript is structured as follows. First, in 2 we state the model problem in weak form. Then, in 3, we introduce the novel numerical stochastic homogenization method. An a posteriori error analysis of the method is performed in 4. 5 derives a worst-case estimate for the quantity appearing in the a posteriori error bound. Practical aspects of implementation are addressed in 6. Finally, 7 provides numerical experiments that underline the theoretical results of this paper.

2 Model problem↩︎

We consider the model problem \[\label{e:modelstrong} \left\{ \begin{align} -\mathrm{div} (\boldsymbol{A}(\omega)(x) \nabla \boldsymbol{u}(\omega)(x) ) &= f(x),&x\in D\\ \boldsymbol{u}(\omega)(x)&=0,&x\in \partial D \end{align} \right\}\quad\text{for almost all }\omega\in\Omega,\tag{1}\] where \((\Omega,\mathcal{F},\mathbb{P})\) denotes the underlying probability space, \(f \in L^2(D)\) is a deterministic right-hand side and \(D\) is a \(d\)-dimensional bounded Lipschitz polytope with \(d\in\{1,2,3\}\). Without loss of generality, we assume that \(D\) is scaled to unit size. Suppose that \(\boldsymbol{A}\) is a \(\mathbb{R}^{d\times d }\)-valued pointwise symmetric Bochner measurable function, which is uniformly elliptic and bounded, i.e., there exist \(0<\alpha \leq \beta < \infty\) such that for almost all \(\omega \in \Omega\) \[\label{eq:coeffspd} \alpha |\xi|^2 \leq \langle \xi,\boldsymbol{A}(\omega)(x)\xi \rangle \leq \beta |\xi|^2\tag{2}\] holds for all \(\xi \in \mathbb{R}^d\) and almost all \(x \in \mathbb{R}^d\), where \(\langle\cdot,\cdot\rangle\) denotes the Euclidean inner product of \(\mathbb{R}^d\) and \(|\cdot|\) its induced norm.

The weak formulation of the model problem 1 seeks a \(H^1_0(D)\)-valued random field \(\boldsymbol{u}\) such that for almost all \(\omega\in\Omega\) it holds that \[\label{e:modelweak} \boldsymbol{a}_\omega(\boldsymbol{u}(\omega),v):={( \boldsymbol{A}(\omega)\nabla\boldsymbol{u}(\omega)\,,\,\nabla v )}_{L^2(D)}= {( f\,,\,v )}_{L^2(D)}\quad\text{for all }v\in H^1_0(D).\tag{3}\] Here, \({\left( \cdot\,,\,\cdot \right)}_{L^2(D)}\) denotes the inner product on \(L^2(D)\) or \((L^2(D))^d\).

Subsequently, we introduce a shorthand notation for norms and inner products of Bochner spaces. Let \(X\) be a Hilbert space equipped with the inner product \((\cdot,\cdot)_X\). In this case, the Bochner space \(L^2(\Omega;X)\), denoting the space of \(X\)-valued random fields with finite second moments, is also a Hilbert space with the inner product \[(\boldsymbol{v},\boldsymbol{w})_{L^2(\Omega;X)}\mathrel{\vcenter{:}}= {\mathbb{E}}\big[(\boldsymbol{v}(\omega),\boldsymbol{w}(\omega))_{X}\big].\] We write \(\|\cdot\|_{L^2(\Omega;X)}^2 \mathrel{\vcenter{:}}= (\cdot,\cdot)_{L^2(\Omega;X)}\) for the induced norm of this inner product.

Under the given assumptions, the bilinear from \(\boldsymbol{a}_\omega\) depends continuously on \(\boldsymbol{A}\) and, in particular, is measurable as a function of \(\omega\). Hence, the above problem can be reformulated in the Hilbert space \(L^2(\Omega;H^1_0(D))\). The Lax–Milgram theorem then proves its well-posedness, i.e., there exists a unique solution \({\boldsymbol{u}\in L^2(\Omega;H^1_0(D))}\) satisfying \[\label{eq:stabilitysol} \|\nabla\boldsymbol{u}\|_{L^2(\Omega;L^2(D))} \leq \alpha^{-1}C_{\mathrm F}\|f\|_{L^2(D)},\tag{4}\] where \({\| \cdot \|}_{L^2(D)}\) denotes the \(L^2(D)\)-norm and \(C_F\) is the Friedrichs constant of \(D\).

3 Numerical stochastic homogenization method↩︎

The construction of the novel stochastic homogenization method is based on ideas of the SLOD introduced in [17]. In the deterministic setting, the SLOD identifies an almost local basis of the space obtained by applying the solution operator to \(\mathbb{P}^0(\mathcal{T}_H)\), the space of piecewise constants with respect to some coarse mesh \(\mathcal{T}_H\) of the domain \(D\). This is achieved by identifying piecewise constant right-hand sides supported on patches of the mesh \(\mathcal{T}_H\) such that their responses under the corresponding localized solution operator have a minimal conormal derivative. These localized responses are then used as the basis functions of the SLOD. In the stochastic setting, an adaptation of this approach is required, which involves identifying deterministic local source terms such that the conormal derivative of the localized responses is small in expectation.

In the following we assume that the considered family of meshes \(\{\mathcal{T}_H\}_{H}\) is quasi-uniform and consists of meshes with closed, convex and shape-regular elements. The parameter \(H>0\) specifies the maximal element diameter of the mesh \(\mathcal{T}_H\). We denote by \(\Pi_H\colon L^2(D) \rightarrow \mathbb{P}^0({\mathcal{T}}_H)\) the \(L^2\)-orthogonal projection onto \(\mathbb{P}^0({\mathcal{T}}_H)\). Let us also give a precise definition of the concept of patches with respect to \(\mathcal{T}_H\). The first-order patch \({\mathsf N}(S) = {\mathsf N}^1(S)\) of a union of elements \(S \subset D\) is defined by \[{\mathsf N}^1(S) \mathrel{\vcenter{:}}= \bigcup \{T\in{\mathcal{T}}_H \,\colon\,T\cap S \neq \emptyset \}.\] For any \(\ell = 2,3,4,\dots\), the \(\ell\)-th order patch \({\mathsf N}^\ell(T)\) of \(T\) is then given recursively by \[\label{eq:patch} {\mathsf N}^\ell(T) \mathrel{\vcenter{:}}= {\mathsf N}^1({\mathsf N}^{\ell-1}(T)).\tag{5}\]

The following derivation of the basis functions considers a fixed element \(T\in {\mathcal{T}}_H\) and oversampling parameter \(\ell \in \mathbb{N}\), where we assume that the patch \(D_T \mathrel{\vcenter{:}}= {\mathsf N}^\ell(T)\) does not coincide with the whole domain. We denote the deterministic source term corresponding to \(T\) by \(g_T \in \mathbb{P}^0({\mathcal{T}}_{H,D_T})\) with the submesh \({\mathcal{T}}_{H,D_T} \mathrel{\vcenter{:}}= \{K \in {\mathcal{T}}_H \,\colon\,K \subset D_T\}\). The global response \(\boldsymbol{\varphi}_T\in L^2(\Omega; H^1_0(D))\) to \(g_T\) is defined for almost all \(\omega \in \Omega\) by \[\label{eq:basisglob} \boldsymbol{a}_\omega(\boldsymbol{\varphi}_T(\omega),v) = {( g_T\,,\,v )}_{L^2(D)}\quad\text{for all }v\in H^1_0(D).\tag{6}\] Its localized version \(\boldsymbol{\varphi}_T^\mathrm{loc}\in L^2(\Omega; H^1_0(D_T))\) is for almost all \(\omega \in \Omega\) defined by \[\label{eq:basisloc} \boldsymbol{a}_\omega(\boldsymbol{\varphi}_T^\mathrm{loc}(\omega),v) = {( g_T\,,\,v )}_{L^2(D_T)}\quad\text{for all }v\in H^1_0(D_T).\tag{7}\] From now on, the dependence of stochastic variables is only indicated by a bold symbol for better readability.

To define the conormal derivative of the localized basis function \(\boldsymbol{\varphi}_T^\mathrm{loc}\), we need to introduce some preliminaries on traces and extensions. We denote by \(H^1_\Gamma(D_T)\) the complete subspace of \(H^1(D_T)\) consisting of functions with trace zero at the boundary segment \(\Gamma \mathrel{\vcenter{:}}= \partial D_T \cap \partial D\). Furthermore, let \[\begin{align} {\operatorname{tr}}\colon H^1_\Gamma(D_T)\to X \mathrel{\vcenter{:}}= \operatorname {im} {\operatorname{tr}}\subset H^{1/2}(\partial D_T) \end{align}\] denote the classical trace operator restricted to \(H^1_\Gamma(D_T)\). As an extension operator, we henceforth consider the \(\boldsymbol{A}\)-harmonic extension operator \(\boldsymbol{\operatorname{tr}}^{-1}\colon L^2(\Omega;X) \to L^2(\Omega; H^1_\Gamma(D_T))\) defined as follows: For almost all \(\omega \in \Omega\) and for any given \(\boldsymbol{b}\in L^2(\Omega;X)\), we set \(({\operatorname{tr}}\boldsymbol{\operatorname{tr}}^{-1}\boldsymbol{b})(\omega) = \boldsymbol{b}(\omega)\) and demand that \[\begin{align} \label{eq:trinv} \boldsymbol{a}(\boldsymbol{\operatorname{tr}}^{-1}\boldsymbol{b}, v ) = 0\quad \text{for all } v\in H^1_0(D_T). \end{align}\tag{8}\] The space of locally \(\boldsymbol{A}\)-harmonic functions satisfying homogeneous Dirichlet boundary conditions on \(\Gamma\) can then be defined as \[\begin{align} \boldsymbol{Y}\mathrel{\vcenter{:}}= \boldsymbol{\operatorname{tr}}^{-1}L^2(\Omega;X) \subset L^2(\Omega; H^1_\Gamma(D_T)). \end{align}\] For more details on trace- and extension operators we refer, e.g., to [27].

Combining 7 8 yields for almost all \(\omega \in \Omega\) and all \(\boldsymbol{v}\in L^2(\Omega;H^1_0(D))\) the identity \[\boldsymbol{a}(\boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) = \boldsymbol{a}(\boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}- \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v}) = (g_T, \boldsymbol{v}- \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})_{L^2(D_T)},\] where we used that \((\boldsymbol{v}- \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})|_{D_T} \in L^2(\Omega;H^1_0(D_T))\). With this identity, the definition of \(\boldsymbol{\varphi}_T\) in 6 , and \(\operatorname{supp} g_T \subset D_T\), it follows that \[\begin{align} \boldsymbol{a}(\boldsymbol{\varphi}_T- \boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) = (g_T, \boldsymbol{v})_{L^2(D_T)} - \boldsymbol{a}(\boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) = (g_T, \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})_{L^2(D_T)}. \end{align}\] Taking the expectation, we obtain for any \(\boldsymbol{v}\in L^2(\Omega;H^1_0(D))\) that \[\begin{align} \label{eq:loc-error} {\mathbb{E}}\big[\boldsymbol{a}(\boldsymbol{\varphi}_T- \boldsymbol{\varphi}_T^\mathrm{loc}, \boldsymbol{v}) \big] = {\mathbb{E}}\big[(g_T, \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v})_{L^2(D_T)}\big] = \left(g_T, \Pi_H\; {\mathbb{E}}[ \boldsymbol{\operatorname{tr}}^{-1}{\operatorname{tr}}\boldsymbol{v}]\right)_{L^2(D_T)}. \end{align}\tag{9}\] As a consequence, the (almost) \(L^2\)-orthogonality of \(g_T\) to the space \({\mathbb{E}}[\boldsymbol{Y}] \subset H^1_\Gamma(D_T)\) leads to a small expected localization error for the basis functions \(\boldsymbol{\varphi}_T^\mathrm{loc}\).

Therefore, we can obtain an optimal choice of \(g_T\) by performing a singular value decomposition (SVD) of the compact operator \((\Pi_{H, D_T} \circ {\mathbb{E}})\vert_{\boldsymbol{Y}}: \boldsymbol{Y}\rightarrow \mathbb{P}^0 ({\mathcal{T}}_{H, D_T})\) restricted to the complete subspace \(\boldsymbol{Y}\). Note that the rank of \((\Pi_{H, D_T} \circ {\mathbb{E}})\vert_{\boldsymbol{Y}}\) is less than or equal to \(N \mathrel{\vcenter{:}}= \# {\mathcal{T}}_{H,D_T}\). Hence, the SVD is given by \[\label{eq:svd} (\Pi_{H, D_T} \circ {\mathbb{E}})\vert_{\boldsymbol{Y}} \; \boldsymbol{v}= \sum_{k=1}^{N} \sigma_k (\boldsymbol{v}, \boldsymbol{w}_k)_{L^2(\Omega; H^1(D_T))} g_k\tag{10}\] with singular values \(\sigma_1 \geq \dots \geq\sigma_N\geq 0\), \(L^2(\Omega; H^1(D_T))\)-orthonormal right singular vectors \(\boldsymbol{w}_1,\dots,\boldsymbol{w}_N\) and \(L^2(D_T)\)-orthonormal left singular vectors \(g_1,\dots, g_N\). The choice \(g_T=g_N\) as the left singular vector corresponding to the smallest singular value \(\sigma_N\) is optimal in the sense that \[g_N \in \underset{g\in\mathbb{P}^0({\mathcal{T}}_{H,D_T})\;:\; \Vert g \Vert_{L^2(D_T)}=1}{\text{argmin}} \quad \underset{\boldsymbol{v}\in \boldsymbol{Y}\; :\; \Vert \boldsymbol{v}\Vert_{L^2(\Omega;H^1(D_T))}=1}{\sup} (g,{\mathbb{E}}[\boldsymbol{v}])_{L^2(D_T)}.\] The corresponding smallest singular value \(\sigma_N\) is a measure of the (quasi-)orthogonality between \(g_T\) and \({\mathbb{E}}[\boldsymbol{Y}]\). We hence define \[\label{eq:32sigma} \sigma_T(H,\varepsilon,\ell) \mathrel{\vcenter{:}}= \sigma_N = \underset{\boldsymbol{v}\in\boldsymbol{Y}\; :\; \Vert \boldsymbol{v}\Vert_{L^2(\Omega;H^1(D_T))}=1}{\sup} (g_N,{\mathbb{E}}[\boldsymbol{v}])_{L^2(D_T)}.\tag{11}\] We emphasize that the practical implementation of the SVD in 10 is difficult due to the stochasticity involved; for a practical implementation based on sampling, see 6. For the error analysis in the following section, we introduce the quantity \[\label{eq:locerrind} \sigma\mathrel{\vcenter{:}}=\sigma(H,\varepsilon,\ell) \mathrel{\vcenter{:}}= \max_{T\in \mathcal{T}_H} \sigma_T(H,\varepsilon,\ell),\tag{12}\] which is an indicator for the overall localization error.

Given that, in expectation, \(\boldsymbol{\varphi}_T^\mathrm{loc}\) closely approximates the response of the global solution operator applied to \(g_T\), it is reasonable to define the approximation of a non-Galerkin, collocation-type numerical stochastic homogenization method by \[\label{eq:stochhommethod} \bar{u}_{H,\ell}\mathrel{\vcenter{:}}= \sum_{T\in{\mathcal{T}}_{H}} c_T \Pi_H\mathbb{E}[\boldsymbol{\varphi}_T^\mathrm{loc}],\tag{13}\] where \((c_T)_{T\in {\mathcal{T}}_{H}}\) are the coefficients of the expansion of \(\Pi_Hf\) in terms of the basis functions \(\{g_T \,\colon\,T \in \mathcal{T}_H\}\). An illustration of the deterministic basis functions \({\mathbb{E}}[\boldsymbol{\varphi}_T^\mathrm{loc}]\) can be found in 1. Error estimates for this method are derived in the following sections.

4 Error analysis↩︎

In this section, we perform an error analysis of the proposed stochastic homogenization method based on results from the theory of quantitative stochastic homogenization. This theory requires structural conditions on the randomness of the coefficient field \(\boldsymbol{A}\). For simplicity, the conditions are formulated for coefficient fields defined on \(\mathbb{R}^d\). Hence, the following assumptions implicitly assume that the coefficient field is defined on the full space \(\mathbb{R}^d\). A random field defined on the bounded domain \(D\) can be obtained by restriction.

We emphasize that these conditions are standard in the theory of quantitative stochastic homogenization; see, e.g., [26].

The error bound presented in this section is an a posteriori bound including the constant \(\sigma\) from 12 and the Riesz stability constant of the local source terms \(\left\{g_T\,\colon\,T\in\mathcal{T}_H\right\}\), which quantifies their linear independence. Both constants can be computed a posteriori as outlined in 6. Additionally, we provide a worst-case a priori upper bound on \(\sigma\) in 5. Note that in a practical implementation, the Riesz stability of the local source terms can be ensured as outlined in 6 or [17].

To handle the stochasticity in the error analysis, we need to estimate the variance of the random variables \((\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}\) for any \(T,K \in \mathcal{T}_H\), where \(\mathbb{1}_K\) denotes the indicator function of the element \(K\). To achieve this, we employ the spectral gap inequality 14 from 1. The following lemma provides a representation of the Fréchet derivative of \((\boldsymbol{\varphi}_T^\mathrm{loc}, \mathbb{1}_K)_{L^2(K)}\), a crucial element for this particular step.

Another ingredient in the error analysis is the following regularity result for the localized basis functions. The result is needed to further estimate the term we get after applying the spectral gap inequality. The proof of this result relies on the condition that the patches take the form of \(d\)-dimensional bricks, cf. 5. This condition can be guaranteed, for example, by considering a brick-shaped domain equipped with a Cartesian mesh.

The following theorem encapsulates the main result of this work, giving an a posteriori error bound for the proposed numerical stochastic homogenization method.

5 Error analysis using LOD techniques↩︎

This section utilizes LOD theory to derive an upper bound for the quantity \(\sigma\) that appears in the error estimate from 2. We further estimate \(C_\mathrm{rb}\) for the choice of LOD basis functions made for the upper bound on \(\sigma\).

5.1 Localization error indicator↩︎

We first derive an upper bound for the localization error \(\sigma\) defined in 12 . The bound is based on the lowest-order LOD from [30]–[32], whose construction uses non-negative bubbles \(\{b_T \,\colon\,T \in \mathcal{T}_H\}\); see also [33]. The bubble function \(b_T\in H_0^1(T)\) is chosen such that \(\Pi_Hb_T = \mathbb{1}_T\) and \[\label{eq:bubble} \pi \Vert b_T \Vert_{L^2(T)} \leq H \Vert \nabla b_T \Vert_{L^2(T)} \leq C_\mathrm{b} \sqrt{\vert T \vert}\tag{22}\] holds, where \(C_\mathrm{b}>0\) is independent of \(H\). Recalling the abbreviation \(D_T = \mathsf N^\ell(T)\) for the \(\ell\)-th order patch around \(T\), cf. 5 , we introduce the space of fine-scale functions supported on \(D_T\) by \(\mathcal{W}_{T,\ell} \mathrel{\vcenter{:}}= \{w\in H_0^1(D_T) \,\colon\,\Pi_H\vert_{D_T} w = 0\}\). The LOD basis function corresponding to the element \(T\in\mathcal{T}_H\) is then defined by \[\label{eq:lodbasis} \boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}\mathrel{\vcenter{:}}= (1-\boldsymbol{\mathcal{C}}_{T,\ell})b_T \in L^2(\Omega; H^1_0(D_T)),\tag{23}\] where \(\boldsymbol{\mathcal{C}}_{T,\ell}b_T \in L^2(\Omega;\mathcal{W}_{T,\ell})\) denotes the fine-scale correction of the bubble \(b_T\), which is defined for almost all \(\omega \in \Omega\) by \[\label{eq:corrc} \boldsymbol{a}(\boldsymbol{\mathcal{C}}_{T,\ell}b_T, w) = \boldsymbol{a}(b_T,w) \quad \text{ for all } w\in\mathcal{W}_{T,\ell}.\tag{24}\]

The bound on \(\sigma\) in the following lemma is based on the observation that the LOD basis function \(\boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}\) possesses a \(\mathcal{T}_H\)-piecewise constant source term \[\boldsymbol{g}_{T,\ell}^\mathrm{LOD}\mathrel{\vcenter{:}}= - \mathrm{div} \boldsymbol{A}\nabla \boldsymbol{\varphi}_{T,\ell}^\mathrm{LOD}\in L^2(\Omega;\mathbb{P}^0(\mathcal{T}_{H,D_T}));\] see, e.g., [17].

Combining this a priori result for \(\sigma\) with 2 yields the error estimate given in the following corollary. The Riesz constant \(C_\mathrm{rb}\) can be computed a posteriori, cf. 6.

5.2 Riesz stability↩︎

In a next step, we show that the local source terms corresponding to the LOD basis functions 23 are Riesz stable in the sense of 2.

6 Practical implementation↩︎

To effectively implement the proposed numerical stochastic homogenization method, it is crucial to employ an efficient sampling strategy for the space \(\boldsymbol{Y}\) and ensure that the local source terms \(\{g_T\,\colon\,T\in \mathcal{T}_H\}\) form a stable basis of \(\mathbb{P}^0(\mathcal{T}_H)\). These aspects will be addressed in the following two subsections.

6.1 Sampling of the space \(\boldsymbol{Y}\)↩︎

We consider an arbitrary patch \(D_T\) and denote the number of coarse elements in this patch by \(N \mathrel{\vcenter{:}}= \#{\mathcal{T}}_{H, D_T}\). In a practical implementation, all local infinite-dimensional problems that appear in the derivation of the basis functions must be replaced by finite-dimensional counterparts. To obtain these finite-dimensional counterparts, we perform a discretization using the \(\mathcal{Q}^1\)-finite element method with respect to the fine mesh \(\mathcal{T}_{h,D_T}\) constructed by uniform refinements of \({\mathcal{T}}_{H,D_T}\). The number of elements of \(\mathcal{T}_{h,D_T}\) is denoted by \(n\).

To handle the stochasticity in the definition of \(\boldsymbol{Y}\), our implementation draws \(M\) samples of the random coefficient \(\boldsymbol{A}\) and, for each sample, closely follows the methodology outlined in [17] for the deterministic case. Specifically, we generate a matrix \(\mathbf{S}_i \in \mathbb{R}^{n\times m}\) for \(i=1,\dots,M\), whose columns represent the coordinate vectors of the discrete \(\boldsymbol{A}(\omega_i)\)-harmonic extensions of \(m\in\mathbb{N}\) samples of random boundary data on \(\partial D_T \backslash \partial D\). Then we compute the matrices \(\mathbf{P}_i\in\mathbb{R}^{N\times m}\) by applying the \(L^2\)-orthogonal projection onto the characteristic functions \(\lbrace \mathbb{1}_T\,\colon\,T\in\mathcal{T}_{H,D_T} \rbrace\) column by column to \(\mathbf{S}_i\). Finally, the SVD of the matrix \(\mathbf{X} \mathrel{\vcenter{:}}= [\mathbf{P}_1,\dots,\mathbf{P}_M]\) is computed, yielding coordinate vectors of potential right-hand sides \(g_T\). For details on the practical realization of this SVD, we refer to [17]. Finally, the localized deterministic basis functions are computed as empirical means, again using \(M\) samples of the random coefficient. In the numerical experiments performed in 7, the number of random boundary samples is set to \(m=3N\). For the number of random coefficient samples, we use \(M = 5000\).

6.2 Stable local source terms↩︎

Next, we discuss how the stability of the local source terms \(\{g_T \,\colon\,T \in \mathcal{T}_H\}\) can be ensured in a practical implementation. Our implementation achieves stability by an additional optimization step, similar to the one used in [20]. Given the singular values \(\sigma_1\geq\sigma_2\geq\dots \geq \sigma_N\geq 0\) of the matrix \(\mathbf{X}\) associated with the patch \(D_T\), we consider all indices \(1\leq i \leq N\) such that \[\begin{align} \label{eq:choicelocrhs} \frac{\sigma_i}{\sigma_1}\leq \max \Big\{\Big(\frac{\sigma_N}{\sigma_1}\Big)^{1/p},10^{-10}\Big\} \end{align}\tag{35}\] and denote the resulting set of indices by \(\mathcal{I}\). Each index in the set \(\mathcal{I}\) corresponds to a potential candidate for a local source term. For the choice \(p = 1\) only the smallest singular value is considered. Since our optimization problem is meaningful whenever multiple functions are considered, we restrict ourselves to the choices \(p > 1\).

Among these candidate functions, we choose the one that maximizes a weighted \(L^2(D_T)\)-norm under the unit mass constraint. The weighted \(L^2(D_T)\)-norm is defined using a piecewise constant weighting function that is zero in the central element \(T\) and grows polynomially as the distance from the center increases. This enforces a concentration of mass in the center of each patch, resulting in linearly independent local source terms \(\{g_T \,\colon\,T \in \mathcal{T}_H\}\) in practice. More specifically, we introduce the distance function \(\mathrm{dist}(T, K)\) between the elements \(T,K \in \mathcal{T}_H\) as \[\begin{align} \mathrm{dist}(T, K) \mathrel{\vcenter{:}}= H^{-1}|m_K - m_T | \in \mathbb{N}^d, \end{align}\] where \(m_T , m_K \in \mathbb{R}^d\) are the midpoints of the elements \(T\) and \(K\), respectively. The weighting function is then defined for each element \(K\in{\mathcal{T}}_{H, D_T}\) as \[\begin{align} w_T(K) \mathrel{\vcenter{:}}= \big\vert \mathrm{dist}(T,K)\big\vert_\infty^{r} \end{align}\] for a parameter \(r\geq1\), where \(|\cdot|_\infty\) denotes the infinity norm on \(\mathbb{R}^d\). 2 provides an illustration of this weighting function in two spatial dimensions. In our numerical experiments in 7, we use \(p=1.5\) and \(r = 6\).

7 Numerical experiments↩︎

The following numerical experiments are intended to demonstrate the effectiveness of the proposed numerical homogenization method. In our implementation, we consider uniform Cartesian meshes of the domain \(D = (0,1)^d\) with \(d\in\{1,2\}\). Note that from now on we use \(H\) to denote the side length of the elements instead of their diameter. For the solution of the local patch problems and the computation of the reference solution \(\boldsymbol{u}_h\) we employ the \(\mathcal{Q}_1\)-finite element method on the fine mesh \({\mathcal{T}}_h\) with \(h=2^{-10}\). We denote by \(\bar{u}_{H,h,\ell}\) the fully discrete numerical approximation to \({\mathbb{E}}[\boldsymbol{u}]\). In the following all expected values are replaced by appropriate empirical means.

The random coefficients \(\boldsymbol{A}\) that are considered in the following numerical experiments are piecewise constant with respect to the uniform Cartesian meshes \(\mathcal{T}_\varepsilon\) with mesh sizes \(\varepsilon \in \lbrace 2^{-5},2^{-6},2^{-7},2^{-8},2^{-9}\rbrace\). These coefficients take independent and identically distributed element values in the interval \([0.1,1]\). We further consider the sequence of coarse meshes \(\mathcal{T}_H\) with mesh sizes \(H\in \lbrace 2^{-3},2^{-4},2^{-5},2^{-6}\rbrace\). Note that we only consider coarse mesh sizes \(H>\varepsilon\) for which the coarse mesh does not resolve the minimal length scale of the random coefficient. We also exclude combinations of \(H\) and \(\ell\) for which a patch coincides with the whole domain \(D\). To calculate the reference solution, we employ \(M=5000\) samples, which is consistent with the number used for the local patch problems. The samples are obtained by a quasi-Monte Carlo sampling strategy in one spatial dimension and a Monte Carlo sampling strategy in two spatial dimensions.

Numerical investigation of \(\sigma\) and \(C_\mathrm{rb}\)↩︎

We first examine the behavior of the localization error indicator \(\sigma\) as a function of the coarse mesh size \(H\) and the correlation length \(\varepsilon\). For this, we consider the case \(d = 2\) and utilize the sequences of coarse meshes and correlation lengths mentioned above. 3 visualizes the values of \(\sigma\) for a fixed correlation length \(\varepsilon\) and varying mesh sizes \(H\) (left) and for fixed \(H\) and varying \(\varepsilon\) (right). In both cases one observes a scaling like \(\tfrac{\varepsilon}{H}\), which numerically validates the upper bound for \(\sigma\) from 3 in the case \(d = 2\). Note that the stochastic errors dominate, and consequently, the first term in 25 , which decays exponentially in \(\ell\), is not visible. Plotting \(\sigma\) as a function of \(\ell\) would give a scaling like \(\ell^{-1/2}\).

Next we examine the behavior of the Riesz stability constant \(C_\mathrm{rb}\) of the local SLOD source terms as a function of \(H\). In 4 we observe that \(C_\mathrm{rb}\) scales like \(H^{-4}\), which is consistent with the results for the stochastically averaged LOD source terms proved in 4. Our numerical experiments indicate no dependency of the Riesz stability constant on \(\varepsilon\) or \(\ell\), which is also in line with the findings from 4.

Numerical validation of convergence↩︎

To numerically verify the convergence of the proposed numerical stochastic homogenization method, we consider the source terms \[f_1(x)= 2\pi^2\sin(x),\quad f(x,y)= 2\pi^2\sin(x)\sin(y)\] in one and two spatial dimensions, respectively. 5 6 show the resulting relative \(L^2\)-errors computed using the reference solution \(\boldsymbol{u}_h\). For fixed \(H\) and varying \(\varepsilon\) we observe the rate \(\varepsilon^{d/2}\), which is in agreement with 1. When considering the converse case, we have to distinguish between one and two spatial dimensions. In one dimension, the expected negative power of \(H\) does not manifest itself, and in our numerical experiments the error remains relatively constant with respect to \(H\) (provided the coarse mesh is sufficiently coarse compared to \(\varepsilon\)). In the two-dimensional case, we observe a negative dependence on \(H\), which is much weaker than the \(H^{-2}\) predicted by 1. The error rather seems to scale like \(H^{-1/3}\).

8 Auxiliary results↩︎

The error analysis of the numerical homogenization method proposed in this paper is based on a regularity theory of quantitative stochastic homogenization. It is established within the framework of an equation on the entire space \({\mathbb{R}}^d\) as outlined in [35], thereby extending earlier results from [36], [37]. In this paper we require its subsequent adaptation to the Dirichlet problem on a bounded box. The results were previously used in [13]. Their proof, which is beyond the scope of this manuscript, is analogous to the full-space case and uses the boundary regularity theory of [38] as well as a classical regularity theory at edges and corners.

In the following, we present two results used in the proof of 3. The first provides a \(L^2\)-representation of the Fréchet derivative, which is needed to apply the spectral gap inequality.

The following result is needed to estimate the terms appearing after applying the spectral gap inequality in the proof of 3.

References↩︎