Hierarchical multiscale finite element method for multi-continuum media

TL;DR

This paper introduces a hierarchical multiscale finite element method for efficiently simulating multi-continuum media with multiple scales and high contrast, achieving optimal computational cost while maintaining accuracy.

Contribution

It develops a hierarchical approach that reduces computation cost by using different resolution levels for cell problems across a dense network of macroscopic points.

Findings

The method achieves the same accuracy as full high-resolution solves.

Numerical results validate the analytical convergence and efficiency.

The approach effectively computes effective permeabilities in multi-continuum systems.

Abstract

Simulation in media with multiple continua where each continuum interacts with every other is often challenging due to multiple scales and high contrast. One needs some types of model reduction. One of the approaches is multi-continuum technique, where every process in each continuum is modeled separately and an interaction term is added. Direct numerical simulation in multi scale media is usually not practicable. For this reason, one constructs the corresponding homogenized equations. The paper develops a hierarchical approach for solving cell problems at a dense network of macroscopic points with an essentially optimal computation cost. The method employs the fact that neighboring representative volume elements (RVEs) share similar features; and effective properties of the neighboring RVEs are close to each other. The hierarchical approach reduces computation cost by using different levels of resolution for cell problems at different macroscopic points. The method requires a hierarchy of macroscopic grid points and a corresponding nested approximation spaces with different levels of resolution. Each level of macroscopic points is assigned to an approximation finite element (FE) space which is used to solve the cell problems at the macroscopic points in that level. We prove rigorously that this hierarchical method achieves the same level of accuracy as that of the full solve where cell problems at every macroscopic point are solved using the FE spaces with the highest level of resolution, but at the essentially optimal computation cost. Numerical implementation that computes effective permeabilities of a two scale multicontinuum system via the numerical solutions of the cell problems supports the analytical results. Finally, we prove the homogenization convergence for our multiscale multi-continuum system.

Tables

Table 1. Table 1 : a = 1, the effective coefficients κ 111 ∗ ​ ( x 1 ) superscript subscript 𝜅 111 subscript 𝑥 1 \kappa_{111}^{*}(x_{1}) and κ 211 ∗ ​ ( x 1 ) superscript subscript 𝜅 211 subscript 𝑥 1 \kappa_{211}^{*}(x_{1}) computed by full mesh reference and hierarchical solve along with percentage relative errors between those.

$x_1$	${\kappa }_{111}^{\ast }\left(x_1\right)$			${\kappa }_{211}^{\ast }\left(x_1\right)$
	Full solve	Hierarchical solve	Relative errors(%)	Full solve	Hierarchical solve	Relative errors(%)
$0$	2.8210	2.8210	0.0000	2.8341	2.8341	0.0000
$\frac{1}{16}$	2.8331	2.8266	0.2305	2.8448	2.8432	0.0560
$\frac{1}{8}$	2.8447	2.8406	0.1425	2.8550	2.8524	0.0919
$\frac{3}{16}$	2.8557	2.8590	0.1160	2.8649	2.8653	0.0163
$\frac{1}{4}$	2.8663	2.8639	0.0823	2.8743	2.8735	0.0304
$\frac{5}{16}$	2.8763	2.8688	0.2621	2.8834	2.8813	0.0726
$\frac{3}{8}$	2.8859	2.8886	0.0936	2.8921	2.8936	0.0519
$\frac{7}{16}$	2.8950	2.8996	0.1606	2.9005	2.9016	0.0376
$\frac{1}{2}$	2.9036	2.9036	0.0000	2.9085	2.9085	0.0000
$\frac{9}{16}$	2.9119	2.9076	0.1449	2.9161	2.9151	0.0336
$\frac{5}{8}$	2.9197	2.9176	0.0712	2.9234	2.9218	0.0531
$\frac{11}{16}$	2.9271	2.9317	0.1580	2.9303	2.9317	0.0470
$\frac{3}{4}$	2.9341	2.9350	0.0295	2.9369	2.9373	0.0167
$\frac{13}{16}$	2.9407	2.9382	0.0851	2.9431	2.9428	0.0096
$\frac{7}{8}$	2.9470	2.9484	0.0479	2.9490	2.9501	0.0387
$\frac{15}{16}$	2.9528	2.9557	0.0980	2.9546	2.9552	0.0220
$1$	2.9583	2.9583	0.0000	2.9598	2.9598	0.0000

Table 2. Table 2 : a = .1, the effective coefficients κ 111 ∗ ​ ( x 1 ) superscript subscript 𝜅 111 subscript 𝑥 1 \kappa_{111}^{*}(x_{1}) and κ 211 ∗ ​ ( x 1 ) superscript subscript 𝜅 211 subscript 𝑥 1 \kappa_{211}^{*}(x_{1}) computed by full mesh reference and hierarchical solve along with percentage relative errors between those.

$x_1$	${\kappa }_{111}^{\ast }\left(x_1\right)$			${\kappa }_{211}^{\ast }\left(x_1\right)$
	Full solve	Hierarchical solve	Relative errors(%)	Full solve	Hierarchical solve	Relative errors(%)
$0$	2.8210	2.8210	0.0000	2.8341	2.8341	0.0000
$\frac{1}{16}$	2.8222	2.8215	0.0241	2.8352	2.8350	0.0059
$\frac{1}{8}$	2.8235	2.8230	0.0161	2.8363	2.8360	0.0102
$\frac{3}{16}$	2.8247	2.8250	0.0125	2.8373	2.8374	0.0020
$\frac{1}{4}$	2.8259	2.8256	0.0112	2.8384	2.8383	0.0037
$\frac{5}{16}$	2.8271	2.8261	0.0347	2.8395	2.8392	0.0095
$\frac{3}{8}$	2.8283	2.8288	0.0154	2.8405	2.8408	0.0081
$\frac{7}{16}$	2.8295	2.8302	0.0232	2.8416	2.8418	0.0056
$\frac{1}{2}$	2.8307	2.8307	0.0000	2.8427	2.8427	0.0000
$\frac{9}{16}$	2.8319	2.8313	0.0230	2.8437	2.8435	0.0056
$\frac{5}{8}$	2.8331	2.8327	0.0150	2.8448	2.8445	0.0096
$\frac{11}{16}$	2.8343	2.8352	0.0327	2.8458	2.8461	0.0090
$\frac{3}{4}$	2.8355	2.8357	0.0100	2.8468	2.8469	0.0035
$\frac{13}{16}$	2.8366	2.8363	0.0124	2.8479	2.8478	0.0019
$\frac{7}{8}$	2.8378	2.8382	0.0144	2.9490	2.8492	0.0093
$\frac{15}{16}$	2.8390	2.8396	0.0221	2.8499	2.8501	0.0054
$1$	2.8401	2.8401	0.0000	2.8510	2.8510	0.0000

Equations

$${\mathcal{C}}_{ii}^{ϵ}(x)\frac{\partial u_i^{ϵ}(t,x)}{\partial t}=\text{div}({\kappa }_i^{ϵ}(x)\nabla u_i^{ϵ}(t,x))+Q_i^{ϵ}(u_1^{ϵ}(t,x),...,u_N^{ϵ}(t,x)),\text{in}\Omega$$

$${\mathcal{C}}_{ii}^{ϵ}(x)\frac{\partial u_i^{ϵ}(t,x)}{\partial t}=\text{div}({\kappa }_i^{ϵ}(x)\nabla u_i^{ϵ}(t,x))+Q_i^{ϵ}(u_1^{ϵ}(t,x),u_i^{ϵ}(t,x))+q,i>1,$$

$${\mathcal{C}}_{ii}(x,y)\ge c,{\kappa }_i(x,y)\ge c,Q(x,y)\ge c.$$

$${\mathcal{C}}_{ii}^{ϵ}(x)={\mathcal{C}}_{ii}(x,\frac{x}{ϵ}),{\kappa }_i^{ϵ}(x)={\kappa }_i(x,\frac{x}{ϵ}),Q^{ϵ}(x)=Q(x,\frac{x}{ϵ}).$$

$$\begin{array}{r}{\mathcal{C}}_{11}^{ϵ}(x)\frac{\partial u_1^{ϵ}(t,x)}{\partial t}=\text{div}({\kappa }_1^{ϵ}(x)\nabla u_1^{ϵ}(t,x))+\frac{1}{{ϵ}^2}{Q}^{ϵ}(x)(u_2^{ϵ}(t,x)-u_1^{ϵ}(t,x))+q,\\ {\mathcal{C}}_{22}^{ϵ}(x)\frac{\partial u_2^{ϵ}(t,x)}{\partial t}=\text{div}({\kappa }_2^{ϵ}(x)\nabla u_2^{ϵ}(t,x))+\frac{1}{{ϵ}^2}{Q}^{ϵ}(x)(u_1^{ϵ}(t,x)-u_2^{ϵ}(t,x))+q,\end{array}$$

$$\begin{array}{r}{u}_1^{ϵ}(t,x)=u_{10}(t,x,\frac{x}{ϵ})+ϵu_{11}(t,x,\frac{x}{ϵ})+\cdots ,u_2^{ϵ}(t,x)=u_{20}(t,x,\frac{x}{ϵ})+ϵu_{21}(t,x,\frac{x}{ϵ})+\cdots ,\end{array}$$

$$\begin{array}{rl}{\mathcal{C}}_{11}& \frac{\partial (u_{10}+ϵu_{11}+\cdots )}{\partial t}\\ & =({\mathrm{div}}_x+\frac{1}{ϵ}{\mathrm{div}}_y)({\kappa }_1({\nabla }_x+\frac{1}{ϵ}{\nabla }_y)(u_{10}+ϵu_{11}+\cdots ))+\frac{1}{{ϵ}^2}Q(u_{20}+ϵu_{21}-u_{10}-ϵu_{11}+\cdots )+q,\\ {\mathcal{C}}_{22}& \frac{\partial (u_{20}+ϵu_{21}+\cdots )}{\partial t}\\ & =({\mathrm{div}}_x+\frac{1}{ϵ}{\mathrm{div}}_y)({\kappa }_2({\nabla }_x+\frac{1}{ϵ}{\nabla }_y)(u_{20}+ϵu_{21}+\cdots ))+\frac{1}{{ϵ}^2}Q(u_{10}+ϵu_{11}-u_{20}-ϵu_{21}+\cdots )+q,\end{array}$$

$$\begin{array}{rl}& {\mathrm{div}}_y({\kappa }_1(x,y){\nabla }_y{u}_{10}(t,x,y))+Q(x,y)(u_{20}(t,x,y)-u_{10}(t,x,y))=0\\ & {\mathrm{div}}_y({\kappa }_2(x,y){\nabla }_y{u}_{20}(t,x,y))+Q(x,y)(u_{10}(t,x,y)-u_{20}(t,x,y))=0.\end{array}$$

$$\begin{array}{r}-{\int }_Y{\kappa }_1{\nabla }_y{u}_{10}\cdot {\nabla }_y{u}_{10}\mathrm{d}y+{\int }_{Y}Q(u_{20}-u_{10})u_{10}\mathrm{d}y=0\\ -{\int }_Y{\kappa }_2{\nabla }_y{u}_{20}\cdot {\nabla }_y{u}_{20}\mathrm{d}y+{\int }_{Y}Q(u_{10}-u_{20})u_{20}\mathrm{d}y=0\end{array}$$

$$\begin{array}{r}{\int }_Y{\kappa }_1{\nabla }_y{u}_{10}\cdot {\nabla }_y{u}_{10}\mathrm{d}y+{\int }_Y{\kappa }_2{\nabla }_y{u}_{20}\cdot {\nabla }_y{u}_{20}\mathrm{d}y+{\int }_{Y}Q(u_{20}-u_{10}{)}^2\mathrm{d}y=0.\end{array}$$

$$\begin{array}{rl}& {\mathrm{div}}_x({\kappa }_1{\nabla }_y{u}_{10})+{\mathrm{div}}_y({\kappa }_1{\nabla }_x{u}_{10})+{\mathrm{div}}_y({\kappa }_1{\nabla }_y{u}_{11})+Q(u_{21}-u_{11})=0\\ & {\mathrm{div}}_x({\kappa }_2{\nabla }_y{u}_{20})+{\mathrm{div}}_y({\kappa }_2{\nabla }_x{u}_{20})+{\mathrm{div}}_y({\kappa }_2{\nabla }_y{u}_{21})+Q(u_{11}-u_{21})=0.\end{array}$$

$$\begin{array}{r}{\mathrm{div}}_y({\kappa }_1{\nabla }_y{u}_{11})+Q(u_{21}-u_{11})=-{\mathrm{div}}_y({\kappa }_1\nabla u_{10})\\ {\mathrm{div}}_y({\kappa }_2{\nabla }_y{u}_{21})+Q(u_{11}-u_{21})=-{\mathrm{div}}_y({\kappa }_2\nabla u_{20})\end{array}$$

$$\begin{array}{r}{\text{div}}_y({\kappa }_1(x,y)(e^i+{\nabla }_y{N}_1^i))+Q(x,y)(N_2^i-N_1^i)=0\\ {\text{div}}_y({\kappa }_2(x,y)(e^i+{\nabla }_y{N}_2^i))+Q(x,y)(N_1^i-N_2^i)=0\end{array}$$

$$\begin{array}{r}{\int }_Y{\mathcal{C}}_{11}\frac{\partial u_{10}}{\partial t}\mathrm{d}y={\int }_Y{\mathrm{div}}_x({\kappa }_1\nabla u_{10})\mathrm{d}y+{\int }_Y{\mathrm{div}}_x({\kappa }_1{\nabla }_y{u}_{11})\mathrm{d}y+{\int }_{Y}Q(u_{22}-u_{12})\mathrm{d}y+{\int }_{Y}q\mathrm{d}y\\ {\int }_Y{\mathcal{C}}_{22}\frac{\partial u_{20}}{\partial t}\mathrm{d}y={\int }_Y{\mathrm{div}}_x({\kappa }_2\nabla u_{20})\mathrm{d}y+{\int }_Y{\mathrm{div}}_x({\kappa }_2{\nabla }_y{u}_{21})\mathrm{d}y+{\int }_{Y}Q(u_{12}-u_{22})\mathrm{d}y+{\int }_{Y}q\mathrm{d}y\end{array}$$

$$\begin{array}{r}\left({\int }_Y{\mathcal{C}}_{11}\mathrm{d}y)+{\int }_Y{\mathcal{C}}_{22}\mathrm{d}y\right)\frac{\partial u_0}{\partial t}=\mathrm{div}({\kappa }_1^{\ast }\nabla u_0)+\mathrm{div}({\kappa }_2^{\ast }\nabla u_0)+{\int }_{Y}2q\mathrm{d}y\text{in}\Omega \end{array}$$

$$\begin{array}{r}{\kappa }_{1ij}^{\ast }(x)={\int }_Y{\kappa }_1(x,y)({\delta }_{ij}+\frac{\partial N_1^j(x,y)}{\partial y_i})dy,{\kappa }_{2ij}^{\ast }(x)={\int }_Y{\kappa }_2(x,y)({\delta }_{ij}+\frac{\partial N_2^j(x,y)}{\partial y_i})dy\end{array}$$

$$u_0(0,x)=\frac{⟨C_{11}⟩g_1(x)+⟨C_{22}⟩g_2(x)}{⟨C_{11}⟩+⟨C_{22}⟩}$$

$$\begin{array}{r}{\int }_Y{\kappa }_1(x,y){\nabla }_y{N}_1^i(x,y)\cdot {\nabla }_y{\varphi }_1(y)\mathrm{d}y-{\int }_{Y}Q(x,y)(N_2^i-N_1^i){\varphi }_1(y)\mathrm{d}y=-{\int }_Y{\kappa }_1(x,y)e^i\cdot {\nabla }_y{\varphi }_1(y)\mathrm{d}y\\ {\int }_Y{\kappa }_2(x,y){\nabla }_y{N}_2^i(x,y)\cdot {\nabla }_y{\varphi }_2(y)\mathrm{d}y-{\int }_{Y}Q(x,y)(N_1^i-N_2^i){\varphi }_2(y)\mathrm{d}y=-{\int }_Y{\kappa }_2(x,y)e^i\cdot {\nabla }_y{\varphi }_2(y)\mathrm{d}y\end{array}$$

$$|||({\varphi }_1,{\varphi }_2)|||=\Vert {\nabla }_y{\varphi }_1{\Vert }_{L^2(Y)}+\Vert {\nabla }_y{\varphi }_2{\Vert }_{L^2(Y)}+\Vert {\varphi }_1-{\varphi }_2{\Vert }_{L^2(Y)}.$$

$$B(x;({\varphi }_1,{\varphi }_2),({\psi }_1,{\psi }_2))$$

$$B(x;({\varphi }_1,{\varphi }_2),({\varphi }_1,{\varphi }_2))\ge c_1|||({\varphi }_1,{\varphi }_2)||{|}^2,\text{mbox}andB(x;({\varphi }_1,{\varphi }_2),({\psi }_1,{\psi }_2))\le c_2|||({\varphi }_1,{\varphi }_2)|||\cdot |||({\psi }_1,{\psi }_2)|||$$

$$\begin{array}{r}B(x;(N_1^i,N_2^i)({\varphi }_1,{\varphi }_2))=-{\int }_Y{\kappa }_1(x,y)e^i\cdot {\nabla }_y{\varphi }_1(y)\mathrm{d}y-{\int }_Y{\kappa }_2(x,y)e^i\cdot {\nabla }_y{\varphi }_2(y)\mathrm{d}y.\end{array}$$

$$\begin{array}{r}\Vert {\kappa }_1(x,\cdot )-{\kappa }_1(x^{\prime },\cdot ){\Vert }_{L^{\infty }(Y)}\le C|x-x^{\prime }|,\Vert {\kappa }_2(x,\cdot )-{\kappa }_2(x^{\prime },\cdot ){\Vert }_{L^{\infty }(Y)}\le C|x-x^{\prime }|,\\ \text{mbox}and\Vert Q(x,\cdot )-Q(x^{\prime },\cdot ){\Vert }_{L^{\infty }(Y)}\le C|x-x^{\prime }|\end{array}$$

$$\underset{\varphi \in {\mathcal{V}}_{L-l}}{\inf }\Vert {\nabla }_y(w-\varphi ){\Vert }_{L^2(Y)}+2^{L-l}\Vert w-\varphi {\Vert }_{L^2(Y)}\le C{2}^{-L+l}\Vert w{\Vert }_{H^2(Y)},$$

$$\begin{split}S_{l}={\cal T}_{l}\Big{\backslash}\bigg{(}{\bigcup_{k<l}S_{k}}\bigg{)}\end{split}$$

$$\begin{array}{r}B(x;({\bar{N}}_1^i,{\bar{N}}_2^i),({\varphi }_1,{\varphi }_2))=-{\int }_Y{\kappa }_1(x,y)e^i\cdot {\nabla }_y{\varphi }_1(y)\mathrm{d}y-{\int }_Y{\kappa }_2(x,y)e^i\cdot {\nabla }_y{\varphi }_2(y)\mathrm{d}y\end{array}$$

$$I_l^x(N_k^i)=\sum _{j=1}^n{c}_j{N}_k^i(x_j,\cdot ),$$

$$\begin{array}{rl}& B(x;({\overline{N_1^i}}^c,{\overline{N_2^i}}^c),({\varphi }_1,{\varphi }_2))\\ & =-\sum _{j=1}^n{c}_j{\int }_Y({\kappa }_1(x,y)-{\kappa }_1(x_j,y)){\nabla }_y\overline{N_1^i}(x_j,y)\cdot {\nabla }_y{\varphi }_1(y)\mathrm{d}y-\sum _{j=1}^n{c}_j{\int }_Y({\kappa }_1(x,y)-{\kappa }_1(x_j,y))e^i\cdot {\nabla }_y{\varphi }_1(y)\mathrm{d}y\\ & -\sum _{j=1}^n{c}_j{\int }_Y({\kappa }_2(x,y)-{\kappa }_2(x_j,y)){\nabla }_y\overline{N_2^i}(x_j,y)\cdot {\nabla }_y{\varphi }_2(y)\mathrm{d}y-\sum _{j=1}^n{c}_j{\int }_Y({\kappa }_2(x,y)-{\kappa }_2(x_j,y))e^i\cdot {\nabla }_y{\varphi }_2(y)\mathrm{d}y\\ & +\sum _{j=1}^n{c}_j{\int }_Y(Q(x_j,y)-Q(x,y))(\overline{N_1^i}(x_j,y)-\overline{N_2^i}(x_j,y))({\varphi }_1(y)-{\varphi }_2(y))\mathrm{d}y,\end{array}$$

$$\overline{N_k^i}(x,\cdot )={\overline{N_k^i}}^c(x,\cdot )+I_l^x(\overline{N_k^i}),$$

$$\begin{array}{r}{I}_l^x({\bar{N}}_k^i)={\bar{N}}_k^i(x^{\prime },\cdot ),k=1,2\end{array}$$

Full text

Key words. multiscale modelling, multi-continuum, homogenization, effective properties, hierarchical finite elements, optimal complexity.

1 Introduction

Media with multiple continua where each continuum interacts with other continua often entail multiple scales and high contrast. For example, fractured media can have multiscale and high contrast due to complex material properties and geography of fractures. Therefore, numerical simulations in this type of media can be expensive and require model reductions. This can be achieved by computing effective properties in each coarse block using the solutions of local representative volume element (RVE) problems. In multi-continuum approaches [2, 14, 12, 16, 13], equations for each continuum are written separately with so called interaction terms. Therefore, one has to deal with a system of equations and compute the effective properties from it.

There have been several methods to solve multiscale equations without computing effective properties and establishing homogenized equations. The multiscale finite element method (MsFEM) [10] solves local cell problems in coarse blocks with fine mesh to obtain basis functions that capture small-scale information. The generalized multiscale finite element method (GMsFEM) [7, 9, 6] follows the outline of MsFEM but adds some degrees of freedom in each coarse block by building the snapshot spaces and solving local spectral problems in the spaces. Using GMsFEM for multi-continuum models is considered in [8]. Although these methods have several advantages, they do not take into account the periodicity or local periodicity of the structures of the media of interest and sometimes introduce high computation cost when they have very small fine mesh sizes. In this paper, we establish an efficient algorithm for obtaining homogenized equations utilizing microscopic periodicity for multiscale multicontiuum systems with optimal computation complexity.

The multi-continuum system depends on small scales. We study the problem via the homogenization approach. For locally periodic multiscale problems such as those considered in this paper, the homogenized problem can be found by the two scale asymptotic expansion [3, 1, 11]. Constructing homogenized equations requires solving local cell problems in microscale representative domain. Computation cost for solving these cell problems can be very high since we have to solve different local cell problems at many macroscopic points.

In this work we use two-scale asymptotic expansion ([3], [1]) to derive the homogenized equations of a two scale multi-continuum system. Effective coefficients of the homogenized equation are established via the solutions of cell problems. Since the coefficients depend on both macro- and micro-scale variables, a different set of cell equations needs to be solved at each macroscopic point. The number of equations to be solved is thus very large. Solving them using the same small mesh at every macroscopic point is extremely expensive. We develop in this paper the hierarchical approach to solve the cell problems for the multi-continuum system for a large number of macroscopic points, using an optimal number of degrees of freedom, without sacrificing the accuracy. The method is developed for two scale elliptic problems in [4] (see also [5]). It solves cell problems for a dense hierarchical network of macroscopc points with different levels of resolution. The problems at those points belonging to a lower level in the hierarchy are solved with a higher level of accuracy. For the solution at a macroscopic point in a higher level in the hierarchy which are obtained with a lower level of accuracy, we use solutions at nearby macroscopic points that are solved with a higher level of accuracy to correct the error. We show that this hierarchical FE approach obtains the same level of accuracy at every macroscopic point as that obtained when every cell problem is solved with the highest level of resolution (we will refer to this as the full reference solve below), but uses only an essentially optimal number of degrees of freedom that is equal to that required to solve only one cell problem at the finest level of resolution (apart from a possible logarithmic factor).

The paper is organized as follows. In the next section, we set up the multiscale multi-continuum system; and we derive the homogenized equations from two scale asymptotic expansion. In Section 3, we outline the hierarchical finite element algorithm for solving the cell problems at a dense network of macroscopic points. We give a rigorous error estimates that show the algorithm has the equivalent accuracy as the full reference solve, at essentially optimal computation cost. In Section 4, we present numerical examples that verify the theoretical results. We compute the effective permeabilities using the hierarchical solve and the full solve. We find that the effective permeabilities obtained from these two approaches are essentially equal, with a very small relative error between each other. Finally, in Section 5, we rigorously prove the homogenization convergence for the two scale multi-continuum system. The paper ends with the conclusions in Section 6.

Throughout the paper, by $\nabla$, we denote the gradient with respect to $x$ of a function that depends only on the spatial variable $x$ and the temporal variable $t$, and by $\nabla_{x}$, we denote the partial gradient with respect to $x$ of a function that depends on $x$, $t$ and also other variables. Repeated indices indicate summation. The notation $\#$ denotes spaces of periodic functions.

2 Problem formulation

2.1 Homogenization of multi-continuum systems

In multi-continuum approaches, equations for each continuum are written separately. We denote by $u_{i}$ the solution for $i$th continuum. In the general case where each continuum interacts with every other continuum, we have the following system of equations introduced in [8]

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ is a domain ($d=2,3$), $\kappa_{i}^{\epsilon}$ are the multiscale permeability and $C_{ii}^{\epsilon}$ are the multiscale porosities, and the functions $Q_{i}^{\epsilon}$ of $(u_{1},...,u_{N})$ are exchange terms (see [2, 14, 12, 16, 13]) that describe the interaction of continua; $\epsilon$ represents the microscopic scale of the local variation. If each continuum only interacts with the background $u_{1}^{\epsilon}$, we have

[TABLE]

where $q$ is the source term. In this paper, we consider a two-continuum system. Let $Y$ be the unit cube in $\mathbb{R}^{d}$. Let ${\mathcal{C}}_{ii}(x,y)$, $\kappa_{i}(x,y)$ ($i=1,2$) be continuous functions on $\Omega\times Y$ which are $Y$-periodic with respect to $y$ and $q$ be a function in $L^{2}(\Omega)$. We assume further that there is a constant $c>0$ such that for all $x\in\Omega,y\in Y$

[TABLE]

We define the twoscale coefficients as

[TABLE]

We consider in this paper the case where the interaction terms are scaled as $O(1/\epsilon^{2})$; this case has the most interesting cell problems in the form of a coupled system. We consider the multiscale multi-continuum system

[TABLE]

with the Dirichlet boundary condition $u_{1}^{\epsilon}(t,x)=u_{2}^{\epsilon}(t,x)=0$ for $x\in\partial\Omega$, and with the initial condition $u_{1}^{\epsilon}(0,x)=g_{1}$, $u_{2}^{\epsilon}(0,x)=g_{2}$ where $g_{1}$ and $g_{2}$ are in $L^{2}(\Omega)$. We consider the following two scale asymptotic expansion of $u_{1}^{\epsilon}$ and $u_{2}^{\epsilon}$.

[TABLE]

where the functions $u_{1i}(t,x,y)$ and $u_{2i}(t,x,y)$ are periodic with respect to $y$. Performing the two scale asymptotic expansion, from (2.2) we obtain

[TABLE]

For the $O(\epsilon^{-2})$ terms, we obtain,

[TABLE]

From this, we have

[TABLE]

Adding these two equations, we obtain

[TABLE]

This implies $\nabla_{y}u_{10}=0$, $\nabla_{y}u_{20}=0$. i.e. $u_{10}$ and $u_{20}$ are independent of $y$, and $u_{10}(t,x)=u_{20}(t,x)$ as $Q(x,y)>c>0\ \forall\,x\in\Omega,y\in Y$. For the $O(\epsilon^{-1})$ terms in (2.3), we have,

[TABLE]

Since $u_{10}$ and $u_{20}$ are independent of $y$, we have

[TABLE]

Thus $u_{11}={\partial u_{10}\over\partial x_{i}}N^{i}_{1}$ and $u_{21}={\partial u_{20}\over\partial x_{i}}N^{i}_{2}$ where $N_{1}^{i}(x,\cdot)\in H^{1}_{\#}(Y)/\mathbb{R}$, and $N_{2}^{i}(x,\cdot)\in H^{1}_{\#}(Y)/\mathbb{R}$ are solutions of the cell problem

[TABLE]

where $e^{i}$ is the $i$th unit vector in the standard basis of ${\mathbb{R}}^{d}$. For the $O(\epsilon^{0})$ terms in (2.3), integrating over $Y$, one has

[TABLE]

Adding these two equations, one obtains the homogenized equation

[TABLE]

where $u_{0}=u_{10}=u_{20}$ and the $x$-dependent permeabilities are defined as

[TABLE]

We will show later that the homogenized matrices $\kappa^{*}_{1ij}(x)$ and $\kappa^{*}_{2ij}(x)$ are positive definite. We will also show that the initial condition for $u_{0}$ is

[TABLE]

where $\langle C_{ii}\rangle=\int_{Y}C_{ii}(y)dy$ for $i=1,2$. Equation (2.5) together with initial condtion (2.7) has a unique solution (see, e.g., [15]).

2.2 Uniqueness of solution to the cell problem

We write the system (2.4) in the variational form as

[TABLE]

where $\phi_{1},\phi_{2}\in H^{1}_{\#}(Y)$. Let $W$ be the space $H^{1}_{\#}(Y)\times H^{1}_{\#}(Y)/(c,c),c\in\mathbb{R}$. The space $W$ is equipped with the norm

[TABLE]

For $x\in\Omega$, we define the bilinear form $B(x;\cdot,\cdot):W\times W\to\mathbb{R}$ as

[TABLE]

for $(\phi_{1},\phi_{2})\in W$ and $(\psi_{1},\psi_{2})\in W$. From (2.1), we deduce that the bilinear form $B$ is uniformly coercive and bounded with respect to $x\in\Omega$, i.e. there are constants $c_{1}>0$ and $c_{2}>0$ such that

[TABLE]

for all $(\phi_{1},\phi_{2})\in W$ and $(\psi_{1},\psi_{2})\in W$. Adding the two equations in (2.8), we obtain

[TABLE]

Theorem 2.1**.**

Problem (2.8) has a unique solution $(N^{i}_{1},N^{i}_{2})\in W$.

Proof.

The conclusion follows from the boundedness and coerciveness of the bilinear form $B$ and the Lax-Milgram lemma. ∎

3 Hierarchical finite element algorithm

Computing effective coefficients $\kappa_{i}^{*}(x)$ requires the solutions of the cell problems (2.4) at many macroscopic points which can be very expensive. We develop in this section the hierarchical FE method which computes the solution of the cell problems at a dense network of macroscopic points using only an essentially optimal number of degrees of freedom which is equal to that for solving one cell problem (apart from a multiplying logarithmic factor). We assume that the coefficients are sufficiently smooth with respect to the macroscopic variable $x$. We make the following assumption.

Assumption 3.1**.**

There is a constant $C>0$ such that for all $x$, $x^{\prime}\in\Omega$,

[TABLE]

3.1 Overview of hierarchical algorithm

We develop an efficient hierarchical finite element algorithm to solve the coupled cell problem (2.4) numerically and to approximate the effective properties $\kappa_{i}^{*}(x)$ in (2.6) for a dense network of macroscopic points $x\in\Omega$. We follow the algorithm introduced in [4].

We outline the algorithm as follows.

Step 1 : Build nested finite element spaces. We employ Gelerkin FE to obtain an approximation of the solution $(N^{i}_{1},N^{i}_{2})\in W$ of (2.4) for each macroscopic point $x\in\Omega$ using FE spaces of different levels of resolution. We assume that there exists a hierarchy of FE spaces ${\cal V}_{0}\subset{\cal V}_{1}\subset\cdots\subset{\cal V}_{L}\subset H^{1}_{\#}(Y)$ where the integer index $L$ denotes the resolution level. We assume further the following approximation properties: for $w\in H^{2}_{\#}(Y)$,

[TABLE]

where the constant $C$ is independent of $L$ and $l$.

Step 2 : Build a hierarchy of macrogrids. We solve the cell equations at different macroscopic points $x\in\Omega$ with different levels of accuracy. We use the solutions solved with a higher accuracy level to correct the solutions obtained with a lower accuracy level. We achieve this by solving the cell problems at different macroscopic points using different FE spaces in the hierarchy in Step 1. This can be done by constructing a hierarchy of macro-grid points. We construct a nested macro-grid, ${\cal T}_{0}\subset{\cal T}_{1}\subset\cdots\subset{\cal T}_{L}\subset\Omega$ as follows. First, we build an initial grid ${\cal T}_{0}$ with a proper grid spacing $H$, the maximal distance between neighboring nodes. We then inductively construct ${\cal T}_{l}$, a refinement of ${\cal T}_{l-1}$, with grid spacing $H2^{-l}$. Then, we define the hierarchy of macro-grids, $\{S_{0},S_{1},\cdots,S_{L}\}$as $S_{0}={\cal T}_{0}$, $S_{1}={\cal T}_{1}\backslash S_{0}$, and for each $l>1$, we have

[TABLE]

We call the nodes in the lowest level grid $S_{0}$ the anchor points. In this way, we obtain a dense hierarchy of the macro-grids. That is, each point $x\in S_{l}$ has at least one point from one of the previous levels, $x^{\prime}\in\bigcup_{k<l}S_{k}$ such that dist($x$, $x^{\prime}$) $<$ $O(H2^{-l})$. Figures 2 and 2 show an example of 3-level hierarchy of macrogrids ${\cal T}_{l}$, $S_{l}$, $l=1,2,3$, constructed in $\Omega=[0,1]^{2}$.

Step 3 : Calculating the correction term. We relate the nested FE spaces and the hierarchy of macrogrids for our algorithm. We first solve the cell problems at anchor points using the standard Galerkin FE with FE space $\mathcal{V}_{L}$. That is, for the points in the coarsest macro-grid $S_{0}$, we solve the cell problems with the finest mesh. More precisely, we find $\bar{N^{i}_{1}}(x,\cdot),\bar{N^{i}_{2}}(x,\cdot)\in\mathcal{V}_{L}$, such that

[TABLE]

for all $\phi_{1},\phi_{2}\in\mathcal{V}_{L}$. Proceeding inductively, for $x\in S_{l}$ ($l=1,\cdots,L$), we choose the points $\{x_{1},x_{2},\cdots,x_{n}\}\in(\bigcup_{l^{\prime}<l}S_{l^{\prime}})$ so that the distance between $x$ and each point in $\{x_{1},x_{2},\cdots,x_{n}\}$ is $O(H2^{-l})$. This is possible from the assumption for the hierarchy of macroscopic points above. We define the $l$-th macro-grid interpolation by

[TABLE]

where the coefficients $c_{j}$ satisfy $\sum_{j=1}^{n}c_{j}=1$ ($k=1,2$). We refer to the $l$-th macro-grid interpolation of Galerkin approximations as $I_{l}^{x}(\bar{N^{i}_{k}})=\sum_{j=1}^{n}c_{j}\bar{N^{i}_{k}}(x_{j},\cdot)$. We solve the following problem: Find ${\bar{N^{i}_{1}}}^{c}(x,\cdot),{\bar{N^{i}_{2}}}^{c}(x,\cdot)\in{\cal V}_{L-l}$ such as

[TABLE]

for all $\phi_{1}$, $\phi_{2}$ $\in\mathcal{V}_{L-l}$. Note that right-hand side data is all known since we have already computed $\{\bar{N^{i}_{k}}(x_{j},\cdot)\}_{j=1}^{n}$ inductively using finer mesh spaces at macro-grid points in $(\bigcup_{l^{\prime}<l}S_{l^{\prime}})$. We let

[TABLE]

be the FE approximation for $N^{i}_{k}(x,\cdot)$ where $k=1,2$. A main goal of this paper is to prove that the approximation (3.3) for $N^{i}_{k}(x,\cdot)$ has the same order of accuracy compared to the approximation we obtain by solving (2.8) using the finest FE space $\mathcal{V}_{L}$ at all macroscopic points. We also prove that we reduce the computation cost with the approximation (3.3) to the optimal level.

Remark. In the following, for simplicity, we use a simple 1-point interpolation to compute the correction term $(\bar{N_{1}^{i}}^{c},\bar{N_{2}^{i}}^{c})$. More precisely, for $x\in S_{l}$ we choose $x^{\prime}\in(\bigcup_{l^{\prime}<l}S_{l^{\prime}})$ such that $\mathrm{dist}(x,x^{\prime})<O(H2^{-l})$. We let

[TABLE]

be the macro-grid interpolation. The FE approximation is

[TABLE]

Remark. Note that as the level $l$ goes higher, we use coarser FE spaces for the corresponding finer macro grids. This balance guarantees that although we use coarser FE spaces, the FE error is still optimal, but with much less computation cost.

3.2 Error estimates

We require that the coefficients $\kappa_{i}$ and $Q$ satisfy Assumption 3.1 and (2.1). We prove that the hierarchical method achieves the same order of accuracy as the full solve. For simplicity, we consider 1-point interpolation for our proof; the proof for the general case is similar.

Lemma 3.1**.**

There exists a positive number $C$ such that $|||(N_{1}^{i}(x,\cdot),N_{2}^{i}(x,\cdot))|||\leq C$ for all $x\in\Omega$.

Proof.

From (2.8), we obtain

[TABLE]

Using the uniform coercivity of the bilinear form $B(x;\cdot,\cdot)$ with respect to $x$, we get

[TABLE]

for $C>0$. From this we get the conclusion. ∎

Let ${N^{i}_{k}}^{c}(x,\cdot)={N^{i}_{k}}(x,\cdot)-{N^{i}_{k}}(x^{\prime},\cdot)$. We have that $({N^{i}_{1}}^{c}(x,\cdot),{N^{i}_{2}}^{c}(x,\cdot))\in W$ satisfies

[TABLE]

$\forall\,(\phi_{1},\phi_{2})\in W$.

Proposition 3.2**.**

There exists $C>0$ such that

[TABLE]

for $x\in{\mathcal{T}}_{L}$.

Proof.

From (3.4), for $(\phi_{1},\phi_{2})=({N_{1}^{i}}^{c}(x,\cdot),{N_{2}^{i}}^{c}(x,\cdot))\in W$ we have

[TABLE]

As $\nabla_{y}N^{i}_{1}(x^{\prime},\cdot)$ and $\nabla_{y}N^{i}_{2}(x^{\prime},\cdot)$ are uniformly bounded in $L^{2}(Y)$ with respect to $x\in\Omega$ by Lemma 3.1. From Assumption 3.1we have

[TABLE]

Thus

[TABLE]

where the constant $C$ is independent of $x$. ∎

Lemma 3.3**.**

There is a positive constant $C$ such that $||\Delta_{y}N^{i}_{1}(x,\cdot)||_{L^{2}(Y)}+||\Delta_{y}N^{i}_{2}(x,\cdot)||_{L^{2}(Y)}\leq C$ for all $x\in\Omega$.

Proof.

We rewrite cell problem (2.4) as

[TABLE]

[TABLE]

Rearranging these equations, we have,

[TABLE]

[TABLE]

By the uniform boundedness of $|||(N_{1}^{i}(x,\cdot),N_{2}^{i}(x,\cdot))|||$ with respect to $x$ and Lemma 3.1, we deduce that $||\Delta_{y}N^{i}_{1}(x,\cdot)||_{L^{2}(Y)}$ and $||\Delta_{y}N^{i}_{2}(x,\cdot)||_{L^{2}(Y)}$ are uniformly bounded for all $x\in\Omega$. ∎

Lemma 3.4**.**

There exists a positive constant $C$ such that

[TABLE]

for all $x\in{\mathcal{T}}_{L}$.

Proof.

From (3.4), we have

[TABLE]

[TABLE]

Therefore,

[TABLE]

[TABLE]

From Lemma 3.1 and Proposition 3.2, we have

[TABLE]

for some constant $C>0$. ∎

We choose $({N_{1}^{i}}^{c},{N_{2}^{i}}^{c})\in W$ such that

[TABLE]

We then have

Lemma 3.5**.**

There is a positive constant $C$ such that $||{N^{i}_{1}}^{c}(x,\cdot)||_{L^{2}(Y)}\leq C|x-x^{\prime}|$ and $||{N^{i}_{2}}^{c}(x,\cdot)||_{L^{2}(Y)}\leq C|x-x^{\prime}|$ for all $x\in{\mathcal{T}}_{L}$.

Proof.

We note that

[TABLE]

Since $\int_{Y}({N^{i}_{1}}^{c}+{N^{i}_{2}}^{c})\mathrm{d}y=0$, by Poincare inequality, and $(\ref{eq:main32})$, the following inequalities hold.

[TABLE]

And then by (3.6),

[TABLE]

∎

Proposition 3.6**.**

There is a constant $C>0$ such that $||{N^{i}_{1}}^{c}||_{H^{2}(Y)}\leq C|x-x^{\prime}|$ and $||{N^{i}_{2}}^{c}||_{H^{2}(Y)}\leq C|x-x^{\prime}|$ for all $x\in{\mathcal{T}}_{L}$.

Proof.

Let $\omega\subset\mathbb{R}^{d}$ be a domain such that $Y\subset\omega$. Let $\phi\in{\cal C}^{\infty}_{0}(\omega)$ be such that $\phi=1$ in $Y$. We have

[TABLE]

Since $\phi{N^{i}_{1}}^{c}=0$ on $\partial\omega$, applying elliptic regularity, we have

[TABLE]

By Proposition 3.2, Lemmas 3.4 and 3.5, and the $Y$-periodicity of ${N^{i}_{1}}^{c}$,

[TABLE]

for all $x\in{\mathcal{T}}_{L}$. Then from (3.7), $||{N^{i}_{1}}^{c}||_{H^{2}(Y)}\leq C|x-x^{\prime}|.$ Similarly, $||{N^{i}_{2}}^{c}||_{H^{2}(Y)}\leq C|x-x^{\prime}|$ for $C>0$. ∎

We consider the problem: Find ${\bar{\bar{N}}^{i}_{1}}^{c}(x,y)\in{\cal V}_{L-l}$ and ${\bar{\bar{N}}^{i}_{2}}^{c}(x,y)\in{\cal V}_{L-l}$ such that

[TABLE]

for all $\phi_{1}\in{\cal V}_{L-l}$ and $\phi_{2}\in{\cal V}_{L-l}$. This is the FE approximation of (3.4). We then have the following result.

Lemma 3.7**.**

There is a positive constant $C^{0}$ such that

[TABLE]

Proof.

It follows from Cea’s Lemma, Proposition 3.6 and (3.1) that

[TABLE]

∎

Proposition 3.8**.**

There is a constant $c_{l}>0$ which only depends on the level $S_{l}$ of $x\in{\cal T}_{L}$ such that

[TABLE]

Proof.

We will prove the proposition by induction. The conclusion holds for $l=0$. We assume that for all $x^{\prime}\in S_{l^{\prime}}$ where $l^{\prime}\leq l-1$.

[TABLE]

From (3.2) and (3.8), we have

[TABLE]

for all $\phi_{1}\in{\cal V}_{L-l}$ and $\phi_{2}\in{\cal V}_{L-l}$. From Assumption 3.1 and the induction hypothesis, we have

[TABLE]

where $\gamma>0$ is independent of $x$ and $l$. By Lemma 3.7 and (3.10),

[TABLE]

Using $\bar{N}^{i}_{k}(x,y)=\bar{N^{i}_{k}}^{c}(x,y)+\bar{N}^{i}_{k}(x^{\prime},y)$, We have

[TABLE]

where

[TABLE]

∎

Theorem 3.9**.**

Under Assumption 3.1 and the uniform boundedness of $\kappa_{i}(x,y)$ and $Q(x,y)$, there is a positive constant $C_{*}$ which depends only on the functions $\kappa_{1}$, $\kappa_{2}$ and $Q$ so that,

[TABLE]

for $x\in S_{l}$.

Proof.

We let $\bar{l}$ be an integer independent of $L$ such that $l2^{-l}<\frac{1}{2\gamma}$ for $l>\bar{l}$. And let

[TABLE]

where $C^{0}$ and $c_{l}$ are the constants in Lemma 3.7 and Proposition 3.8. Now we prove

[TABLE]

by induction. From (3.14), this holds for all $l\leq\bar{l}$. Suppose that (3.15) holds for all $l^{\prime}\leq l$. Then from (3.12), we obtain

[TABLE]

∎

Theorem 3.10**.**

The total number of degrees of freedom required to solve (2.8) for all points in $S_{0},S_{1},\cdots,S_{L}$ is $\mathcal{O}((L+1)2^{dL})$ for the hierarchical solve while it is $\mathcal{O}((2^{dL})^{2})$ in the full solve where cell problems are solved with the finest mesh level at all macrogrid points.

Proof.

Since the number of macroscopic points in $S_{l}$ is $\mathcal{O}(2^{dl})$, and the space $\mathcal{V}_{L-l}$ is of dimension $\mathcal{O}(2^{d(L-l)})$, the total number of degrees of freedom for solving (2.8) for all points in $S_{l}$ is $\mathcal{O}(2^{dl})\mathcal{O}(2^{d(L-l)})=\mathcal{O}(2^{dL})$. Therefore, the total number of degrees of freedom required to solve (2.8) for all points in $S_{0},S_{1},\cdots,S_{L}$ is $\mathcal{O}((L+1)2^{dL})$. ∎

4 Numerical example

In this section, we apply the hierarchical finite element algorithm to a numerical example for computing the effective coefficients of a multiscale multi-continuum system at a dense network of macrogrid points. To show the accuracy of the algorithm, we compare the results to the approximations to the effective coefficients obtained from solving the cell problems using the finest meshes at all macroscopic points.

4.1 Numerical Implementation

We let $\Omega=[0,1]^{2}$ be the macroscopic domain and $Y=[0,1]^{2}$ be the unit cell. We consider the locally periodic coefficients

[TABLE]

where the constant $a$ is chosen below. We use 4 square meshes in $[0,1]^{2}$ to construct a nested sequence of FE spaces, $\{{{\cal V}_{3-l}}\}^{3}_{l=0}$ so that the mesh size of each space is $h_{l}=2^{l}(2^{-4})$ for $l=0,1,2,3$. Since $\kappa_{1}$, $\kappa_{2}$ and $Q$ are independent of $x_{2}$, we only consider 1-dimensional macrogrids in $[0,1]$. The nested macrogrids $\{{\cal T}_{l}\}^{L}_{l=0}\subset[0,1]$ and the subsequent macrogrid hierarchy, $\{S_{l}\}^{3}_{l=0}$ are constructed as follows. We first let ${\cal T}_{0}=S_{0}=\{0,\frac{1}{2},1\}$. Considering that our macrogrids have grid spacing $H2^{-l}$ for $l=0,1,2,3$, where $H=\frac{1}{2}$ in this case, we have following hierarchy of macrogrids.

[TABLE]

Figure 3 indicates how these macrogrids and the approximation spaces are related in numerical implementation.

We implement the algorithm as follows. For $x^{\prime}\in S_{0}=\{0,\frac{1}{2},1\}$, we solve (2.8) for $\bar{N^{i}_{1}}(x^{\prime},\cdot)$, $\bar{N^{i}_{2}}(x^{\prime},\cdot)\in{\cal V}_{3}$, for all $\phi_{1},\ \phi_{2}\in{\cal V}_{3}$ by the standard Galerkin FEM. We then use a simple 1-point interpolation to compute the correction terms. That is, for $x\in S_{l}$ we choose $x^{\prime}\in(\bigcup_{k<l}S_{k})$ such that $|x^{\prime}-x|\leq 2^{-l}$. We let the $l$th macrogrid interpolation be

[TABLE]

We find $\bar{N^{i}_{1}}^{c}(x,y)$ and $\bar{N^{i}_{2}}^{c}(x,y)$ in $\mathcal{V}_{L-l}$ such that

[TABLE]

and

[TABLE]

for $\forall\phi_{1},\phi_{2}\in\mathcal{V}_{L-l}$. We let

[TABLE]

be the approximation to $N_{k}^{i}(x,\cdot)$. We continue inductively. For example, for $x=\frac{1}{2}\in S_{0}$, we compute $\bar{N^{i}_{1}}(\frac{1}{2},\cdot)$, $\bar{N^{i}_{2}}(\frac{1}{2},\cdot)$ using the standard Galerkin FEM. Then for $\frac{3}{8}\in S_{1}$, we find the correction terms $\bar{N^{i}_{1}}^{c}(\frac{3}{8},\cdot)$, $\bar{N^{i}_{2}}^{c}(\frac{3}{8},\cdot)\in{\cal V}_{L-1}$ that satisfy (4.1) and (4.2), where $x^{\prime}=\frac{1}{2}$. And we let the solutions at $x=\frac{3}{8}$ be

[TABLE]

We continue this procedure based on Figure 3.

Tables 1 and 2 indicate $\kappa^{*}_{111}$ and $\kappa^{*}_{211}$ obtained by both the hierarchical solve and the full solve where the finest mesh is used for all cell problems, at each $x_{1}$ and the relative errors between them, where relative errors are calculated by $\frac{100|\kappa^{*}_{full}-\kappa^{*}_{hier}|}{\kappa^{*}_{full}}$ with obvious notations for $a=1$ and $a=0.1$ respectively. The results show clearly that the effective coefficients obtained from hierarchical algorithm are very closed to the reference effective coefficients. We can see from the tables that relatively large errors occur at the highest level macroscopic points where more than one layer of corrections is performed, i.e. the corrector itself is corrected by the solution at a macroscopic point belonging to a lower level. We note that the error for the case $a=0.1$ is much smaller as the change of $\kappa_{i}$ in $x$ is much smaller. That is, large Lipschitz constants in Assumption 3.1 tend to result in large errors. The results in Tables 1 and 2 are obtained when only one corrector point is employed. If we use more corrector points, the error can be reduced significantly. In Table 3 we show the relative errors, in comparison to the coefficients obtained from the full solve where the finest mesh is used for all the cell problems, for the effective coefficients obtained from the hierarchical solve for the two cases where one point and two point interpolations are used. The table shows that the result can be improved by employing two point interpolation.

5 Proof of homogenization convergence

In this section, we prove rigorously the homogenization convergence, i.e. the convergence of the solution of the two scale equation (2.2) to the solution of the homogenized equation (2.5). Throughout this section, we denote the spaces $L^{2}(\Omega)$ and $H_{0}^{1}(\Omega)$ as $H$ and $V$ respectively. We recall the two-scale multi-continuum system

[TABLE]

[TABLE]

We have the following theorem.

Lemma 5.1**.**

The solution $(u_{1}^{\epsilon}$, $u_{2}^{\epsilon})$ of (5.1) and (5.2) are uniformly bounded in $L^{\infty}(0,T;H)$ and $L^{2}(0,T;V)$.

Proof.

Multiplying $\phi_{1}$ and $\phi_{2}$ $\in V$ to (5.1) and (5.2) respectively and integrating over $\Omega$, one has

[TABLE]

Summing these equations, we get

[TABLE]

$\forall\phi_{1},\phi_{2}\in V$. Substituting $u_{1}^{\epsilon}$ and $u_{2}^{\epsilon}$ into $\phi_{1}$ and $\phi_{2}$ in (5.4) respectively, we have

[TABLE]

Integrating this equation over $(0,\tau)$, we get

[TABLE]

Therefore,

[TABLE]

Using the uniform boundedness from below of $C_{11}^{\epsilon}$ and $C_{22}^{\epsilon}$, we have

[TABLE]

Choosing $\delta$ sufficiently small, we deduce that $u_{1}^{\epsilon}$ and $u_{2}^{\epsilon}$ are uniformly bounded in $L^{\infty}(0,T;H)$ and $L^{2}(0,T;V)$. ∎

Note that because of the 5th term of equation (5.5), $\displaystyle\lim_{\epsilon\to 0}u_{1}^{\epsilon}=\lim_{\epsilon\to 0}u_{2}^{\epsilon}$. Thus, there exist subsequences of $u_{1}^{\epsilon}$ and $u_{2}^{\epsilon}$, which we still denote by $u_{1}^{\epsilon}$ and $u_{2}^{\epsilon}$ , and $u_{0}$ such that

[TABLE]

Recall that $({N^{i}_{1}},{N^{i}_{2}})\in W$ is the solution of cell problem.

[TABLE]

We assume that $N_{1}^{i}$ and $N_{2}^{i}$ are sufficiently smooth with respect to both $x$ and $y$. Let $\omega_{1}(x)=\frac{x_{i}}{\epsilon}+N^{i}_{1}(x,\frac{x}{\epsilon})$ and $\omega_{2}(x)=\frac{x_{i}}{\epsilon}+N^{i}_{2}(x,\frac{x}{\epsilon})$. We define $\omega_{1}^{\epsilon}$ and $\omega_{2}^{\epsilon}$ as

[TABLE]

Assuming that $\kappa_{1}$, $\kappa_{2}$, $N_{1}^{i}$ and $N_{2}^{i}$ are sufficiently smooth, for all $\psi_{1},\psi_{2}\in V$ we have

[TABLE]

and

[TABLE]

due to (5.6). Let $\phi_{1}(x)=\phi(x)\omega_{1}^{\epsilon}(x)$, $\phi_{2}(x)=\phi(x)\omega_{2}^{\epsilon}(x)$ where $\phi\in\mathcal{C}^{\infty}_{0}(\Omega)$ in (5.3), we have

[TABLE]

Let $\psi_{1}(x)$ and $\psi_{2}(x)$ in (5.7) and (5.8) be $\phi u_{1}^{\epsilon}$ and $\phi u_{2}^{\epsilon}$ respectively. We have

[TABLE]

Let $\psi\in C^{\infty}_{0}(0,T)$. We multiply (5.9) and (5.10) by $\psi$ and intergrate over $(0,T)$ with respect to $t$. After subtracting the resulting equations by each other, we obtain

[TABLE]

We have the following lemma.

Lemma 5.2**.**

The functions $\int_{0}^{T}\psi(t)u_{1}^{\epsilon}(x,t)dt$ and $\int_{0}^{T}\psi(t)u_{2}^{\epsilon}(x,t)dt$ converge strongly in $H$ to $\int_{0}^{T}\psi(t)u_{0}(x,t)dt$.

*Proof  * This is the standard result in Jikov et al. [11]. As $u_{1}^{\epsilon}$ is uniformly bounded in $L^{2}(0,T;V)$, $\int_{0}^{T}\psi(t)u_{1}^{\epsilon}(x,t)\mathrm{d}t$ is uniformly bounded in $V$ when $\epsilon\to 0$. Thus we can extract a subsequence which converges weakly in $V$ and strongly in $H$. As for all $\phi\in C^{\infty}_{0}(\Omega)$,

[TABLE]

the limit is $\int_{0}^{T}\psi(t)u_{0}(x,t)\mathrm{d}t$. $\Box$

We have

[TABLE]

As $C_{11}^{\epsilon}$ converges weakly to $\int_{Y}C_{11}(x,y)\mathrm{d}y$ in $H$, $\int_{0}^{T}u_{1}^{\epsilon}{\partial\psi\over\partial t}\mathrm{d}t$ converges weakly to $\int_{0}^{T}u_{0}{\partial\psi\over\partial t}\mathrm{d}t$ in $V$, we have

[TABLE]

We note that

[TABLE]

Also, note that due to $Y$-periodicity of $\kappa$ and $N^{i}$, we have

[TABLE]

We observe that $x_{i}+\epsilon N^{i}_{1}\rightarrow x_{i}$ strongly in $H$. Passing to the limit in the left hand side of (5.11), we obtain from Lemma 5.2,

[TABLE]

Note that $\int_{\Omega}q\phi\omega_{k}^{\epsilon}\mathrm{d}x\rightarrow\int_{\Omega}q\phi x_{i}\mathrm{d}x$ since $\omega_{k}^{\epsilon}\phi\rightarrow x_{i}\phi$ in $H$ weakly. Thus from (5.11) and (5.12), we have

[TABLE]

Let $\phi_{1}$ and $\phi_{2}$ in (5.3) be $\phi x_{i}$ where $\phi\in C^{\infty}_{0}(\Omega)$. Adding the two equations, we have

[TABLE]

Passing to the limit, we obtain

[TABLE]

Using (5.13) and (5.14), one obtains

[TABLE]

Since $\kappa_{1}$, $\kappa_{2}$, $N_{1}^{i}$ and $N_{2}^{i}$ are independent of $t$, by Lemma 5.2, we have

[TABLE]

Therefore, we have

[TABLE]

From this, we deduce

[TABLE]

For consistency with formula (2.6), we note the following result.

Lemma 5.3**.**

$\int_{Y}\kappa_{1}{\partial N^{j}_{1}(x,y)\over\partial y_{i}}\mathrm{d}y+\int_{Y}\kappa_{2}{\partial N^{j}_{2}(x,y)\over\partial y_{i}}\mathrm{d}y=\int_{Y}\kappa_{1}{\partial N^{i}_{1}(x,y)\over\partial y_{j}}\mathrm{d}y+\int_{Y}\kappa_{2}{\partial N^{i}_{2}(x,y)\over\partial y_{j}}\mathrm{d}y$**

Proof.

From the cell problem, we have

[TABLE]

Thus,

[TABLE]

Now we show

[TABLE]

From the cell problem, we know that

[TABLE]

Thus, by (5.16) and (5.17), we have the result. ∎

Theorem 5.4**.**

Assume that the solution $N_{1}^{i}$ and $N_{2}^{i}$ of cell problem (2.8) belong to $C^{2}(\bar{\Omega},C^{2}(\bar{Y}))$ and the coefficients $\kappa_{1}$ and $\kappa_{2}$ belong to $C^{1}(\bar{\Omega},C^{1}(\bar{Y}))$. The limit function $u_{0}$ of the sequences $u_{1}^{\epsilon}$, $u_{2}^{\epsilon}$ is the unique solution of the homogenized equation (2.5) with the initial condition (2.7).

Proof.

Note that from the equation (5.1), we obtain

[TABLE]

for all $\phi\in\mathcal{C}^{\infty}_{0}(\Omega)$ and $\psi\in C^{\infty}_{0}((0,T))$. Passing to the limit, from (5.15), Lemmas 5.2 and 5.3, we have

[TABLE]

where

[TABLE]

We now show the initial condition. Adding (5.1) and (5.2), we have

[TABLE]

As $u_{1}^{\epsilon}$ and $u_{2}^{\epsilon}$ are bounded in $L^{2}(0,T;V)$, we deduce that ${\mathcal{C}}_{11}^{\epsilon}{\partial u_{1}^{\epsilon}\over\partial t}+{\mathcal{C}}_{22}^{\epsilon}{\partial u_{2}^{\epsilon}\over\partial t}$ is bounded in $L^{2}(0,T;V^{\prime})$. Let $\psi(t,x)\in{\mathcal{C}}_{0}^{\infty}(0,T;V)$, i.e. $\psi(0,x)=\psi(T,x)=0$. We have

[TABLE]

This shows that the weak limit of ${\mathcal{C}}_{11}^{\epsilon}{\partial u_{1}^{\epsilon}\over\partial t}+{\mathcal{C}}_{22}^{\epsilon}{\partial u_{2}^{\epsilon}\over\partial t}$ in $L^{2}(0,T;V^{\prime})$ is \big{(}\langle{\mathcal{C}}_{11}\rangle+\langle{\mathcal{C}}_{22}\rangle\big{)}\frac{\partial u_{0}}{\partial t}. Now we choose $\psi\in{\mathcal{C}}^{\infty}(0,T;V)$ so that $\psi(T,x)=0$. Then

[TABLE]

On the other hand

[TABLE]

This shows that \big{(}\langle{\mathcal{C}}_{11}\rangle+\langle{\mathcal{C}}_{22}\rangle\big{)}u_{0}(0,x)=\langle{\mathcal{C}}_{11}\rangle g_{1}(x)+\langle{\mathcal{C}}_{22}\rangle g_{2}(x). i.e. the initial condition of $u_{0}$ is

[TABLE]

∎

6 Conclusions

In this paper, we developed an efficient algorithm for computing the effective coefficients of a coupled multiscale multi-continuum system. We derived the coupled cell problems and the homogenized equation from twosale asymptotic expansion. We solved the cell problems using hierarchical FE algorithm and use the solutions to compute the effective coefficients. To establish the hierarchical FE algorithm, we first constructed a dense hierarchy of macrogrids and the corresponding nested FE spaces. Based on the hierarchy, we solve the cell problems using different resolution FE spaces at different macroscopic points. We use solutions solved with a higher level of accuracy to correct solutions obtained with a lower level of accuracy at nearby macroscopic points. We rigorously showed that this hierarchical FE method achieves the same order of accuracy as the reference full solve where cell problems at every macroscopic point are solved with the highest level of accuracy, at a significantly reduced computational cost, using an essentially optimal number of degrees of freedom. For numerical example, we applied this algorithm to a multi-continuum model in a two dimensional domain. The algorithm was implemented on macroscopic points in a one dimensional domain. The numerical results strongly support the error estimates we provided in section 3.

Acknowledgment A part of this work is conducted when Jun Sur Richard Park was a visiting PhD student at Nanyang Technological University (NTU) under East Asia and Pacific Summer Institutes (EAPSI) programme organized by the US National Science Foundation (NSF) and Singapore National Research Foundation (NRF) under Grant No. 1713805. Jun Sur Richard Park thanks US NSF and Singapore NRF for the financial support and NTU for hospitality. Viet Ha Hoang is supported by the MOE AcRF Tier 1 grant RG30/16 and the MOE Tier 2 grant MOE2017-T2-2-144.