An Efficient Online-Offline Method for Elliptic Homogenization Problems

We consider a prototypical elliptic boundary value problem

where $\varepsilon$ is a small parameter that signifies explicitly the multiscale nature of the problem. We assume that the coefficient $a^{\;\varepsilon}$, which is not necessarily symmetric, belongs to a set $\mathcal{M}(\alpha,\beta,D)$ that is defined by

where $D$ is a bounded domain in $\mathbb{R}^{d}$ and $(\cdot,\cdot)$ denotes the inner product in $\mathbb{R}^{d}$ , while $\lvert\cdot\rvert$ is the corresponding norm.

In the sense of H-convergence [53, 54], for every sequence $a^{\;\varepsilon}\in\mathcal{M}(\alpha,\beta,D)$ and $f\in H^{-1}(D)$, the sequence $u^{\,\varepsilon}$ of the solution to (1.1) satisfies

where $u_{0}$ is the solution of the homogenization problem

and $\mathcal{A}\in\mathcal{M}(\alpha,\beta,D)$. Here $H_{0}^{1}(D),L^{2}(D)$ and $H^{-1}(D)$ are standard Sobolev spaces [4].

The quantities of interest for Problem (1.1) and Problem (1.4) are the homogenized solution $u_{0}$ over the whole domain and the solution $u^{\,\varepsilon}$ at certain critical local region. The former stands for the information at the large scale, and the later mimics the information at small scale. There are lots of work devoted to efficiently compute such quantities during the last several decades; see, e.g., [6, 20, 17], among many others. Presently we are interested in the efficient way to compute $u_{0}$. A typical way that towards this is provided by the heterogeneous multiscale method (HMM) [18, 3], and the FE${}^{2}-$method [40] commonly used in the engineering community that is also in the same spirit of HMM. The underlying idea of this approach is to extract $\mathcal{A}$ by solving the cell problems posed on the sampling points of the macoscopic solver. At each point, one needs to solve $d$ cell problems with $d$ the dimensionality. Therefore, the main computational cost comes from solving all these cell problems. The number of the cell problems grows rapidly when higher-order macroscopic solvers are employed. To reduce the cost, certain nonconventional quadrature schemes were proposed in [16] when finite element method is used as the macroscopic solver. The number of the cell problems reduces to one third compared to the standard mid-point quadrature scheme when $\mathbb{P}_{2}$ Lagrange finite element method is employed as the macroscopic solver. Unfortunately, it does not seem easy to extend such idea to even higher order macroscopic solvers because the quadrature nodes tend to accumulate in the interior of the element [51, 48, 50].

In [35], the authors presented a local least-squares reconstruction of the effective matrix using the solution of the cell problems posed on the vertices of the triangulation, which was dubbed as HMM-LS. The total number of the cell problems equals to the total number of the interior vertices of the triangulation, which is of $\mathcal{O}(h^{-d})$ with $h$ the mesh size of the macroscopic solver. This method achieves higher-order accuracy with almost the same cost of HMM with $\mathbb{P}_{1}$ Lagrange finite element method as the macroscopic solver [19]. A drawback of this method is that the number of the cell problems is still quite large when mesh refinement is necessary. Moreover, if the adaptive strategy is used in the macroscopic solver, then one has to solve many cell problems around the regions with mesh refinement.

In this work, we propose an offline-online method to compute $u_{0}$ efficiently. The main idea is to separate the microscopic solver from the macroscopic solver. In the offline stage, we firstly solve the cell problems posed on a sampling point set and obtain the effective matrix at all these points, then we reconstruct an effective matrix locally by solving the discrete least-squares. The sampling set is constructed from a triangulation of domain, which is usually coarser than the online triangulation. In the online stage, we solve the macroscopic problem with the effective matrix prepared in the offline stage. Such decoupled strategy brings more flexibilities to reduce the number of the cell problems, we may either refine the offline triangulation mesh or increase the reconstruction order, which is guided by the a priori error estimate. The offline computation bears certain similarity with h-p finite element method [7, 47]. Moreover, the offline computation is no longer linked to the macroscopic solvers, which is particularly attractive to higher order macroscopic solver and three dimensional multiscale problems. With the aid of the theoretical results proved in [36, 35], we study the accuracy and the stability of the reconstruction procedure, which is crucial to prove the optimal error estimate of the proposed method. As illustrated by the numerical tests in § 4, the offline-online method converges with optimal order while the cost is smaller than both HMM and HMM-LS method.

The reduced basis HMM proposed in [1, 2] also employed the offline-online idea. The difference between reduced basis HMM and our method lies in the following points: Firstly they employed reduced basis idea in the offline stage while we construct an independent triangulation upon which the cell problems are solved. Secondly they used the empirical interpolation method [8] while we resort to a local least-squares to reconstruct the effective matrix. Finally, a thorough analysis of the least-squares reconstruction is conducted in our work, which concerns the approximation accuracy and the stability of the reconstruction, such rigorous theoretical results put the method on a firm footing. In addition, the analysis of the least-squares reconstruction is of independent interest for other problems such as the construction of the optimal polynomial admissible meshes [13, 45], the discrete norm for polynomials [46], the approximation of the Fekete points [11] and the discontinuous Galerkin method based on patch reconstruction [36, 37].

The rest of the paper is organized as follows. In § 2, we introduce the offline-online method that is based on a local discrete least-squares reconstruction. We derive the optimal error estimate in § 3, in particular, we prove the discrete least-squares is stable with respect to small perturbation. In § 4, we report numerical examples in two and three dimensions, the coefficient $a^{\;\varepsilon}$ may be locally periodic, almost periodic and random checker-board. To demonstrate the efficiency of the offline-online method, we compare it with HMM and HMM-LS method. In addition, we solve a problem posed on L-shape domain with nonsmooth solution. The conclusions are drawn in the last section.

Throughout this paper, we shall use the standard notations for the Sobolev spaces, norms and semi-norms, cf., [4], e.g.,

For any measurable set $E$, we define the mean of an integrable function $g$ over $E$ as

We shall also use the discrete $\ell_{p}$ norm for any $x\in\mathbb{R}^{n}$ as

Throughout the paper the generic constant $C$ may be different from line to line, while it is independent of $\varepsilon$ and the mesh size parameters $h,H$.

The macroscopic solver is chosen as the standard $\mathbb{P}_{l}$ Lagrange finite element method ($\mathbb{P}_{l}$ FEM) [14]. $\mathbb{P}_{l}$ is defined as the set of polynomials with degree less than $l$ for the sum of all variables. The finite element space is denoted by $V_{h}$ corresponding to the triangulation $\tau_{h}$ with mesh size $h$ that is the maximum of the element size $h_{\tau}$ for all elements $\tau\in\tau_{h}$, where $h_{\tau}$ is the diameter of $\tau$. We assume that all the elements $\tau$ in $\tau_{h}$ satisfy the shape-regular condition in the sense of Ciarlet and Raviart [14], i.e., there exists a constant $\sigma_{0}$ such that $h_{\tau}/\rho_{\tau}\leq\sigma_{0}$, where $\rho_{\tau}$ is the diameter of the smallest ball inscribed into $\tau$, and $\sigma_{0}$ is the so-called chunkiness parameter [12].

The method consists of offline part and online part. In the offline part, we approximate the effective matrix $\mathcal{A}$ as follows.

At each sampling point $x_{K}$, the effective matrix $A_{H}(x_{K})$ is defined by averaging the flux arising from the cell problems:

where the cell $I_{\delta}(x_{K}){:}=x_{K}+\delta Y$ with $Y{:}=(-1/2,1/2)^{d}$ and $\delta$ the cell size. Here for $i=1,\cdots,d$, $v_{i}^{\varepsilon}$ satisfies

In what follows, we supplement some details in the algorithm. The first thing is the construction of the element patch $S(K)$ for any element $K\in\mathcal{T}_{H}$. We start from assigning a threshold value $N_{\text{lowest}}$ that is used to control the size of $S(K)$. There are several different ways to find $S(K)$. One way is to define $S(K)$ in an recursive way as in [35]: For any $t\in\mathbb{N}$, we let

Once $\#S_{t}(K)\geq N_{\text{lowest}}$, we stop the construction and let $S(K)=S_{t}(K)$. This means that we add the Moore neighbors [52] to $S(K)$ in a recursive way. Another way is using the Von Neumann neighbor [52], i.e., we include the adjacent edge-neighboring elements into the element patch instead of the Moore neighbor. We refer to [36] and [37, Appendix A] for a detailed description for such construction, while $S(K)$ in all the tests in § 4 are defined as in (2.5). We denote by $\mathcal{I}_{t}(K)$ the set containing all the sampling points that belong to $S_{t}(K)$. In all the tests in § 4, we use the barycentric of each element $K$ as the sampling point. There are also other choices for construction $\mathcal{I}_{t}(K)$. However, the vertices are not preferred because the communication cost is quite high for three dimensional problems, though this is a good choice for two dimensional problems; cf., [35]. In what follows, we may drop the subscripts $t$ in $S_{t}(K)$ and $\mathcal{I}_{t}(K)$ when there is no confusion may occur.

We remark that there are many other variants for the definitions of the cell problems and the effective matrix in the literatures; see e.g., [59, 24]. The periodic cell problem will be discussed in § 4. As illustrated by the examples in § 4, the overhead caused by the local least-squares is small compared to the computational cost of solving all the cell problems, though the number of the local least-squares is also the same with the cell problems.

To study the convergence of the method, we define

where $\|\,B\,\|_{F}$ is the Frobenius norm of a $d-$by$-d$ matrix $B$.

Similar to [35, Lemma 3.1, Lemma 3.2], the following lemma gives the error estimates of the proposed method.

The proof of the above lemma follows the same line of [19, Theorem 1.1] except the explicit constants in (3.1), we omit the proof.

To further elucidate the error structure of the method, we define the reconstruction operator for any piecewise constant function $v$ defined on $\mathcal{T}_{H}$ as follows. Let $\mathcal{R}_{K}v$ be the solution of the least-squares

We imbed $\mathcal{R}_{K}$ into a global operator as $\mathcal{R}|_{K}=\mathcal{R}_{K}$.

Given the reconstruction operator $\mathcal{R}$, we decompose $e(\text{MOD})$ as

where $\overline{A}_{H}(x)$ is a piecewise matrix posed on $\mathcal{T}_{H}$ that is defined by

with $A_{H}(x_{K})$ given by (2.2) and (2.3).

The first term in the right-hand side of (3.4) is the reconstruction error while the second one is the so-called estimation error. To quantify these two terms, we need some properties of the reconstruction operator, which is by now well-understood by virtue of [35] and [36]. To state such properties, we make two assumptions on $S(K)$ and $\mathcal{I}(K)$.

This assumption concerns the geometry of $S(K)$, which is crucial for the uniform boundedness of $\mathcal{R}$. The motivation for this assumption lies in the following Markov inequality [38]:

Here $\|\,\nabla g\,\|_{L^{\infty}(S(K))}{:}=\max_{x\in\overline{S(K)}}\|\,\nabla g(x)\,\|_{\ell_{2}}$. This inequality is proved in [36, Lemma 5], which is a combination of [57, Proposition 11.6] and the fact that $S(K)$ satisfies the uniform interior cone condition [4], which is a direct consequence of Assumption A.

In [35], the authors make the following assumption on $S(K)$.

By [26, Lemma 1.2.2.2 and Corollary 1.2.2.3], Assumption A’ implies Assumption A. Under Assumption A’, Wilhelmsen [58] proved the following Markov inequality:

where $w(K)$ is the width of $S(K)$, which is the minimum distance between parallel supporting hyperplanes of $S(K)$.

The above Markov inequalities may be viewed as a type of inverse inequality. By the classical inverse inequality for p-finite element method [47, Theorem 4.76] and a simple scaling argument, we may conclude that there exists $C$ independent of the diameter of $S(K)$ but depends on the shape of $S(K)$ such that for all $g\in\mathbb{P}_{m}(S(K))$,

The index $2$ is sharp if $S(K)$ is a locally Lipschitz domain. However, for a cuspidal domain, Kroó and Szabodos [34] proved that if $S(K)$ is a Lip$\gamma$-domain with $0<\gamma<1$¹¹1A typical Lip$\gamma$-domain is the unit $\ell_{\gamma}-$ball, i.e.,$\{\,x\in\mathbb{R}^{d}\,\mid\,\lvert x_{1}\rvert^{\gamma}+\cdots+\lvert x_{d}\rvert^{\gamma}\leq 1\,\}$ with $0<\gamma<1$., then there exists $C$ depends on the shape of $S(K)$ such that

where the index $2/\gamma$ is sharp. This would require a larger patch to ensure the reconstruction accuracy, and the reconstruction is less stable than the patch satisfying either Assumption A or Assumption A’. Moreover, the constant $C$ in (3.7) is only known for $S(K)$ with special shape in the literature, e.g., $S(K)$ is an interval [46], and $S(K)$ is a simplex [56, 33]. All the prefactors are important for us to derive the realistic conditions that ensure the uniform boundedness of $\Lambda(m,\mathcal{I}(K))$; cf., Lemma 3.3.

On the other hand, The authors in [35] derived an explicit expression of $C$ that depends on the recursion depth $t$ and the chunkiness parameter $\sigma$ under Assumption A’. It seems that the convexity of the patch is rather restrictive in implementation, particularly for an $L$-shape domain. Assumption A is less restrictive and is easy to check in practice.

This assumption concerns the cardinality of the sampling set $\mathcal{I}(K)$, which gives the uniqueness and hereby the existence of the solution of the discrete least-squares (2.1). Assumption B requires that the cardinality of $\mathcal{I}(K)$ is at least $\binom{m+d}{d}$ to ensure the unisolvence of the discrete least-squares. A quantitative version of this assumption is

with

In practical implementation, the positions of the sampling nodes may be slightly perturbed due to the measure error or certain uncertainties [15]. A natural question arises whether the reconstruction is robust with respect to the uncertainties. We shall prove that the reconstruction is stable with respect to small perturbation.

Next we prove some properties of the reconstruction operator $\mathcal{R}_{K}$.

The above lemma except the perturbation estimate (3.12) is proved in [35, Theorem 3.3], which is crucial for the accuracy of the reconstruction operator $\mathcal{R}$, more refined estimates on the accuracy of $\mathcal{R}$ may be found in [36, Lemma 4]. The perturbation estimate (3.12) shows that the set of the sampling nodes $\mathcal{I}(K)$ perturbed a little bit remains a norming set with a slightly bigger upper bound, i.e., for all $\delta\in(0,1)$,

This perturbation estimate has been encapsulated in an abstract form in [36, Lemma 2], while there is no proof for the perturbed reconstruction operator $\widetilde{R}$.

The following lemma ensures the uniform boundedness of $\Lambda(m,\mathcal{I}(K))$.

The first statement is proved in [36, Lemma 5], we only prove the second statement (3.15), which improves [35, Lemma 3.5].

If Assumption A holds, then we usually have $R\simeq tH_{K}$, the estimate (3.14) suggests that $r\simeq m\sqrt{t}H_{K}$ implies the uniform boundedness of $\Lambda(m,\mathcal{I}(K))$. This means that $S(K)$ cannot be too narrow in certain directions. If Assumption A’ holds, then the estimate (3.15) shows that $w(K)\simeq m^{2}H_{K}$, which immediately implies $t\simeq m^{2}$. Both conditions show that a relative large patch is required for the reconstruction. Furthermore, if $S(K)$ is a cuspidal domain, then we may use (3.8) to prove that if

with a constant $c_{0}$ depending on the shape of $S(K)$, then the bound (3.14) is also valid. This condition indicates that $S(K)$ has a recursive depth $t\simeq m^{\gamma}$ with $\gamma>2$. An even larger patch is required to ensure the stability. This may occur for domain $D$ with complicated boundary or rough boundary.

It remains to find an upper bound for $\#\mathcal{I}(K)$. This is a direct consequence of Assumption A or Assumption A’ and the shape regularity of $\mathcal{T}_{H}$.

It is clear that the upper bound (3.16) is independent of the way for construction $S(K)$. For the two ways based on either the Moore neighbor or the von Neumann neighbor, we have $\#\mathcal{I}(K)\simeq t^{d}$ with $t$ the recursion depth, which is consistent with the corresponding upper bound proved in [35, Lemma 3.4] and [36, Lemma 6], in which the two-dimensional problem has been dealt with. Both Lemma 3.3 and Lemma 3.4 require that $\#\mathcal{I}(K)$ should be quite large, equivalently, $S(K)$ is large, so that the uniform boundedness of $\Lambda(m,\mathcal{I}(K))$ is valid. In numerical tests below, we observe that the method still works quite well even when $\#\mathcal{I}(K)$ is far less than the theoretical threshold. We refer to [35] for a list of the size of $S(K)$ and the upper bound of $\max_{K\in\mathcal{T}_{H}}\Lambda(m,\mathcal{I}(K))$.

Based on the above three lemmas, we are ready to estimate (3.4).

It is worthwhile to mention that $e_{1}(\text{MOD})$ is the so-called estimating error in HMM [18]. There are many works devoted to bounding $e_{1}(\text{MOD})$ and developing new algorithms to improve the estimates; see; e.g., [19, 59, 10, 23, 5, 24, 43] and the references therein.

Substituting the estimate (3.17) into Lemma 3.1, we obtain the main result of this part.