# Adaptive Double-Diffusion Model and Comparison to a Highly Heterogeneous Micro-Model

Authors: Viviane Klein; Malgorzata Peszynska
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Double-diffusion model is used to simulate slightly compressible fluid flow in periodic porous media as a macro-model in place of the original highly heterogeneous micro-model. In this paper, we formulate an adaptive two-grid numerical finite element discretization of the double-diffusion system and perform a comparison between the micro- and macro-model. Our numerical results show that the micro-model solutions appear to converge to the macro-model linearly with the parameterε of periodic geometry. For the two-grid discretization, the a priori and a posteriori error estimates are proved, and we show how to adapt the grid for each component independently.


---

## Body

## 1. Introduction

When modeling phenomena in highly heterogeneous media, one frequently finds that the coefficients of differential equations describing these phenomena vary by several orders of magnitude between close-by locations. Consequently, the solutions to the exact, that is,micro-models, vary at multiple separate spatial and time scales, and it is convenient to work with their spatial and temporal averages, that is, with the solutions to macro-models. Various multiscale modeling techniques have been introduced with the aim to derive, analyze, and approximate micro-models by the macro-models. In particular, the mathematical framework for construction and analysis of the average, that is, homogenized, models for multiscale media is well understood, see, for example, [1–3]. However, there are few results devoted to the comparison between the micro- and macro-models for evolution equations and to their adaptive numerical discretizations.Consider the following transient diffusion model which describes the densityρ of single-phase slightly compressible flow in a porous domain Ω:
(1.1)ϕ∂ρ∂t-∇⋅(k∇ρ)=f,inΩ,t>0,
where the coefficients ϕ, k denote, respectively, the porosity and conductivity.In this paper, we propose a two-grid adaptive finite element approximation of the following multiscale averaged model for (1.1) in highly heterogeneous media such as porous media with fractures(1.2a)ϕ~1∂ρ~1∂t-∇⋅(k~1∇ρ~1)+c(ρ~1-ρ~2)=f~,inΩ,t>0,(1.2b)ϕ~2∂ρ~2∂t-∇⋅(k~2∇ρ~2)+c(ρ~2-ρ~1)=0,inΩ,t>0,proposed in [4, 5]. The system (1.2a) and (1.2b) describes the average behavior of fluids in media with two or more subregions of distinct features, where the subregions may or may not be connected globally at micro-scale. At macro-scale the interaction of fluids associated with these subregions is modeled with the exchange termc(ρ~2-ρ~1), and the coefficients k~i, ϕ~i are computed from micro-scale geometry and coefficients. The system (1.2a) and (1.2b) is useful in applications in porous media as well as in description of other multiscale evolution phenomena in heterogeneous media, for example, heat conduction. In addition, its special, degenerate, case of (1.2a) and (1.2b) with k~2=0 known as the Warren-Root model [6](1.3a)ϕ~1∂ρ~1∂t-∇⋅(k~1∇ρ~1)+c(ρ~1-ρ~2)=f~,inΩ,t>0,(1.3b)ϕ~2∂ρ~2∂t+c(ρ~2-ρ~1)=0,inΩ,t>0,is a popular way of modeling of subscale diffusion accompanying advection and adsorption, see [7–10].The macro-model (1.2a) and (1.2b) of interest to this paper was introduced in [4] based on heuristic arguments. The model, while intuitively clear and useful for practical purposes, for a long time was lacking a rigorous multiscale analysis and interpretation via an associated appropriate micro-model. Only recently in [11] it was shown rigorously by two-scale convergence that (1.2a) and (1.2b) is indeed a homogenized limit of (1.1) in certain highly heterogeneous media. Most interesting is that the original microstructure which leads to (1.2a) and (1.2b) is composed of three rather than two media as in [12]. Also, the derivation in [11] demonstrates clearly that in two-dimensional geometries and more generally for quasistatic models in disconnected microscale subregions we have k~2=0.However, the analysis in [11] does not include a direct comparison of (1.1) and (1.2a) and (1.2b), and this motivated our research. We compare the computational solutions to (1.2a) and (1.2b) and (1.1) quantitatively in function of the microstructure parameters. In addition, since typically 0≤k~2≪k~1, it is natural to suggest that in numerical simulation ρ~1 and ρ~2 should be approximated on separate grids, since they evolve differently and since this can lead to more efficient numerical algorithms. In fact, the spatial grids should be chosen adaptively guided by some a posteriori estimators and should be allowed to vary in time. The idea of using two grids extends that in [13] where we considered a stationary analogue of (1.2a) and (1.2b) and proved upper and lower bounds for residual type a posteriori estimators. The estimators proposed in [13] were shown to be robust, that is, insensitive to k~2≈0 and to other parameters. In this paper we extend those to the evolution system (1.2a) and (1.2b) and present an a posteriori error estimator of residual type based partly on the work in [14–16] for scalar equations.The aim of this paper is thus twofold, and the paper is organized as follows. In Section2 we introduce the notation and the generic flow model as well as its numerical approximation. In Section 3 we develop the details of micro- and macro- models and compare their discrete solutions in function of a characteristic parameter of the microstructure ε. In Section 4, we define the numerical approximations to these models and formulate the a-priori and a-posteriori estimates for the approximation error. We then illustrate the use of a-posteriori estimators with an adaptive example. The appendix provides the technical proofs of some of the results derived in this paper.


## 2. Notation and Preliminaries on the Flow Model

Here we introduce the notation and preliminaries. We follow, for example, [17, 18]. Let the domain of flow Ω⊂ℝd, d=2,3, be a bounded polygonal region with piecewise smooth boundary ∂Ω=ΓD∪ΓN¯, where ΓD, ΓN are disjoint regions of the boundary ∂Ω where Dirichlet and Neumann boundary conditions are prescribed, respectively.For any subsetω⊆Ω, we use the standard notation for Lebesgue L2(ω) and Sobolev spaces Hk(ω), k∈ℕ. These are equipped with the usual seminorms |·|k,ω, norms ∥·∥k,ω:=∥·∥Hk(ω) and the usual scalar product (s) (f,g)ω:=(f,g)L2(ω). We denote ∥f∥ω:=(f,f)ω1/2. If ω=Ω, the subscript ω will be omitted. We also let HD1(Ω):={v∈H1(Ω):v|ΓD=0} under the standard norm in H1(Ω), ∥v∥1:=(∥v∥2+∥∇v∥2)1/2.For anyT>0 and any vector space S, the space L2(0,T;S) consists of all square-integrable functions with values in S such that ∥v∥L2(0,T;S):=(∫0T‍∥v∥S2)1/2 is finite. The space C([0,T];S) is defined similarly.


### 2.1. Variational Form of Micro-Model

Consider (1.1) as a model for single-phase slightly compressible flow in porous media in which the gravity terms have been ignored; see [19] for derivation. Here ϕ=ϕ(x) and k=k(x),x∈Ω denote the nonnegative coefficients of porosity and conductivities, respectively, and f accounts for external sources. We assume for simplicity that k is a scalar quantity.It is well known that (1.1) is well posed if appropriate boundary and initial data are imposed. We will use
(2.1)ρ=ρ0,inΩ×{0},ρ=ρD,on⁡ΓD,k∇ρ⋅n→=0,on⁡ΓN.
For simplicity we assume that ρD≡0 and that the flow is driven by f. However, in simulation examples we use f≡0,  ρD≢0.The variational formulation of problem (1.1) and (2.1) readsFind
ρ
∈
C
(
[
0
,
T
]
;
L
2
(
Ω
)
)
∩
L
2
(
0
,
T
;
H
D
1
(
Ω
)
)  so that(2.2)(ϕ∂ρ∂t,ψ)+(k∇ρ,ψ)=(f,ψ),∀ψ∈HD1(Ω).
For any ρ0∈L2(Ω), f∈C([0,T];L2(Ω)), there is a unique solution ρ∈C([0,T];L2(Ω))∩L2(0,T;HD1(Ω)) to (2.7); see [20, Theorem  11.1.1, page 366]).


### 2.2. Variational Form of Macro-Model

Now consider (1.2a) and (1.2b) and assume that ϕ~1, ϕ~2, k~1, k~2, and c∈L∞(Ω) are bounded uniformly from above and below by positive constants and are independent of time. We consider separately the particular case of (1.3a) and (1.3b) when k~2≡0.Our variational formulation for (1.2a) and (1.2b) uses the theory developed in [17, 18] for evolution problems with m-accretive operators by the Hille-Yosida Theorem. For analysis of classical solutions see [21, 22].For the macro-model (1.2a) and (1.2b) and if k~2>0, we assume the following initial and boundary conditions defined for i=1,2:
(2.3)ρ~i=ρ~i,0,inΩ×{0},ρ~i=ρD,on⁡ΓD,k~i∇ρ~i⋅n→=0,on⁡ΓN.
As before, we assume that ρD≡0 for now.Letρ~i(·,0)=ρ~i,0∈L2(Ω), f∈C([0,T],L2(Ω)), uD≡0. The variational form of (1.2a), (1.2b), and (2.3) is as follows:Find
ρ
~
i
∈
C
(
[
0
,
T
]
;
L
2
(
Ω
)
)
∩
L
2
(
0
,
T
;
H
D
1
(
Ω
)
),  i=1,2, such that(2.4)(ϕ~1∂ρ~1∂t,ψ)+(k~1∇ρ~1,∇ψ)+(c(ρ~1-ρ~2),ψ)=(f,ψ),∀ψ∈HD1(Ω),(ϕ~2∂ρ~2∂t,ξ)+(k~2∇ρ~2,∇ξ)+(c(ρ~2-ρ~1),ξ)=0,∀ξ∈HD1(Ω).
The well posedness of this system (2.4) follows immediately by showing that the operator ℬ:𝒟⊂HD1(Ω)×HD1(Ω)→L2(Ω)×L2(Ω)(2.5)B=[-∂k~1∂+cI-cI-cI-∂k~2∂+cI],whereIis the identity  operator,
is m-accretive. The accretiveness is straightforward by showing that (Bu,u)≥0 in the product space. The maximal property follows from that of the operators -∂k~i∂, I.Now, ifk~2≡0, we do not impose boundary conditions on ρ~2. The variational formulation of (3) and (2.3), readsFind
ρ
~
1
∈
C
(
[
0
,
T
]
;
L
2
(
Ω
)
)
∩
L
2
(
0
,
T
;
H
D
1
(
Ω
)
),  ρ~2∈C([0,T];L2(Ω))∩L2(0,T;L2(Ω))   such that(2.6)(ϕ~1∂ρ~1∂t,ψ)+(k~1∇ρ~1,∇ψ)+(c(ρ~1-ρ~2),ψ)=(f,ψ),∀ψ∈HD1(Ω),(ϕ~2∂ρ~2∂t,ξ)+(c(ρ~2-ρ~1),ξ)=0,∀ξ∈L2(Ω).
The well posedness of (2.6) follows similarly from the Hille-Yosida theorem. Also, ρ~2 has the minimum of regularity of its initial data and of ρ~1.


### 2.3. Notation on Numerical Discretization

For computations of solutions we use conforming (Galerkin) finite elements for spatial discretization and implicit Euler time stepping. We partition the time interval[0,T] into subintervals [tn-1,tn], 1≤n≤N, such that 0=t0<t1<⋯<tN=T. We denote Δtn=tn-tn-1 and z(tn)=zn for any function z.For spatial discretization, we adopt standard finite element nomenclature that can be found in textbooks such as [23–25]. At each time step n, we denote by 𝒯hn, h>0, a family of admissible and shape-regular partitions of Ω into a finite number of elements. For any element T∈𝒯hn, we let hT be the diameter of T and denote h=maxT∈𝒯hhT. Also, ℰhn is the set of all edges of 𝒯hn. For any edge E∈ℰhn, we let hE denote the diameter of the edge E. Finally, we denote by 𝒫k(T) the space of polynomials of degree k in ℝd.Numerical Approximation of Micro-Model
At each time stepn, we define
(2.7)Vhn={vh∈C(Ω¯):∀T∈Thn,vh∣T∈Pk(T),v∣ΓD=0},
and seek an approximation ρhn∈Vhn such that for all ψh∈Vhn(2.8a)(ϕρhn-ρhn-1Δtn,ψh)+(k∇ρhn,∇ψh)=(f,ψh),(2.8b)ρh0=Ihρ0inΩ.Here ℐh denotes an interpolation or projection operator into Vhn.
It is standard [20] that there is a unique solution ρhn∈Vhn for (2.8a) at each n, 1≤n≤N. A priori estimates for the error between the solution to (2.2) and that of (2.8a) and (2.8b) can be found in the literature, see, for example, [20, 25, 26]. For a posteriori error estimates for parabolic problems see [15, 27–29]; the methods in the latter two papers are most relevant to our paper.Numerical Approximation of Macro-Model
Consider now a fully discrete approximation to the solutions to (2.4) or (2.6). The formulation with a fixed grid immediately extends (2.8a) and (2.8b), and appropriate error estimates can be proven. More generally, each of ρ~i can be associated with its own grid 𝒯hin; see details in Section 4.


## 2.1. Variational Form of Micro-Model

Consider (1.1) as a model for single-phase slightly compressible flow in porous media in which the gravity terms have been ignored; see [19] for derivation. Here ϕ=ϕ(x) and k=k(x),x∈Ω denote the nonnegative coefficients of porosity and conductivities, respectively, and f accounts for external sources. We assume for simplicity that k is a scalar quantity.It is well known that (1.1) is well posed if appropriate boundary and initial data are imposed. We will use
(2.1)ρ=ρ0,inΩ×{0},ρ=ρD,on⁡ΓD,k∇ρ⋅n→=0,on⁡ΓN.
For simplicity we assume that ρD≡0 and that the flow is driven by f. However, in simulation examples we use f≡0,  ρD≢0.The variational formulation of problem (1.1) and (2.1) readsFind
ρ
∈
C
(
[
0
,
T
]
;
L
2
(
Ω
)
)
∩
L
2
(
0
,
T
;
H
D
1
(
Ω
)
)  so that(2.2)(ϕ∂ρ∂t,ψ)+(k∇ρ,ψ)=(f,ψ),∀ψ∈HD1(Ω).
For any ρ0∈L2(Ω), f∈C([0,T];L2(Ω)), there is a unique solution ρ∈C([0,T];L2(Ω))∩L2(0,T;HD1(Ω)) to (2.7); see [20, Theorem  11.1.1, page 366]).


## 2.2. Variational Form of Macro-Model

Now consider (1.2a) and (1.2b) and assume that ϕ~1, ϕ~2, k~1, k~2, and c∈L∞(Ω) are bounded uniformly from above and below by positive constants and are independent of time. We consider separately the particular case of (1.3a) and (1.3b) when k~2≡0.Our variational formulation for (1.2a) and (1.2b) uses the theory developed in [17, 18] for evolution problems with m-accretive operators by the Hille-Yosida Theorem. For analysis of classical solutions see [21, 22].For the macro-model (1.2a) and (1.2b) and if k~2>0, we assume the following initial and boundary conditions defined for i=1,2:
(2.3)ρ~i=ρ~i,0,inΩ×{0},ρ~i=ρD,on⁡ΓD,k~i∇ρ~i⋅n→=0,on⁡ΓN.
As before, we assume that ρD≡0 for now.Letρ~i(·,0)=ρ~i,0∈L2(Ω), f∈C([0,T],L2(Ω)), uD≡0. The variational form of (1.2a), (1.2b), and (2.3) is as follows:Find
ρ
~
i
∈
C
(
[
0
,
T
]
;
L
2
(
Ω
)
)
∩
L
2
(
0
,
T
;
H
D
1
(
Ω
)
),  i=1,2, such that(2.4)(ϕ~1∂ρ~1∂t,ψ)+(k~1∇ρ~1,∇ψ)+(c(ρ~1-ρ~2),ψ)=(f,ψ),∀ψ∈HD1(Ω),(ϕ~2∂ρ~2∂t,ξ)+(k~2∇ρ~2,∇ξ)+(c(ρ~2-ρ~1),ξ)=0,∀ξ∈HD1(Ω).
The well posedness of this system (2.4) follows immediately by showing that the operator ℬ:𝒟⊂HD1(Ω)×HD1(Ω)→L2(Ω)×L2(Ω)(2.5)B=[-∂k~1∂+cI-cI-cI-∂k~2∂+cI],whereIis the identity  operator,
is m-accretive. The accretiveness is straightforward by showing that (Bu,u)≥0 in the product space. The maximal property follows from that of the operators -∂k~i∂, I.Now, ifk~2≡0, we do not impose boundary conditions on ρ~2. The variational formulation of (3) and (2.3), readsFind
ρ
~
1
∈
C
(
[
0
,
T
]
;
L
2
(
Ω
)
)
∩
L
2
(
0
,
T
;
H
D
1
(
Ω
)
),  ρ~2∈C([0,T];L2(Ω))∩L2(0,T;L2(Ω))   such that(2.6)(ϕ~1∂ρ~1∂t,ψ)+(k~1∇ρ~1,∇ψ)+(c(ρ~1-ρ~2),ψ)=(f,ψ),∀ψ∈HD1(Ω),(ϕ~2∂ρ~2∂t,ξ)+(c(ρ~2-ρ~1),ξ)=0,∀ξ∈L2(Ω).
The well posedness of (2.6) follows similarly from the Hille-Yosida theorem. Also, ρ~2 has the minimum of regularity of its initial data and of ρ~1.


## 2.3. Notation on Numerical Discretization

For computations of solutions we use conforming (Galerkin) finite elements for spatial discretization and implicit Euler time stepping. We partition the time interval[0,T] into subintervals [tn-1,tn], 1≤n≤N, such that 0=t0<t1<⋯<tN=T. We denote Δtn=tn-tn-1 and z(tn)=zn for any function z.For spatial discretization, we adopt standard finite element nomenclature that can be found in textbooks such as [23–25]. At each time step n, we denote by 𝒯hn, h>0, a family of admissible and shape-regular partitions of Ω into a finite number of elements. For any element T∈𝒯hn, we let hT be the diameter of T and denote h=maxT∈𝒯hhT. Also, ℰhn is the set of all edges of 𝒯hn. For any edge E∈ℰhn, we let hE denote the diameter of the edge E. Finally, we denote by 𝒫k(T) the space of polynomials of degree k in ℝd.Numerical Approximation of Micro-Model
At each time stepn, we define
(2.7)Vhn={vh∈C(Ω¯):∀T∈Thn,vh∣T∈Pk(T),v∣ΓD=0},
and seek an approximation ρhn∈Vhn such that for all ψh∈Vhn(2.8a)(ϕρhn-ρhn-1Δtn,ψh)+(k∇ρhn,∇ψh)=(f,ψh),(2.8b)ρh0=Ihρ0inΩ.Here ℐh denotes an interpolation or projection operator into Vhn.
It is standard [20] that there is a unique solution ρhn∈Vhn for (2.8a) at each n, 1≤n≤N. A priori estimates for the error between the solution to (2.2) and that of (2.8a) and (2.8b) can be found in the literature, see, for example, [20, 25, 26]. For a posteriori error estimates for parabolic problems see [15, 27–29]; the methods in the latter two papers are most relevant to our paper.Numerical Approximation of Macro-Model
Consider now a fully discrete approximation to the solutions to (2.4) or (2.6). The formulation with a fixed grid immediately extends (2.8a) and (2.8b), and appropriate error estimates can be proven. More generally, each of ρ~i can be associated with its own grid 𝒯hin; see details in Section 4.


## 3. Comparison of Micro- and Macro-Models

The difficulties associated with numerical approximation (2.8a) and (2.8b) to (1.1) arise when ϕ, k are highly varying coefficients. In particular, consider a porous domain Ω composed of three disjoint regions Ωi, i=1,…,3, see Figure 1, with positive constants (ϕi,ki)i=13 so that
(3.1)ϕ(x)=ϕ1χ1(x)+ϕ2χ2(x)+ϕ3χ3(x),k(x)=k1χ1(x)+k2χ2(x)+k3χ3(x).
Here χi denotes the characteristic function of Ωi.Figure 1
Periodic heterogeneous media. (a)Ω1: white, Ω2: grey, and Ω3: dark grey. (b) Zoom of one cell Y: Y1: white, Y2: grey, and Y3: dark grey. (c) Y3 in grey with the internal and external boundaries Γ2,3,Γ1,3, respectively. (d) Illustration of three regions as in Example 3.1, where Ω2 is delimited by the white boundaries, and Ω3 is the region between by the black and the white boundaries.An accurate numerical approximation to (1.1) with (3.1) requires that the grid 𝒯h is very fine so that it lines up with the interfaces ∂Ωi∩∂Ωk. The computational complexity of the associated numerical implementation is very large whenever the geometry of interfaces is complex, but this can be overcome by a modeling approximation. In [11] it is shown that (1.2a) and (1.2b) provides solutions close to those of (1.1), with the coefficients ϕ~i,k~i which do not vary locally. Thus a grid for (1.2a) and (1.2b) can be chosen depending on the global macrodynamics rather than local microdynamics and, consequently, it can be much coarser than 𝒯h.


### 3.1. Derivation of Macro-Model

One successful avenue to reduce the complexity of a highly heterogeneous model is to consider ahomogenized or upscaled problem, whose coefficients are derived from the original coefficients [1, 2, 30], and which can be simulated on a coarse grid quite accurately. However, for time-dependent problems, this leads to inaccurate representation of the global dynamics, especially if k1≫k2 and/or k1≫k3.Thedouble-diffusion models proposed in [4], the relaxation model from [6], and double-porosity models derived in [12] are families of successful modeling strategies for parabolic PDEs with periodic highly heterogeneous media. The recent general derivation in [11] is most general and includes those from [4, 6, 12] as special cases.The double-diffusion model developed in [4] was shown in [11] to be a limit of a certain micro-model with three regions arranged periodically. The rigorous derivation in [11] allows for each of the regions Ωi to be connected. Furthermore, it assumes that Ω3 separates Ω1 and Ω2, that is, that ∂Ω1∩∂Ω2=∅. The model in [11] develops a general nonlocal exchange term acting across Ω3. Its special case equivalent to the model in [4] assumes ϕ3≈0 thereby making a quasistatic assumption and allows to compute the coefficient c explicitly.Double-porosity models proposed in [12] assume that only one of the regions Ω1 is connected and so effectively one can set Ω2=∅. However, ϕ3>0. In the macro-model the second equation (1.2b) is replaced by a local micro-model equation, the exchange follows through macro-micro boundary conditions and is nonlocal in nature. We refer to [31–35] for various analyses, models, approximations, and extensions of double-porosity models.Of course, all three regionsΩ1, Ω2, and Ω3 can be connected only in d=3. In d=2 only one of these regions, say Ω1, can be connected; see Figure 1(a). As we point out below, this implies k~2=0, and the macro-model (1.2a) and (1.2b) in d=2 becomes the Warren-Root relaxation model (3) as in [6].In this paper we are interested in the double-diffusion models with a constantc from [4, 6] and its corresponding quasistatic micro-model from [11]. The exchange region Ω3 is assumed to be small and acts like a thick skin separating Ω1 from Ω2. In consequence, the fluid flows effectively only in two regions Ω1 and Ω2; this is represented by the two sub-equations in (1.2a) and (1.2b). Details, comparison, and simulations are shown in what follows.


#### 3.1.1. Macro-Model Parameters

To compute parametersϕ~1, ϕ~2, k~1, k~2, and c of the macro-model (1.2a) and (1.2b), one has to calculate local averages and compute auxiliary solutions of differential equations. These depend on ϕi, ki and on the geometry of Ωε. We follow closely [11].First we formalize the notion of periodicity and heterogeneity inΩ. Without loss of generality we assume that Ω1ε is globally connected and ϕ3=0; all other coefficients in (3.1) are positive constants. We follow the usual structure [1, 11, 12]. Let the unit cube Y=(0,1)d⊂ℝd be divided into three distinct regions Y1, Y2, and Y3, such that ∂Y1∩∂Y2=∅ and Y¯=∪i=13Y¯i as illustrated by Figure 1(b). Denote by Ci the characteristic function of Yi, i=1,2,3, extended Y-periodically to all ℝd. We assume that the domains {x∈ℝd:Ci(x)=1}, i=1,2,3, have a smooth boundary. Also, we define the ε-periodic characteristic functions Ciε(x)=Ci(x/ε), i=1,2,3,  x∈Ω. Thus, for a fixed ε, Ω is subdivided into three distinct regions Ωiε={x∈Ω:Ciε(x)=1}, i=1,2,3. Another subdivision of Ω is into O(ε-d)-periodic cells.Macroporositiesϕ~i
From  [11,equation 16, page 206] we have ϕ~i=∫Yi‍ϕi(y)dy.Macroconductivitiesk~i
From  [11, equations 16, 14, page 205-206] we have
(3.2)k~p(i,j)=∫Yp‍kp(y)(e→i+∇Wip(y))⋅(e→j+∇Wjp(y))dy,p=1,2;i,j=1,…d,
where Wip, for p=1,2;i=1,…,d, is the solution to an auxiliary PDE(3.3a)∇⋅[kp(e→i+∇Wip(y))]=0,inYp,(3.3b)kp(e→i+∇Wip(y))⋅n→p=0,on⁡Γp,3,(3.3c)Wip,kp∇Wip⋅n→pare  periodic.One can prove [1, 3, 11] that k~p is a symmetric matrix.
Now considerd=2. Here the region Ω2ε cannot be connected. Thus, the boundary condition (3.3c) does not apply and (3.3a), (3.3b), and (3.3c) admit the trivial solution e→i+∇Wi2(y)≡0 for i=1,2. Consequently, k~2 is the 2×2 null matrix, and the model (1.2a) and (1.2b) becomes the Warren-Root model (3) proposed in [6].Exchange Term Parameterc
In general, the flow betweenΩ1 and Ω2acrossΩ3 has transient character, and an appropriate term describing it must be nonlocal in nature.
However, for very smallϕ3, the flow is of quasistatic nature as discussed in [4]. The exchange term then has the form c(ρ~1-ρ~2). To compute c, we consider the solution U of(3.4a)-∇⋅[k3∇U]=0,inY3,(3.4b)U=1,on⁡Γ1,3,(3.4c)U=0,on⁡Γ2,3,(3.4d)U,k3∇U⋅n→3are  periodic onΓ3,3.The parameter c see  [11, equation 19, page 209] is given by
(3.5)c=∫Γ1,3‍k3∇yU⋅n→3ds.
We note that the general non-quasistatic case is when ϕ3(∂U/∂t) is included in (3.4d) thereby changing the character of exchange term to nonlocal. Conversely, (3.4a), (3.4b), (3.4c), and (3.4d) can be considered as its special case when ϕ3≈0.


#### 3.1.2. Numerical Computation of Macro-Model Parameters

The solutions of the auxiliary PDEs (3.3a), (3.3b), and (3.3c) and (3.4a), (3.4b), (3.4c), and (3.4d)  are approximated using finite elements with hk1-1=40 and hc-1=20. See [36] for treatment of periodic boundary conditions.In what follows we ignore the distinction between the exact and numerical values of macro-model parametersk~i, c and provide examples of their calculations for various values of k3 and choices of geometry of Y.In all examplesd=2 hence k~2=[0000]. Also, in all examples k1≡const and the geometry of Y is axisymmetric, and thus k~1 is a scalar constant.Example 3.1.
HereY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2,  and k=1χ1+10-4χ2+10-7χ3. We obtain
(3.6)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-6.Example 3.2.
LetY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2. ϕ=1χ1+10-4χ2, and  k=1χ1+10-1χ2+10-4χ3. We obtain
(3.7)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-3.Example 3.3.
HereY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2,  and k=1χ1+10-1χ2+10-2χ3. We obtain
(3.8)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-1.Example 3.4.
LetY2=(0.3,0.7)2,  Y3=(0.1,0.9)2∖Y2,  ϕ=1χ1+10-4χ2, and k=1χ1+10-1χ2+10-4χ3. We get
(3.9)ϕ~1=0.36,ϕ~2=1.6×10-5,k~1=0.2130,c=1.0412×10-3.Example 3.5.
Y
2
=
(
0.25,0.75
)
2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2, and k=1χ1+10-1χ2+10-4χ3. We get
(3.10)ϕ~1=0.64,ϕ~2=2.5×10-5,k~1=0.4519,c=4.3000×10-3.


### 3.2. Comparison of the Models

Now we can compare the solutions of the micro-model (1.1) to those of the macro-model (1.2a) and (1.2b). In [11] it was shown that the latter two-scale converge to the former as ε→0. In addition, the analysis in [11] suggests that it is appropriate to compare ρ~i to χiρε. However, this notion of convergence does not give information about the rate at which ρ~i-χiρε may converge, and it involves special periodic test functions.In this paper we estimate this rate by comparing their numerical approximations directly without the use of any test functions but rather in a certain metric of interest. To the best of our knowledge such comparison or convergence rate was not discussed elsewhere.Strictly speaking, the two-scale convergence proof considered in [11], and a similar proof for the double-porosity model in [12], includes scaling of k3 with ε2. This scaling is a formal device necessary to preserve certain parts of the boundary value problem under investigation in the limit as ε→0. However, the homogenization limit is intended to serve only as an approximation of the true model, which has a given fixed set of parameters. We do not include this scaling in our computations. Rather, we treat each example, for a given ε, as a data set in its own merit, rather than as an element of a sequence intended to two-scale converge to the limit.


#### 3.2.1. Setup of Computational Experiments

We set up simulations for comparison using compatible data for micro- and macro-models; we setf≡0 and choose initial and boundary data driving the flow with interesting dynamics. For porosities and conductivities, we use the values computed in Section 3.1.2.LetΩ=(0,1)2 and ΓD={0}×[0,1]∪{1}×[0,1]. On ΓD we define ρD(0,y)=1,ρD(1,y)=0,y∈[0,1]. On the lateral sides of Ω the Neumann no-flow condition is imposed. Also, let ρ0≡1. Thus the flow in the micro-model goes from left to right, and the solution evolves towards the stationary solution 〈1-x,1-x〉. Due to the incompatibility between initial and boundary data, we have high gradients of the solution close to x=1 for small t.We fix the final timeT=0.05 and use uniform time stepping with Δt=10-4. The solution of the micro-model ρ depends on the number of cells in Ω≡Ωε, with ε∈(0,1]. For example, if ε=1/2, the domain Ω is composed of 2×2=4 cells. We denote the numerical solution of (2.8a) and (2.8b) at tn by ρε,hn where h denotes the grid parameter for the micro-model.For the macro-model, we useρ~i,0=ρ0≡1. Also, we use ρ~D,i=ρD for i=1,2 when k~2≢0. If k~2≡0, then the boundary condition for ρ~2 is not prescribed.


#### 3.2.2. Qualitative Comparison

In our first comparison we use parameters from Example3.1. We solve the macro-model (4.2) for {〈ρ~1,h1n,ρ~2,h2n〉}n=1N and the micro-model (2.8a) and (2.8b) for ρε,hn. The plots of ρhn and ρ~i,hn are in Figure 2, and a different view is shown in Figure 1(d).Numerical approximation attn=3×10-2 of (a)-(b) the solutions of the two components of the macro-model (4.2), and (c) the solutions of the micro-model (2.8a) and (2.8b). Figure (d) shows the side view of (a, b, c) at y=0.4 and y=0.5. Data from Example 3.1 with 1/ε=5 and h=h1=h2=1/50. Note that (a) and (b) have different vertical scales.
(a)
(b)
(c)
(d)The heterogeneous structure ofΩε is well visible from the behavior of the micro-model solution ρhn. Large gradients of solution are visible on cell boundaries due to the large difference between k1, k2, and k3.However, the behavior in the fast region is well approximated byρ~1,h1n which envelopes ρhnχ1, while ρ~2,h2n envelopes ρhnχ2 well. This is very well seen in a side view in Figure 2(d).Quite interesting is the behavior ofρ2|x=1. Since in our examples k~2≡0, no boundary condition at x=1 is imposed on ρ2. Thus, for small t, ρ2(1,t) evolves from its initial constant value of 1 according to (1.3b) with the input from ρ~1. The latter satisfies, however, the homogeneous boundary condition at x=1. In particular, with constant ϕ~2, c, we have
(3.11)ρ~2(x,y,t)=ρ~2,0(x,y)e-(c/ϕ~2)t+∫0t‍ρ~1(x,y,s)e-(c/ϕ~2)(t-s)ds.
Thus, for small t and small c/ϕ~2, ρ~2(1,y,t) is away from 0, but, as time increases, its magnitude decreases proportionally to e-(c/ϕ~2)t.Next, we use parameters from Example3.2 and plot them in Figure 3. In this example the conductivity k1 is much larger than that of previous case, the local gradients in the micro-model are smaller, and the solutions to the macro-model achieve a smoother profile faster.Numerical approximation attn=3×10-2 of (a) and (b) the solutions of the two components of the macro-model (4.2) and (c) the solutions of the micro-model (2.8a) and (2.8b). Figure (d) shows the side view of (a, b, c) at y=0.4 and y=0.5. Data from Example 3.2 with 1/ε=5 and h=h1=h2=1/50.
(a)
(b)
(c)
(d)


#### 3.2.3. Quantitative Comparison—Fine Grid in Macro-Model

Now we are ready to discuss a quantitative comparison between solutions to the micro- and macro-models. For this, we set up a family of cases withε=1/(2l+1),  l≥1. We use uniform mesh with h=h1=h2=ε/20 for now. See Section 3.2.4 for hi≫h and Section 4 for adaptive gridding.From  [11, Theorem  18, page 208]  and from our plots it follows that we can actually compare the micro- and macrosolutions, as long as the characteristic functions χi are involved. While the results in [11] use two-scale convergence with ε, that is, rely on special periodic test functions, we compare ρε,hn and ρ~i,hn directly. We see that it is easy to do so only in the connected region Ω1.We define a quantityeεn=(ρε,hn-ρ~1,h1n)χ1 to be used in comparison. We also define, for a fixed t∈(0,T] and some function z,
(3.12)‖z(x,y;t)‖*,L1:=∫01‍|z(x,0.5;t)|dx,‖z(x,y;t)‖*,L2:={∫01‍|z(x,0.5;t)|2dx}1/2,
which will be applied to eεn. Clearly ∥·∥*,W is not a norm but is a useful quantity of interest.In Tables1, 2, 3, and 4 we present ∥eεn∥*,W for different ε and W=L2,L1. For each case we consider convergence of ∥eεn∥*,W with ε↓0.Table 1
Comparison of solutions to the macro-model with those to the micro-model for different values ofε with data as in Example 3.2.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1832
×
1
0
-
2
—
3.2071
×
1
0
-
2
—
5
×
5
1
/
5
7.0915
×
1
0
-
3
1.0188
1.9044
×
1
0
-
2
1.0390
7
×
7
1
/
7
5.0522
×
1
0
-
3
1.0077
1.3528
×
1
0
-
2
1.0164
9
×
9
1
/
9
3.9285
×
1
0
-
3
1.0011
1.0502
×
1
0
-
2
1.0075
11
×
11
1
/
11
3.2162
×
1
0
-
3
0.9969
8.5890
×
1
0
-
3
1.0021Table 2
Comparison of the macro-model to the micro-model for different values ofε as in Example 3.3.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1446
×
1
0
-
2
—
3.0715
×
1
0
-
2
—
5
×
5
1
/
5
7.0065
×
1
0
-
3
1.0408
1.8387
×
1
0
-
2
0.9956
7
×
7
1
/
7
5.1541
×
1
0
-
3
1.0958
1.3363
×
1
0
-
2
1.0543
9
×
9
1
/
9
4.1704
×
1
0
-
3
1.1867
1.0684
×
1
0
-
2
1.1232
11
×
11
1
/
11
3.5794
×
1
0
-
3
1.3131
9.0460
×
1
0
-
3
1.2058Table 3
Comparison of the macro-model to the micro-model for different values ofε Example 3.4.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.3758
×
1
0
-
2
—
3.5390
×
1
0
-
2
—
5
×
5
1
/
5
8.2810
×
1
0
-
3
1.0330
2.0998
×
1
0
-
2
0.9786
7
×
7
1
/
7
5.7341
×
1
0
-
3
0.9155
1.4741
×
1
0
-
2
0.9510
9
×
9
1
/
9
4.3987
×
1
0
-
3
0.9479
1.1336
×
1
0
-
2
0.9568Table 4
Comparison of the macro-model to the micro-model for different values ofε Example 3.5.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1826
×
1
0
-
2
—
3.2048
×
1
0
-
2
—
5
×
5
1
/
5
7.0877
×
1
0
-
3
0.9978
1.9029
×
1
0
-
2
0.9800
7
×
7
1
/
7
5.0497
×
1
0
-
3
0.9924
1.3518
×
1
0
-
2
0.9840
9
×
9
1
/
9
3.9281
×
1
0
-
3
1.0006
1.0502
×
1
0
-
2
0.9955From homogenization theory it is not clear what order of convergence should be expected. We observe that∥eεn∥*,W decreases linearly with ε for both *=L1,L2 -like quantities in all cases.Next we discuss the dependence of the results on conductivity and other data. Consider Example3.2 as a reference example with its corresponding Table 1. Compare it now with that from Table 2 where k3 is larger than that in the reference case to see how this influences the errors. As expected, and as shown in Example 3.3, c changes by the same factor, and this means that the coupling between ρ~1 and ρ~2 is stronger which corresponds to faster flow across Y3 in the micro-model. However, the influence on the convergence of the micro-macro difference, from Table 2 is rather weak.Next consider Example3.4 which uses a different geometry than that in the reference case, with thicker region Y3. This produces a slightly smaller c since the gradient of U is smaller. Also, ϕ~1 and k~1 are predictably different. The error quantity is influenced insignificantly.The opposite effect is seen in Table4 when Y3 is selected as thinner than in the reference Example 3.2.


#### 3.2.4. Coarse Macrogrid and Parameter Grid

In the examples presented in Section3.2.3 we used h=h1=h2. This choice of compatible grids is convenient for visualization purposes. However, for the idea of the macro-model to be useful, the computational complexity of the macro-model needs to be lower than that of the micro-model. Since the macro-model is a system of two equations, its lower complexity can be only achieved if the grid in the macro-model is significantly coarser than that of micro-model. Below we show that one can use h≪h1 in the macro-model.Example 3.6.
We use data from Example3.2 with ε=0.2. Consider a fixed microgrid with h-1=100 and a few cases of macrogrid with h1-1=h2-1=10,50,100. Now we compare the micro- and macrosolutions ρε,hn and ρ1,h1, see Table 5. We also compare the run time of the macro-model to that of the reference micro-model which is 120 seconds.Table 5
Comparison of the micro-model to the macro-model as in Example3.6.
h
1
-
1
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
run time(s)
100
7.0915
×
1
0
-
3
1.9044
×
1
0
-
2
511
50
7.0937
×
1
0
-
3
1.9048
×
1
0
-
2
56
10
7.2380
×
1
0
-
3
1.9173
×
1
0
-
2
3We see that, as predicted above, the macro-model run time athi=h is not competitive with that of the micro-model; this is exacerbated by the overhead of various adaptive procedures to be discussed below. However, when hi≫h, the ability of the macro-model solutions to approximate those of micro-model appears reasonable and is associated with a much lower cost. This is true even when hi is by a factor of 10 coarser than h resulting in a computational time decreased by two orders of magnitude.Clearly, the macro-model can run much faster than the micro-model especially when adaptive nonuniform grids withh1≠h2 are implemented, see Section 4.Influence of Parameter Grid
Last we address the effect of the computed parametersc, k~1 on the closeness of the macro-model solution to that of micro-model. These coefficients are precomputed as discussed in Section 3.1.1.Example 3.7.
Letε=1/3, and let the micro-model grid be h-1=60. Use h1-1=h2-1=60. Compute numerical parameters chc, k~1,hk1 by approximating solutions to (3.4a), (3.4b), (3.4c), and (3.4d)  and (3.3a), (3.3b), and (3.3c), respectively.
Example3.7 and Table 6 show that the grid used for computing these parameters only mildly affects the quality of the solution, thus one can use coarse grid with confidence.Table 6
Sensitivity of the macro-model solution to the grid used to computechc≈c and k~1,hk1≈k~1 as Example 3.7.
h
c
-
1
c
h
c
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
10
1.9000
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
50
1.8370
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
100
1.8292
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
200
1.8258
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
h
k
-
1
k
~
1
(
1,1
)
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
10
4.6306
×
1
0
-
1
3.5392
×
1
0
-
2
5.7835
×
1
0
-
2
50
4.5141
×
1
0
-
1
3.6222
×
1
0
-
2
5.8989
×
1
0
-
2
100
4.5046
×
1
0
-
1
3.6301
×
1
0
-
2
5.9093
×
1
0
-
2
200
4.5009
×
1
0
-
1
3.6332
×
1
0
-
2
5.9134
×
1
0
-
2


## 3.1. Derivation of Macro-Model

One successful avenue to reduce the complexity of a highly heterogeneous model is to consider ahomogenized or upscaled problem, whose coefficients are derived from the original coefficients [1, 2, 30], and which can be simulated on a coarse grid quite accurately. However, for time-dependent problems, this leads to inaccurate representation of the global dynamics, especially if k1≫k2 and/or k1≫k3.Thedouble-diffusion models proposed in [4], the relaxation model from [6], and double-porosity models derived in [12] are families of successful modeling strategies for parabolic PDEs with periodic highly heterogeneous media. The recent general derivation in [11] is most general and includes those from [4, 6, 12] as special cases.The double-diffusion model developed in [4] was shown in [11] to be a limit of a certain micro-model with three regions arranged periodically. The rigorous derivation in [11] allows for each of the regions Ωi to be connected. Furthermore, it assumes that Ω3 separates Ω1 and Ω2, that is, that ∂Ω1∩∂Ω2=∅. The model in [11] develops a general nonlocal exchange term acting across Ω3. Its special case equivalent to the model in [4] assumes ϕ3≈0 thereby making a quasistatic assumption and allows to compute the coefficient c explicitly.Double-porosity models proposed in [12] assume that only one of the regions Ω1 is connected and so effectively one can set Ω2=∅. However, ϕ3>0. In the macro-model the second equation (1.2b) is replaced by a local micro-model equation, the exchange follows through macro-micro boundary conditions and is nonlocal in nature. We refer to [31–35] for various analyses, models, approximations, and extensions of double-porosity models.Of course, all three regionsΩ1, Ω2, and Ω3 can be connected only in d=3. In d=2 only one of these regions, say Ω1, can be connected; see Figure 1(a). As we point out below, this implies k~2=0, and the macro-model (1.2a) and (1.2b) in d=2 becomes the Warren-Root relaxation model (3) as in [6].In this paper we are interested in the double-diffusion models with a constantc from [4, 6] and its corresponding quasistatic micro-model from [11]. The exchange region Ω3 is assumed to be small and acts like a thick skin separating Ω1 from Ω2. In consequence, the fluid flows effectively only in two regions Ω1 and Ω2; this is represented by the two sub-equations in (1.2a) and (1.2b). Details, comparison, and simulations are shown in what follows.


### 3.1.1. Macro-Model Parameters

To compute parametersϕ~1, ϕ~2, k~1, k~2, and c of the macro-model (1.2a) and (1.2b), one has to calculate local averages and compute auxiliary solutions of differential equations. These depend on ϕi, ki and on the geometry of Ωε. We follow closely [11].First we formalize the notion of periodicity and heterogeneity inΩ. Without loss of generality we assume that Ω1ε is globally connected and ϕ3=0; all other coefficients in (3.1) are positive constants. We follow the usual structure [1, 11, 12]. Let the unit cube Y=(0,1)d⊂ℝd be divided into three distinct regions Y1, Y2, and Y3, such that ∂Y1∩∂Y2=∅ and Y¯=∪i=13Y¯i as illustrated by Figure 1(b). Denote by Ci the characteristic function of Yi, i=1,2,3, extended Y-periodically to all ℝd. We assume that the domains {x∈ℝd:Ci(x)=1}, i=1,2,3, have a smooth boundary. Also, we define the ε-periodic characteristic functions Ciε(x)=Ci(x/ε), i=1,2,3,  x∈Ω. Thus, for a fixed ε, Ω is subdivided into three distinct regions Ωiε={x∈Ω:Ciε(x)=1}, i=1,2,3. Another subdivision of Ω is into O(ε-d)-periodic cells.Macroporositiesϕ~i
From  [11,equation 16, page 206] we have ϕ~i=∫Yi‍ϕi(y)dy.Macroconductivitiesk~i
From  [11, equations 16, 14, page 205-206] we have
(3.2)k~p(i,j)=∫Yp‍kp(y)(e→i+∇Wip(y))⋅(e→j+∇Wjp(y))dy,p=1,2;i,j=1,…d,
where Wip, for p=1,2;i=1,…,d, is the solution to an auxiliary PDE(3.3a)∇⋅[kp(e→i+∇Wip(y))]=0,inYp,(3.3b)kp(e→i+∇Wip(y))⋅n→p=0,on⁡Γp,3,(3.3c)Wip,kp∇Wip⋅n→pare  periodic.One can prove [1, 3, 11] that k~p is a symmetric matrix.
Now considerd=2. Here the region Ω2ε cannot be connected. Thus, the boundary condition (3.3c) does not apply and (3.3a), (3.3b), and (3.3c) admit the trivial solution e→i+∇Wi2(y)≡0 for i=1,2. Consequently, k~2 is the 2×2 null matrix, and the model (1.2a) and (1.2b) becomes the Warren-Root model (3) proposed in [6].Exchange Term Parameterc
In general, the flow betweenΩ1 and Ω2acrossΩ3 has transient character, and an appropriate term describing it must be nonlocal in nature.
However, for very smallϕ3, the flow is of quasistatic nature as discussed in [4]. The exchange term then has the form c(ρ~1-ρ~2). To compute c, we consider the solution U of(3.4a)-∇⋅[k3∇U]=0,inY3,(3.4b)U=1,on⁡Γ1,3,(3.4c)U=0,on⁡Γ2,3,(3.4d)U,k3∇U⋅n→3are  periodic onΓ3,3.The parameter c see  [11, equation 19, page 209] is given by
(3.5)c=∫Γ1,3‍k3∇yU⋅n→3ds.
We note that the general non-quasistatic case is when ϕ3(∂U/∂t) is included in (3.4d) thereby changing the character of exchange term to nonlocal. Conversely, (3.4a), (3.4b), (3.4c), and (3.4d) can be considered as its special case when ϕ3≈0.


### 3.1.2. Numerical Computation of Macro-Model Parameters

The solutions of the auxiliary PDEs (3.3a), (3.3b), and (3.3c) and (3.4a), (3.4b), (3.4c), and (3.4d)  are approximated using finite elements with hk1-1=40 and hc-1=20. See [36] for treatment of periodic boundary conditions.In what follows we ignore the distinction between the exact and numerical values of macro-model parametersk~i, c and provide examples of their calculations for various values of k3 and choices of geometry of Y.In all examplesd=2 hence k~2=[0000]. Also, in all examples k1≡const and the geometry of Y is axisymmetric, and thus k~1 is a scalar constant.Example 3.1.
HereY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2,  and k=1χ1+10-4χ2+10-7χ3. We obtain
(3.6)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-6.Example 3.2.
LetY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2. ϕ=1χ1+10-4χ2, and  k=1χ1+10-1χ2+10-4χ3. We obtain
(3.7)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-3.Example 3.3.
HereY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2,  and k=1χ1+10-1χ2+10-2χ3. We obtain
(3.8)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-1.Example 3.4.
LetY2=(0.3,0.7)2,  Y3=(0.1,0.9)2∖Y2,  ϕ=1χ1+10-4χ2, and k=1χ1+10-1χ2+10-4χ3. We get
(3.9)ϕ~1=0.36,ϕ~2=1.6×10-5,k~1=0.2130,c=1.0412×10-3.Example 3.5.
Y
2
=
(
0.25,0.75
)
2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2, and k=1χ1+10-1χ2+10-4χ3. We get
(3.10)ϕ~1=0.64,ϕ~2=2.5×10-5,k~1=0.4519,c=4.3000×10-3.


## 3.1.1. Macro-Model Parameters

To compute parametersϕ~1, ϕ~2, k~1, k~2, and c of the macro-model (1.2a) and (1.2b), one has to calculate local averages and compute auxiliary solutions of differential equations. These depend on ϕi, ki and on the geometry of Ωε. We follow closely [11].First we formalize the notion of periodicity and heterogeneity inΩ. Without loss of generality we assume that Ω1ε is globally connected and ϕ3=0; all other coefficients in (3.1) are positive constants. We follow the usual structure [1, 11, 12]. Let the unit cube Y=(0,1)d⊂ℝd be divided into three distinct regions Y1, Y2, and Y3, such that ∂Y1∩∂Y2=∅ and Y¯=∪i=13Y¯i as illustrated by Figure 1(b). Denote by Ci the characteristic function of Yi, i=1,2,3, extended Y-periodically to all ℝd. We assume that the domains {x∈ℝd:Ci(x)=1}, i=1,2,3, have a smooth boundary. Also, we define the ε-periodic characteristic functions Ciε(x)=Ci(x/ε), i=1,2,3,  x∈Ω. Thus, for a fixed ε, Ω is subdivided into three distinct regions Ωiε={x∈Ω:Ciε(x)=1}, i=1,2,3. Another subdivision of Ω is into O(ε-d)-periodic cells.Macroporositiesϕ~i
From  [11,equation 16, page 206] we have ϕ~i=∫Yi‍ϕi(y)dy.Macroconductivitiesk~i
From  [11, equations 16, 14, page 205-206] we have
(3.2)k~p(i,j)=∫Yp‍kp(y)(e→i+∇Wip(y))⋅(e→j+∇Wjp(y))dy,p=1,2;i,j=1,…d,
where Wip, for p=1,2;i=1,…,d, is the solution to an auxiliary PDE(3.3a)∇⋅[kp(e→i+∇Wip(y))]=0,inYp,(3.3b)kp(e→i+∇Wip(y))⋅n→p=0,on⁡Γp,3,(3.3c)Wip,kp∇Wip⋅n→pare  periodic.One can prove [1, 3, 11] that k~p is a symmetric matrix.
Now considerd=2. Here the region Ω2ε cannot be connected. Thus, the boundary condition (3.3c) does not apply and (3.3a), (3.3b), and (3.3c) admit the trivial solution e→i+∇Wi2(y)≡0 for i=1,2. Consequently, k~2 is the 2×2 null matrix, and the model (1.2a) and (1.2b) becomes the Warren-Root model (3) proposed in [6].Exchange Term Parameterc
In general, the flow betweenΩ1 and Ω2acrossΩ3 has transient character, and an appropriate term describing it must be nonlocal in nature.
However, for very smallϕ3, the flow is of quasistatic nature as discussed in [4]. The exchange term then has the form c(ρ~1-ρ~2). To compute c, we consider the solution U of(3.4a)-∇⋅[k3∇U]=0,inY3,(3.4b)U=1,on⁡Γ1,3,(3.4c)U=0,on⁡Γ2,3,(3.4d)U,k3∇U⋅n→3are  periodic onΓ3,3.The parameter c see  [11, equation 19, page 209] is given by
(3.5)c=∫Γ1,3‍k3∇yU⋅n→3ds.
We note that the general non-quasistatic case is when ϕ3(∂U/∂t) is included in (3.4d) thereby changing the character of exchange term to nonlocal. Conversely, (3.4a), (3.4b), (3.4c), and (3.4d) can be considered as its special case when ϕ3≈0.


## 3.1.2. Numerical Computation of Macro-Model Parameters

The solutions of the auxiliary PDEs (3.3a), (3.3b), and (3.3c) and (3.4a), (3.4b), (3.4c), and (3.4d)  are approximated using finite elements with hk1-1=40 and hc-1=20. See [36] for treatment of periodic boundary conditions.In what follows we ignore the distinction between the exact and numerical values of macro-model parametersk~i, c and provide examples of their calculations for various values of k3 and choices of geometry of Y.In all examplesd=2 hence k~2=[0000]. Also, in all examples k1≡const and the geometry of Y is axisymmetric, and thus k~1 is a scalar constant.Example 3.1.
HereY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2,  and k=1χ1+10-4χ2+10-7χ3. We obtain
(3.6)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-6.Example 3.2.
LetY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2. ϕ=1χ1+10-4χ2, and  k=1χ1+10-1χ2+10-4χ3. We obtain
(3.7)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-3.Example 3.3.
HereY2=(0.3,0.7)2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2,  and k=1χ1+10-1χ2+10-2χ3. We obtain
(3.8)ϕ~1=0.64,ϕ~2=1.6×10-5,k~1=0.4519,c=1.8634×10-1.Example 3.4.
LetY2=(0.3,0.7)2,  Y3=(0.1,0.9)2∖Y2,  ϕ=1χ1+10-4χ2, and k=1χ1+10-1χ2+10-4χ3. We get
(3.9)ϕ~1=0.36,ϕ~2=1.6×10-5,k~1=0.2130,c=1.0412×10-3.Example 3.5.
Y
2
=
(
0.25,0.75
)
2,  Y3=(0.2,0.8)2∖Y2,  ϕ=1χ1+10-4χ2, and k=1χ1+10-1χ2+10-4χ3. We get
(3.10)ϕ~1=0.64,ϕ~2=2.5×10-5,k~1=0.4519,c=4.3000×10-3.


## 3.2. Comparison of the Models

Now we can compare the solutions of the micro-model (1.1) to those of the macro-model (1.2a) and (1.2b). In [11] it was shown that the latter two-scale converge to the former as ε→0. In addition, the analysis in [11] suggests that it is appropriate to compare ρ~i to χiρε. However, this notion of convergence does not give information about the rate at which ρ~i-χiρε may converge, and it involves special periodic test functions.In this paper we estimate this rate by comparing their numerical approximations directly without the use of any test functions but rather in a certain metric of interest. To the best of our knowledge such comparison or convergence rate was not discussed elsewhere.Strictly speaking, the two-scale convergence proof considered in [11], and a similar proof for the double-porosity model in [12], includes scaling of k3 with ε2. This scaling is a formal device necessary to preserve certain parts of the boundary value problem under investigation in the limit as ε→0. However, the homogenization limit is intended to serve only as an approximation of the true model, which has a given fixed set of parameters. We do not include this scaling in our computations. Rather, we treat each example, for a given ε, as a data set in its own merit, rather than as an element of a sequence intended to two-scale converge to the limit.


### 3.2.1. Setup of Computational Experiments

We set up simulations for comparison using compatible data for micro- and macro-models; we setf≡0 and choose initial and boundary data driving the flow with interesting dynamics. For porosities and conductivities, we use the values computed in Section 3.1.2.LetΩ=(0,1)2 and ΓD={0}×[0,1]∪{1}×[0,1]. On ΓD we define ρD(0,y)=1,ρD(1,y)=0,y∈[0,1]. On the lateral sides of Ω the Neumann no-flow condition is imposed. Also, let ρ0≡1. Thus the flow in the micro-model goes from left to right, and the solution evolves towards the stationary solution 〈1-x,1-x〉. Due to the incompatibility between initial and boundary data, we have high gradients of the solution close to x=1 for small t.We fix the final timeT=0.05 and use uniform time stepping with Δt=10-4. The solution of the micro-model ρ depends on the number of cells in Ω≡Ωε, with ε∈(0,1]. For example, if ε=1/2, the domain Ω is composed of 2×2=4 cells. We denote the numerical solution of (2.8a) and (2.8b) at tn by ρε,hn where h denotes the grid parameter for the micro-model.For the macro-model, we useρ~i,0=ρ0≡1. Also, we use ρ~D,i=ρD for i=1,2 when k~2≢0. If k~2≡0, then the boundary condition for ρ~2 is not prescribed.


### 3.2.2. Qualitative Comparison

In our first comparison we use parameters from Example3.1. We solve the macro-model (4.2) for {〈ρ~1,h1n,ρ~2,h2n〉}n=1N and the micro-model (2.8a) and (2.8b) for ρε,hn. The plots of ρhn and ρ~i,hn are in Figure 2, and a different view is shown in Figure 1(d).Numerical approximation attn=3×10-2 of (a)-(b) the solutions of the two components of the macro-model (4.2), and (c) the solutions of the micro-model (2.8a) and (2.8b). Figure (d) shows the side view of (a, b, c) at y=0.4 and y=0.5. Data from Example 3.1 with 1/ε=5 and h=h1=h2=1/50. Note that (a) and (b) have different vertical scales.
(a)
(b)
(c)
(d)The heterogeneous structure ofΩε is well visible from the behavior of the micro-model solution ρhn. Large gradients of solution are visible on cell boundaries due to the large difference between k1, k2, and k3.However, the behavior in the fast region is well approximated byρ~1,h1n which envelopes ρhnχ1, while ρ~2,h2n envelopes ρhnχ2 well. This is very well seen in a side view in Figure 2(d).Quite interesting is the behavior ofρ2|x=1. Since in our examples k~2≡0, no boundary condition at x=1 is imposed on ρ2. Thus, for small t, ρ2(1,t) evolves from its initial constant value of 1 according to (1.3b) with the input from ρ~1. The latter satisfies, however, the homogeneous boundary condition at x=1. In particular, with constant ϕ~2, c, we have
(3.11)ρ~2(x,y,t)=ρ~2,0(x,y)e-(c/ϕ~2)t+∫0t‍ρ~1(x,y,s)e-(c/ϕ~2)(t-s)ds.
Thus, for small t and small c/ϕ~2, ρ~2(1,y,t) is away from 0, but, as time increases, its magnitude decreases proportionally to e-(c/ϕ~2)t.Next, we use parameters from Example3.2 and plot them in Figure 3. In this example the conductivity k1 is much larger than that of previous case, the local gradients in the micro-model are smaller, and the solutions to the macro-model achieve a smoother profile faster.Numerical approximation attn=3×10-2 of (a) and (b) the solutions of the two components of the macro-model (4.2) and (c) the solutions of the micro-model (2.8a) and (2.8b). Figure (d) shows the side view of (a, b, c) at y=0.4 and y=0.5. Data from Example 3.2 with 1/ε=5 and h=h1=h2=1/50.
(a)
(b)
(c)
(d)


### 3.2.3. Quantitative Comparison—Fine Grid in Macro-Model

Now we are ready to discuss a quantitative comparison between solutions to the micro- and macro-models. For this, we set up a family of cases withε=1/(2l+1),  l≥1. We use uniform mesh with h=h1=h2=ε/20 for now. See Section 3.2.4 for hi≫h and Section 4 for adaptive gridding.From  [11, Theorem  18, page 208]  and from our plots it follows that we can actually compare the micro- and macrosolutions, as long as the characteristic functions χi are involved. While the results in [11] use two-scale convergence with ε, that is, rely on special periodic test functions, we compare ρε,hn and ρ~i,hn directly. We see that it is easy to do so only in the connected region Ω1.We define a quantityeεn=(ρε,hn-ρ~1,h1n)χ1 to be used in comparison. We also define, for a fixed t∈(0,T] and some function z,
(3.12)‖z(x,y;t)‖*,L1:=∫01‍|z(x,0.5;t)|dx,‖z(x,y;t)‖*,L2:={∫01‍|z(x,0.5;t)|2dx}1/2,
which will be applied to eεn. Clearly ∥·∥*,W is not a norm but is a useful quantity of interest.In Tables1, 2, 3, and 4 we present ∥eεn∥*,W for different ε and W=L2,L1. For each case we consider convergence of ∥eεn∥*,W with ε↓0.Table 1
Comparison of solutions to the macro-model with those to the micro-model for different values ofε with data as in Example 3.2.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1832
×
1
0
-
2
—
3.2071
×
1
0
-
2
—
5
×
5
1
/
5
7.0915
×
1
0
-
3
1.0188
1.9044
×
1
0
-
2
1.0390
7
×
7
1
/
7
5.0522
×
1
0
-
3
1.0077
1.3528
×
1
0
-
2
1.0164
9
×
9
1
/
9
3.9285
×
1
0
-
3
1.0011
1.0502
×
1
0
-
2
1.0075
11
×
11
1
/
11
3.2162
×
1
0
-
3
0.9969
8.5890
×
1
0
-
3
1.0021Table 2
Comparison of the macro-model to the micro-model for different values ofε as in Example 3.3.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1446
×
1
0
-
2
—
3.0715
×
1
0
-
2
—
5
×
5
1
/
5
7.0065
×
1
0
-
3
1.0408
1.8387
×
1
0
-
2
0.9956
7
×
7
1
/
7
5.1541
×
1
0
-
3
1.0958
1.3363
×
1
0
-
2
1.0543
9
×
9
1
/
9
4.1704
×
1
0
-
3
1.1867
1.0684
×
1
0
-
2
1.1232
11
×
11
1
/
11
3.5794
×
1
0
-
3
1.3131
9.0460
×
1
0
-
3
1.2058Table 3
Comparison of the macro-model to the micro-model for different values ofε Example 3.4.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.3758
×
1
0
-
2
—
3.5390
×
1
0
-
2
—
5
×
5
1
/
5
8.2810
×
1
0
-
3
1.0330
2.0998
×
1
0
-
2
0.9786
7
×
7
1
/
7
5.7341
×
1
0
-
3
0.9155
1.4741
×
1
0
-
2
0.9510
9
×
9
1
/
9
4.3987
×
1
0
-
3
0.9479
1.1336
×
1
0
-
2
0.9568Table 4
Comparison of the macro-model to the micro-model for different values ofε Example 3.5.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1826
×
1
0
-
2
—
3.2048
×
1
0
-
2
—
5
×
5
1
/
5
7.0877
×
1
0
-
3
0.9978
1.9029
×
1
0
-
2
0.9800
7
×
7
1
/
7
5.0497
×
1
0
-
3
0.9924
1.3518
×
1
0
-
2
0.9840
9
×
9
1
/
9
3.9281
×
1
0
-
3
1.0006
1.0502
×
1
0
-
2
0.9955From homogenization theory it is not clear what order of convergence should be expected. We observe that∥eεn∥*,W decreases linearly with ε for both *=L1,L2 -like quantities in all cases.Next we discuss the dependence of the results on conductivity and other data. Consider Example3.2 as a reference example with its corresponding Table 1. Compare it now with that from Table 2 where k3 is larger than that in the reference case to see how this influences the errors. As expected, and as shown in Example 3.3, c changes by the same factor, and this means that the coupling between ρ~1 and ρ~2 is stronger which corresponds to faster flow across Y3 in the micro-model. However, the influence on the convergence of the micro-macro difference, from Table 2 is rather weak.Next consider Example3.4 which uses a different geometry than that in the reference case, with thicker region Y3. This produces a slightly smaller c since the gradient of U is smaller. Also, ϕ~1 and k~1 are predictably different. The error quantity is influenced insignificantly.The opposite effect is seen in Table4 when Y3 is selected as thinner than in the reference Example 3.2.


### 3.2.4. Coarse Macrogrid and Parameter Grid

In the examples presented in Section3.2.3 we used h=h1=h2. This choice of compatible grids is convenient for visualization purposes. However, for the idea of the macro-model to be useful, the computational complexity of the macro-model needs to be lower than that of the micro-model. Since the macro-model is a system of two equations, its lower complexity can be only achieved if the grid in the macro-model is significantly coarser than that of micro-model. Below we show that one can use h≪h1 in the macro-model.Example 3.6.
We use data from Example3.2 with ε=0.2. Consider a fixed microgrid with h-1=100 and a few cases of macrogrid with h1-1=h2-1=10,50,100. Now we compare the micro- and macrosolutions ρε,hn and ρ1,h1, see Table 5. We also compare the run time of the macro-model to that of the reference micro-model which is 120 seconds.Table 5
Comparison of the micro-model to the macro-model as in Example3.6.
h
1
-
1
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
run time(s)
100
7.0915
×
1
0
-
3
1.9044
×
1
0
-
2
511
50
7.0937
×
1
0
-
3
1.9048
×
1
0
-
2
56
10
7.2380
×
1
0
-
3
1.9173
×
1
0
-
2
3We see that, as predicted above, the macro-model run time athi=h is not competitive with that of the micro-model; this is exacerbated by the overhead of various adaptive procedures to be discussed below. However, when hi≫h, the ability of the macro-model solutions to approximate those of micro-model appears reasonable and is associated with a much lower cost. This is true even when hi is by a factor of 10 coarser than h resulting in a computational time decreased by two orders of magnitude.Clearly, the macro-model can run much faster than the micro-model especially when adaptive nonuniform grids withh1≠h2 are implemented, see Section 4.Influence of Parameter Grid
Last we address the effect of the computed parametersc, k~1 on the closeness of the macro-model solution to that of micro-model. These coefficients are precomputed as discussed in Section 3.1.1.Example 3.7.
Letε=1/3, and let the micro-model grid be h-1=60. Use h1-1=h2-1=60. Compute numerical parameters chc, k~1,hk1 by approximating solutions to (3.4a), (3.4b), (3.4c), and (3.4d)  and (3.3a), (3.3b), and (3.3c), respectively.
Example3.7 and Table 6 show that the grid used for computing these parameters only mildly affects the quality of the solution, thus one can use coarse grid with confidence.Table 6
Sensitivity of the macro-model solution to the grid used to computechc≈c and k~1,hk1≈k~1 as Example 3.7.
h
c
-
1
c
h
c
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
10
1.9000
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
50
1.8370
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
100
1.8292
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
200
1.8258
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
h
k
-
1
k
~
1
(
1,1
)
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
10
4.6306
×
1
0
-
1
3.5392
×
1
0
-
2
5.7835
×
1
0
-
2
50
4.5141
×
1
0
-
1
3.6222
×
1
0
-
2
5.8989
×
1
0
-
2
100
4.5046
×
1
0
-
1
3.6301
×
1
0
-
2
5.9093
×
1
0
-
2
200
4.5009
×
1
0
-
1
3.6332
×
1
0
-
2
5.9134
×
1
0
-
2


## 3.2.1. Setup of Computational Experiments

We set up simulations for comparison using compatible data for micro- and macro-models; we setf≡0 and choose initial and boundary data driving the flow with interesting dynamics. For porosities and conductivities, we use the values computed in Section 3.1.2.LetΩ=(0,1)2 and ΓD={0}×[0,1]∪{1}×[0,1]. On ΓD we define ρD(0,y)=1,ρD(1,y)=0,y∈[0,1]. On the lateral sides of Ω the Neumann no-flow condition is imposed. Also, let ρ0≡1. Thus the flow in the micro-model goes from left to right, and the solution evolves towards the stationary solution 〈1-x,1-x〉. Due to the incompatibility between initial and boundary data, we have high gradients of the solution close to x=1 for small t.We fix the final timeT=0.05 and use uniform time stepping with Δt=10-4. The solution of the micro-model ρ depends on the number of cells in Ω≡Ωε, with ε∈(0,1]. For example, if ε=1/2, the domain Ω is composed of 2×2=4 cells. We denote the numerical solution of (2.8a) and (2.8b) at tn by ρε,hn where h denotes the grid parameter for the micro-model.For the macro-model, we useρ~i,0=ρ0≡1. Also, we use ρ~D,i=ρD for i=1,2 when k~2≢0. If k~2≡0, then the boundary condition for ρ~2 is not prescribed.


## 3.2.2. Qualitative Comparison

In our first comparison we use parameters from Example3.1. We solve the macro-model (4.2) for {〈ρ~1,h1n,ρ~2,h2n〉}n=1N and the micro-model (2.8a) and (2.8b) for ρε,hn. The plots of ρhn and ρ~i,hn are in Figure 2, and a different view is shown in Figure 1(d).Numerical approximation attn=3×10-2 of (a)-(b) the solutions of the two components of the macro-model (4.2), and (c) the solutions of the micro-model (2.8a) and (2.8b). Figure (d) shows the side view of (a, b, c) at y=0.4 and y=0.5. Data from Example 3.1 with 1/ε=5 and h=h1=h2=1/50. Note that (a) and (b) have different vertical scales.
(a)
(b)
(c)
(d)The heterogeneous structure ofΩε is well visible from the behavior of the micro-model solution ρhn. Large gradients of solution are visible on cell boundaries due to the large difference between k1, k2, and k3.However, the behavior in the fast region is well approximated byρ~1,h1n which envelopes ρhnχ1, while ρ~2,h2n envelopes ρhnχ2 well. This is very well seen in a side view in Figure 2(d).Quite interesting is the behavior ofρ2|x=1. Since in our examples k~2≡0, no boundary condition at x=1 is imposed on ρ2. Thus, for small t, ρ2(1,t) evolves from its initial constant value of 1 according to (1.3b) with the input from ρ~1. The latter satisfies, however, the homogeneous boundary condition at x=1. In particular, with constant ϕ~2, c, we have
(3.11)ρ~2(x,y,t)=ρ~2,0(x,y)e-(c/ϕ~2)t+∫0t‍ρ~1(x,y,s)e-(c/ϕ~2)(t-s)ds.
Thus, for small t and small c/ϕ~2, ρ~2(1,y,t) is away from 0, but, as time increases, its magnitude decreases proportionally to e-(c/ϕ~2)t.Next, we use parameters from Example3.2 and plot them in Figure 3. In this example the conductivity k1 is much larger than that of previous case, the local gradients in the micro-model are smaller, and the solutions to the macro-model achieve a smoother profile faster.Numerical approximation attn=3×10-2 of (a) and (b) the solutions of the two components of the macro-model (4.2) and (c) the solutions of the micro-model (2.8a) and (2.8b). Figure (d) shows the side view of (a, b, c) at y=0.4 and y=0.5. Data from Example 3.2 with 1/ε=5 and h=h1=h2=1/50.
(a)
(b)
(c)
(d)


## 3.2.3. Quantitative Comparison—Fine Grid in Macro-Model

Now we are ready to discuss a quantitative comparison between solutions to the micro- and macro-models. For this, we set up a family of cases withε=1/(2l+1),  l≥1. We use uniform mesh with h=h1=h2=ε/20 for now. See Section 3.2.4 for hi≫h and Section 4 for adaptive gridding.From  [11, Theorem  18, page 208]  and from our plots it follows that we can actually compare the micro- and macrosolutions, as long as the characteristic functions χi are involved. While the results in [11] use two-scale convergence with ε, that is, rely on special periodic test functions, we compare ρε,hn and ρ~i,hn directly. We see that it is easy to do so only in the connected region Ω1.We define a quantityeεn=(ρε,hn-ρ~1,h1n)χ1 to be used in comparison. We also define, for a fixed t∈(0,T] and some function z,
(3.12)‖z(x,y;t)‖*,L1:=∫01‍|z(x,0.5;t)|dx,‖z(x,y;t)‖*,L2:={∫01‍|z(x,0.5;t)|2dx}1/2,
which will be applied to eεn. Clearly ∥·∥*,W is not a norm but is a useful quantity of interest.In Tables1, 2, 3, and 4 we present ∥eεn∥*,W for different ε and W=L2,L1. For each case we consider convergence of ∥eεn∥*,W with ε↓0.Table 1
Comparison of solutions to the macro-model with those to the micro-model for different values ofε with data as in Example 3.2.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1832
×
1
0
-
2
—
3.2071
×
1
0
-
2
—
5
×
5
1
/
5
7.0915
×
1
0
-
3
1.0188
1.9044
×
1
0
-
2
1.0390
7
×
7
1
/
7
5.0522
×
1
0
-
3
1.0077
1.3528
×
1
0
-
2
1.0164
9
×
9
1
/
9
3.9285
×
1
0
-
3
1.0011
1.0502
×
1
0
-
2
1.0075
11
×
11
1
/
11
3.2162
×
1
0
-
3
0.9969
8.5890
×
1
0
-
3
1.0021Table 2
Comparison of the macro-model to the micro-model for different values ofε as in Example 3.3.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1446
×
1
0
-
2
—
3.0715
×
1
0
-
2
—
5
×
5
1
/
5
7.0065
×
1
0
-
3
1.0408
1.8387
×
1
0
-
2
0.9956
7
×
7
1
/
7
5.1541
×
1
0
-
3
1.0958
1.3363
×
1
0
-
2
1.0543
9
×
9
1
/
9
4.1704
×
1
0
-
3
1.1867
1.0684
×
1
0
-
2
1.1232
11
×
11
1
/
11
3.5794
×
1
0
-
3
1.3131
9.0460
×
1
0
-
3
1.2058Table 3
Comparison of the macro-model to the micro-model for different values ofε Example 3.4.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.3758
×
1
0
-
2
—
3.5390
×
1
0
-
2
—
5
×
5
1
/
5
8.2810
×
1
0
-
3
1.0330
2.0998
×
1
0
-
2
0.9786
7
×
7
1
/
7
5.7341
×
1
0
-
3
0.9155
1.4741
×
1
0
-
2
0.9510
9
×
9
1
/
9
4.3987
×
1
0
-
3
0.9479
1.1336
×
1
0
-
2
0.9568Table 4
Comparison of the macro-model to the micro-model for different values ofε Example 3.5.
# cells
ε
∥
e
ε
N
∥
*
,
L
1
α
1
∥
e
ε
N
∥
*
,
L
2
α
2
3
×
3
1
/
3
1.1826
×
1
0
-
2
—
3.2048
×
1
0
-
2
—
5
×
5
1
/
5
7.0877
×
1
0
-
3
0.9978
1.9029
×
1
0
-
2
0.9800
7
×
7
1
/
7
5.0497
×
1
0
-
3
0.9924
1.3518
×
1
0
-
2
0.9840
9
×
9
1
/
9
3.9281
×
1
0
-
3
1.0006
1.0502
×
1
0
-
2
0.9955From homogenization theory it is not clear what order of convergence should be expected. We observe that∥eεn∥*,W decreases linearly with ε for both *=L1,L2 -like quantities in all cases.Next we discuss the dependence of the results on conductivity and other data. Consider Example3.2 as a reference example with its corresponding Table 1. Compare it now with that from Table 2 where k3 is larger than that in the reference case to see how this influences the errors. As expected, and as shown in Example 3.3, c changes by the same factor, and this means that the coupling between ρ~1 and ρ~2 is stronger which corresponds to faster flow across Y3 in the micro-model. However, the influence on the convergence of the micro-macro difference, from Table 2 is rather weak.Next consider Example3.4 which uses a different geometry than that in the reference case, with thicker region Y3. This produces a slightly smaller c since the gradient of U is smaller. Also, ϕ~1 and k~1 are predictably different. The error quantity is influenced insignificantly.The opposite effect is seen in Table4 when Y3 is selected as thinner than in the reference Example 3.2.


## 3.2.4. Coarse Macrogrid and Parameter Grid

In the examples presented in Section3.2.3 we used h=h1=h2. This choice of compatible grids is convenient for visualization purposes. However, for the idea of the macro-model to be useful, the computational complexity of the macro-model needs to be lower than that of the micro-model. Since the macro-model is a system of two equations, its lower complexity can be only achieved if the grid in the macro-model is significantly coarser than that of micro-model. Below we show that one can use h≪h1 in the macro-model.Example 3.6.
We use data from Example3.2 with ε=0.2. Consider a fixed microgrid with h-1=100 and a few cases of macrogrid with h1-1=h2-1=10,50,100. Now we compare the micro- and macrosolutions ρε,hn and ρ1,h1, see Table 5. We also compare the run time of the macro-model to that of the reference micro-model which is 120 seconds.Table 5
Comparison of the micro-model to the macro-model as in Example3.6.
h
1
-
1
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
run time(s)
100
7.0915
×
1
0
-
3
1.9044
×
1
0
-
2
511
50
7.0937
×
1
0
-
3
1.9048
×
1
0
-
2
56
10
7.2380
×
1
0
-
3
1.9173
×
1
0
-
2
3We see that, as predicted above, the macro-model run time athi=h is not competitive with that of the micro-model; this is exacerbated by the overhead of various adaptive procedures to be discussed below. However, when hi≫h, the ability of the macro-model solutions to approximate those of micro-model appears reasonable and is associated with a much lower cost. This is true even when hi is by a factor of 10 coarser than h resulting in a computational time decreased by two orders of magnitude.Clearly, the macro-model can run much faster than the micro-model especially when adaptive nonuniform grids withh1≠h2 are implemented, see Section 4.Influence of Parameter Grid
Last we address the effect of the computed parametersc, k~1 on the closeness of the macro-model solution to that of micro-model. These coefficients are precomputed as discussed in Section 3.1.1.Example 3.7.
Letε=1/3, and let the micro-model grid be h-1=60. Use h1-1=h2-1=60. Compute numerical parameters chc, k~1,hk1 by approximating solutions to (3.4a), (3.4b), (3.4c), and (3.4d)  and (3.3a), (3.3b), and (3.3c), respectively.
Example3.7 and Table 6 show that the grid used for computing these parameters only mildly affects the quality of the solution, thus one can use coarse grid with confidence.Table 6
Sensitivity of the macro-model solution to the grid used to computechc≈c and k~1,hk1≈k~1 as Example 3.7.
h
c
-
1
c
h
c
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
10
1.9000
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
50
1.8370
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
100
1.8292
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
200
1.8258
×
1
0
-
3
3.6251
×
1
0
-
2
5.9026
×
1
0
-
2
h
k
-
1
k
~
1
(
1,1
)
∥
e
ε
N
∥
*
,
L
1
∥
e
ε
N
∥
*
,
L
2
10
4.6306
×
1
0
-
1
3.5392
×
1
0
-
2
5.7835
×
1
0
-
2
50
4.5141
×
1
0
-
1
3.6222
×
1
0
-
2
5.8989
×
1
0
-
2
100
4.5046
×
1
0
-
1
3.6301
×
1
0
-
2
5.9093
×
1
0
-
2
200
4.5009
×
1
0
-
1
3.6332
×
1
0
-
2
5.9134
×
1
0
-
2


## 4. Error Estimates and Two-Grid Adaptivity

We present now the details of discrete formulation of (2.4) using independent grids for each component and adaptive gridding. The adaptive two-grid algorithm further reduces the computational cost of the macro-model. The formalism of the two-grid solution connects loosely to [37, 38] where their benefit was considered for a nonlinear scalar equation.In what follows we assume thatk~1≫k~2, that is, that the first component varies in space faster than the second. As a special case, this includes the case k~2=0. We consider two triangulations 𝒯hi, i=1,2, which will be used for ρ~hi, respectively. In principle, these can be chosen independently. In fact, we allow the triangulations to vary in time and thus consider 𝒯hin and the associated spaces Vhin as in (2.7). To avoid the loss of accuracy due to excessive interpolation and intergrid projections and because k~1≫k~2, we assume that 𝒯h1n is a refinement of the partition 𝒯h2n.We need intergrid operators to handle two components that live on separate grids. LetΠ:Vh2n→Vh1n be the interpolation operator and Π′:Vh1n→Vh2n the L2 projection defined by
(4.1)(Π′ϕh1,ψh2):=(ϕh1,Πψh2),∀ψh2∈Vh2.Now we define the discrete solutions. At each time step1≤n≤N, we seek ρ~i,hin∈Vhin, i=1,2, which satisfy the discrete problem, that is, (2.4) restricted to the finite dimensional subspaces Vhin so that, for all ψh1∈Vh1n, for all ξh2∈Vh2n(4.2)(ϕ~1ρ~1,h1n-ρ~1,h1n-1Δtn,ψh1)+(k~1∇ρ~1,h1n,∇ψh1)+(c(ρ~1,h1n-Πρ~2,h2n),ψh1)=(f,ψh1),(ϕ~2ρ~2,h2n-ρ~2,h2n-1Δtn,ξh2)+(k~2∇ρ~2,h2n,∇ξh2)+(c(ρ~2,h2n-Π′ρ~1,h1n),ξh2)=0,〈ρ~1,h10,ρ~2,h20〉=〈Ih1ρ~1,0,Ih2ρ~2,0〉inΩ.
If k~2=0, the system (4.2) is modified appropriately so that ρ~2(·,t)∈L2(Ω) instead of HD1(Ω), and we do not impose boundary conditions on this component. In fact, instead of (2.7) we define
(4.3)Vh2n={vh∈C(Ω¯):∀T∈Thn,vh∣T∈Pk(T)}.
(It holds that vh∈C(Ω¯) but it is not necessary). It is straightforward that (4.2) is uniquely solvable.


### 4.1.A Priori Error Estimate

We have the following convergence result proved for the error in energy norm. In what follows1≤n≤N and i=1,2. We denote ei,hi(tn)=ρ~i,hin-ρ~in and define the energy norm for the product space
(4.4)|‖〈z,w〉‖|2(tn):=ϕ~1‖z(tn)‖2+ϕ~2‖w(tn)‖2+∑m=1n‍Δtm(k~1|z(tm)|12+k~2|w(tm)|12).Assuming that〈ρ~1,ρ~2〉 is sufficiently smooth we have the following a-priori estimate proven in the appendix.Theorem 4.1.
Let〈ρ~1,ρ~2〉,{〈ρ~1,h1n,ρ~2,h2n〉}n=1N be the smooth solutions of (2.4) and (4.2), respectively. Then one has
(4.5)|‖〈e1,h1,e2,h2〉‖|(tn)≤k~11/2C1h1+k~21/2C2h2+C3max⁡m=1:nΔtm+higher  order  terms(h.o.t.),
where C1, C2, and C3 are independent of h.
Furthermore, ifk~2=0, then
(4.6)|‖〈e1,h1,e2,h2〉‖|(tn)≤k~11/2C1h1+C2h22+C3max⁡m=1:nΔtm+higher  order  terms(h.o.t.).Now we notice that (4.5) suggests the following choice of grid parameters h1, h2. If k~1≫k~2, then to balance the components of the error one can use h1≪h2. Similarly, if k~2=0, then (4.6) implies that h2 can be chosen to be of the order of h11/2. The use of coarse grid h2≫h1 for the second component reduces the size of the linear system to be solved and decreases the computational time of the macro-model.Also, the a-priori estimates (4.5), (4.6) are global. One can further reduce the computational cost while maintaining accuracy by using local grid adaptivity guided by a posteriori estimators. This is pursued below.


### 4.2.A-Posteriori Error Estimates

Estimators for time-dependent problems can be defined in many ways including the now classical space-time element and adjoint approaches [27, 28]. In this paper we follow the residual estimator framework extending [15, 29] to the double-diffusion system using the ideas in [13, 39, 40] originally formulated for an elliptic system.The estimatorηn=(∑m=1n‍[ΔtmSm2+Tm2])1/2 for (4.2) is composed of the temporal part Tn which adapts the time stepping and the spatial part Sn which guides the spatial grid adaptivity. We define
(4.7)Tn2≔Δtn3(‖k~11/2(ρ~1,h1n-ρ~1,h1n-1)‖2+‖k~21/2(ρ~2,h2n-ρ~2,h2n-1)‖2)+Δtn3‖c1/2(ρ~1,h1n-ρ~1,h1n-1)-c1/2(ρ~2,h2n-ρ~2,h2n-1)‖2,
and Sn2=S1,n+S2,n, where Si,n=∑Ti∈𝒯hi‍Si,Tin, and where, as usual (see [23]), one defines
(4.8)Si,Tin:=θi,Ti2‖Ri,Tin‖2+12∑E∈ETi‍γi,E2‖Ri,En‖2.
The element and edge residuals in Si,Tin,i=1,2, are given by
(4.9)Ri,Tin=fi-ϕ~iρ~i,hin-ρ~i,hin-1Δtn+∇⋅(k~i∇ρ~i,hin)+(-1)ic(ρ~1,h1n-ρ~2,h2n),(4.10)Ri,En={[k~i∂νρ~i,hin]EifE⊄∂Ω,-k~i∂νρ~i,hinifE⊂ΓN,0,otherwise.
Here [z]E denotes the jump of the flux z across the edge E of an element. Usually f1≡f, f2≡0. Also, if k~2=0, then R2,En≡0.The scaling constantsθi,T,γi,E take into account the contribution of the exchange term c(ρ~1-ρ~2) and are defined so that the estimators remain robust when the coefficients of the problem change substantially:
(4.11)θ1,T=min⁡{hT1k~1-1/2,max⁡{c-1/2,hT1k~2-1/2}},T1∈Th1n∪Eh1n,(4.12)θ2,T=min⁡{hT2k~2-1/2,max⁡{c-1/2,hT2k~1-1/2}},T2∈Th2n∪Eh2n,(4.13)γi,E=2hE-1/2θi,E,i=1,2.
This a posteriori estimator works well with two-grid discretizations as shown in the examples to follow. It is well formulated also for k~2≡0 when ρ~2∈L2(Ω). We have the following result on reliability of the estimator proven in the appendix.Theorem 4.2.
Let〈ρ~1,ρ~2〉, {〈ρ~1,h1n,ρ~2,h2n〉}n=1N be the solution of (1.2a) and (1.2b), (2.3) and (4.2), respectively. Then
(4.14)|‖〈e1,h1,e2,h2〉‖|(tn)≤ηn+h.o.t.;∀n≥1


### 4.3. Adaptive Two-Grid Discretization and Implementation

Consider a fixed time stepn=no for which some triangulation is chosen, the solution is found, and the spatial error indicators Si,n are computed. To apply the two-grid spatial adaptivity we use the following refinement algorithm.Adaptive Two-Grid AlgorithmA
(1)
Select the triangulation(s) to be refined: IfS1,no>3S2,no, refine only 𝒯h1no. If S2,no>3S1,no, refine only 𝒯h2no. Otherwise, refine both.
(2)
In𝒯hi, i=1,2 selected in Step 1, refine any T∈𝒯h,i for which Si,Tno≥0.5maxT∈𝒯hiSi,Tno.
(3)
Enforce the requirement that𝒯h1n is a refinement of 𝒯h2n by adding extra elements as needed.
The steps in𝒜 should be repeated a certain number N𝒜 of times until an ideal grid is found. See [41] for analysis of whether the iterative process of refinement and coarsening is, in general, convergent.Example of Two-Grid Adaptivity
Now we show an example on how this adaptive algorithm works for the double-diffusion system. We use data and setup from Example3.2, with T=0.02 and uniform time stepping with Δt=10-4. Other data is as in Section 3.2.1.
Consider a fixedn0=200 with a uniform triangulation 𝒯h1no=𝒯h2no with h1=h2=0.1. Then apply NA=6 times the steps in 𝒜, see the solution and the refined meshes in Figure 4. The process works as expected: since at the time tn0 the solution has a high gradient near x=1, the refinement occurs there. The refinement affects the grid for the first component only, because the gradients of the second component are not included due to k~2=0. This can be compared to the a-priori estimate in Theorem 4.1 from which we know the convergence is O(h1+h22+Δtn).
In Table7 we provide the details on 𝒯hino and compare the effectiveness of the local two-grid refinement by algorithm 𝒜 with uniform grid refinement. With adaptive refinement we get Sno=0.35834 with 435+121=556 unknowns, while the uniform refinement needs 441+441=882 unknowns to achieve comparable Sno=0.39221.Table 7
Adaptive versus uniform refinement.
Uniform refinement
S
n
o
# elem𝒯h1no
# elem𝒯h1no
dim
(
V
h
1
n
o
)
dim
(
V
h
2
n
o
)
Initial mesh
0.81803
200
200
121
121
1st refinement
0.39221
800
800
441
441
Adaptive refinement
S
n
o
# elem𝒯h1no
# elem𝒯h2no
dim
(
V
h
1
n
o
)
dim
(
V
h
1
n
o
)
Initial mesh
0.8180
200
200
121
121
1st level
0.8031
240
200
141
121
2nd level
0.7076
292
200
167
121
3rd level
0.6499
316
200
179
121
4th level
0.5525
462
200
254
121
5th level
0.4205
720
200
385
121
6th level
0.3583
820
200
435
121Sixth level refinement iteration attno=2×10-2. Figures (a) and (b): plots of ρ~1,h1no and 𝒯h1no. Figures (c), (e), and (d): plots of ρ~2,h2no and 𝒯h2no.
(a)
(b)
(c)
(d)Last, we describe the process with which𝒯hin,  i=1,2, may vary between time steps. For a fixed n denote by 𝒯hin the final triangulation obtained by the adaptive algorithm 𝒜. For the new time step n+1, we use 𝒯hin as the initial triangulation. If the algorithm 𝒜 suggests that 𝒯hin+1 is modified, we need to project ρ~i,hin to the new grid so it can be used as initial data for the new step, but ρ~i,hin differs from their interpolation only in the elements where the triangulation is coarsened. Another way to deal with different triangulations 𝒯hin+1, 𝒯hin is to compute the a posteriori error estimator in a common refinement of the two triangulations 𝒯hin+1, 𝒯hin [29, 42], but this will not be pursued further.Implementation
The decision to use two grids implies that we have to interpolate between the finite element spacesVh1n and Vh2n, see (4.1). Suppose that {ζi,j}j,  i=1,2, is a basis for Vhin. Then the operator Π in (4.1) is represented by the matrix 𝕀h1,h2 whose components are given by 𝕀h1,h2(j,k)=(ζ1,j,ζ2,k).
Consider now and compare the computational cost associated with one-grid and two-grid approaches. To solve (4.2) with one-grid Vh1n we need to assemble the stiffness matrix S1,h1(j,k)=(k~1∇ζ1,j,∇ζ1,k), the mass matrix M1,h1(j,k)=(ζ1,j,ζ1,k), and the block matrix
(4.15)Bone-grid=[-ΔtS1,h1+(ϕ~1+cΔt)M1,h1-cΔtM1,h1-cΔtM1,h1(ϕ~2+cΔt)M1,h1].
At each time step we solve the linear system with Bone-grid; this is referred to as 𝕊one-grid.
To solve (4.2) with two grids and two spaces Vh1n,Vh2n, we need to assemble the stiffness matrix S1,h1(j,k)=(k~1∇ζ1,j,∇ζ1,k), the mass matrices M1,h1(j,k)=(ζ1,j,ζ1,k), M2,h2(j,k)=(ζ2,j,ζ2,k), the interpolation matrix 𝕀h1,h2(j,k)=(ζ1,j,ζ2,k), and the block matrix
(4.16)Btwo-grid=[-ΔtS1,h1+(ϕ~1+cΔt)M1,h1-cΔtIh1,h2T-cΔtIh1,h2(ϕ~2+cΔt)M2,h2].
To finish we solve the linear system which we refer to as 𝕊two-grid.Suppose now we solve (4.2) and wish to maintain Sn≤δ for some tolerance δ over 1,000 time steps. With data as in Section 3.1.2 we run the simulations using one grid and find that we need h1=h2=1/80. Using the two-grid approach we only need h1=1/80,h2=1/20.The computational time needed to assemble the matrices and solve the systems in our MATLAB implementation is displayed in Table8. Clearly the assembly process takes more time for two-grid case than for one-grid discretization. However, the economy in the cost of solving the linear system makes up for that cost.Table 8
Cost of one-grid versus two-grid discretization in a simulation with 1,000 time steps.
One-grid withh1=h2=1/80
Matrix
Cost (seconds)
S
1
,
h
1 + M1,h1
3.6547
B
one-grid
0.0123
Solve𝕊one-grid 1,000 times
0.2761× 1,000
Total
279.7670
Two-grid withh1=1/80,h2=1/20
Matrix
Cost (seconds)
S
1
,
h
1 + M1,h1
3.6547
M
2
,
h
2
0.0438
𝕀
h
1
,
h
2
0.2121
B
two-grid
0.0063
Solve𝕊two-grid 1,000 times
0.1030× 1,000
Total
106.9169


## 4.1.A Priori Error Estimate

We have the following convergence result proved for the error in energy norm. In what follows1≤n≤N and i=1,2. We denote ei,hi(tn)=ρ~i,hin-ρ~in and define the energy norm for the product space
(4.4)|‖〈z,w〉‖|2(tn):=ϕ~1‖z(tn)‖2+ϕ~2‖w(tn)‖2+∑m=1n‍Δtm(k~1|z(tm)|12+k~2|w(tm)|12).Assuming that〈ρ~1,ρ~2〉 is sufficiently smooth we have the following a-priori estimate proven in the appendix.Theorem 4.1.
Let〈ρ~1,ρ~2〉,{〈ρ~1,h1n,ρ~2,h2n〉}n=1N be the smooth solutions of (2.4) and (4.2), respectively. Then one has
(4.5)|‖〈e1,h1,e2,h2〉‖|(tn)≤k~11/2C1h1+k~21/2C2h2+C3max⁡m=1:nΔtm+higher  order  terms(h.o.t.),
where C1, C2, and C3 are independent of h.
Furthermore, ifk~2=0, then
(4.6)|‖〈e1,h1,e2,h2〉‖|(tn)≤k~11/2C1h1+C2h22+C3max⁡m=1:nΔtm+higher  order  terms(h.o.t.).Now we notice that (4.5) suggests the following choice of grid parameters h1, h2. If k~1≫k~2, then to balance the components of the error one can use h1≪h2. Similarly, if k~2=0, then (4.6) implies that h2 can be chosen to be of the order of h11/2. The use of coarse grid h2≫h1 for the second component reduces the size of the linear system to be solved and decreases the computational time of the macro-model.Also, the a-priori estimates (4.5), (4.6) are global. One can further reduce the computational cost while maintaining accuracy by using local grid adaptivity guided by a posteriori estimators. This is pursued below.


## 4.2.A-Posteriori Error Estimates

Estimators for time-dependent problems can be defined in many ways including the now classical space-time element and adjoint approaches [27, 28]. In this paper we follow the residual estimator framework extending [15, 29] to the double-diffusion system using the ideas in [13, 39, 40] originally formulated for an elliptic system.The estimatorηn=(∑m=1n‍[ΔtmSm2+Tm2])1/2 for (4.2) is composed of the temporal part Tn which adapts the time stepping and the spatial part Sn which guides the spatial grid adaptivity. We define
(4.7)Tn2≔Δtn3(‖k~11/2(ρ~1,h1n-ρ~1,h1n-1)‖2+‖k~21/2(ρ~2,h2n-ρ~2,h2n-1)‖2)+Δtn3‖c1/2(ρ~1,h1n-ρ~1,h1n-1)-c1/2(ρ~2,h2n-ρ~2,h2n-1)‖2,
and Sn2=S1,n+S2,n, where Si,n=∑Ti∈𝒯hi‍Si,Tin, and where, as usual (see [23]), one defines
(4.8)Si,Tin:=θi,Ti2‖Ri,Tin‖2+12∑E∈ETi‍γi,E2‖Ri,En‖2.
The element and edge residuals in Si,Tin,i=1,2, are given by
(4.9)Ri,Tin=fi-ϕ~iρ~i,hin-ρ~i,hin-1Δtn+∇⋅(k~i∇ρ~i,hin)+(-1)ic(ρ~1,h1n-ρ~2,h2n),(4.10)Ri,En={[k~i∂νρ~i,hin]EifE⊄∂Ω,-k~i∂νρ~i,hinifE⊂ΓN,0,otherwise.
Here [z]E denotes the jump of the flux z across the edge E of an element. Usually f1≡f, f2≡0. Also, if k~2=0, then R2,En≡0.The scaling constantsθi,T,γi,E take into account the contribution of the exchange term c(ρ~1-ρ~2) and are defined so that the estimators remain robust when the coefficients of the problem change substantially:
(4.11)θ1,T=min⁡{hT1k~1-1/2,max⁡{c-1/2,hT1k~2-1/2}},T1∈Th1n∪Eh1n,(4.12)θ2,T=min⁡{hT2k~2-1/2,max⁡{c-1/2,hT2k~1-1/2}},T2∈Th2n∪Eh2n,(4.13)γi,E=2hE-1/2θi,E,i=1,2.
This a posteriori estimator works well with two-grid discretizations as shown in the examples to follow. It is well formulated also for k~2≡0 when ρ~2∈L2(Ω). We have the following result on reliability of the estimator proven in the appendix.Theorem 4.2.
Let〈ρ~1,ρ~2〉, {〈ρ~1,h1n,ρ~2,h2n〉}n=1N be the solution of (1.2a) and (1.2b), (2.3) and (4.2), respectively. Then
(4.14)|‖〈e1,h1,e2,h2〉‖|(tn)≤ηn+h.o.t.;∀n≥1


## 4.3. Adaptive Two-Grid Discretization and Implementation

Consider a fixed time stepn=no for which some triangulation is chosen, the solution is found, and the spatial error indicators Si,n are computed. To apply the two-grid spatial adaptivity we use the following refinement algorithm.Adaptive Two-Grid AlgorithmA
(1)
Select the triangulation(s) to be refined: IfS1,no>3S2,no, refine only 𝒯h1no. If S2,no>3S1,no, refine only 𝒯h2no. Otherwise, refine both.
(2)
In𝒯hi, i=1,2 selected in Step 1, refine any T∈𝒯h,i for which Si,Tno≥0.5maxT∈𝒯hiSi,Tno.
(3)
Enforce the requirement that𝒯h1n is a refinement of 𝒯h2n by adding extra elements as needed.
The steps in𝒜 should be repeated a certain number N𝒜 of times until an ideal grid is found. See [41] for analysis of whether the iterative process of refinement and coarsening is, in general, convergent.Example of Two-Grid Adaptivity
Now we show an example on how this adaptive algorithm works for the double-diffusion system. We use data and setup from Example3.2, with T=0.02 and uniform time stepping with Δt=10-4. Other data is as in Section 3.2.1.
Consider a fixedn0=200 with a uniform triangulation 𝒯h1no=𝒯h2no with h1=h2=0.1. Then apply NA=6 times the steps in 𝒜, see the solution and the refined meshes in Figure 4. The process works as expected: since at the time tn0 the solution has a high gradient near x=1, the refinement occurs there. The refinement affects the grid for the first component only, because the gradients of the second component are not included due to k~2=0. This can be compared to the a-priori estimate in Theorem 4.1 from which we know the convergence is O(h1+h22+Δtn).
In Table7 we provide the details on 𝒯hino and compare the effectiveness of the local two-grid refinement by algorithm 𝒜 with uniform grid refinement. With adaptive refinement we get Sno=0.35834 with 435+121=556 unknowns, while the uniform refinement needs 441+441=882 unknowns to achieve comparable Sno=0.39221.Table 7
Adaptive versus uniform refinement.
Uniform refinement
S
n
o
# elem𝒯h1no
# elem𝒯h1no
dim
(
V
h
1
n
o
)
dim
(
V
h
2
n
o
)
Initial mesh
0.81803
200
200
121
121
1st refinement
0.39221
800
800
441
441
Adaptive refinement
S
n
o
# elem𝒯h1no
# elem𝒯h2no
dim
(
V
h
1
n
o
)
dim
(
V
h
1
n
o
)
Initial mesh
0.8180
200
200
121
121
1st level
0.8031
240
200
141
121
2nd level
0.7076
292
200
167
121
3rd level
0.6499
316
200
179
121
4th level
0.5525
462
200
254
121
5th level
0.4205
720
200
385
121
6th level
0.3583
820
200
435
121Sixth level refinement iteration attno=2×10-2. Figures (a) and (b): plots of ρ~1,h1no and 𝒯h1no. Figures (c), (e), and (d): plots of ρ~2,h2no and 𝒯h2no.
(a)
(b)
(c)
(d)Last, we describe the process with which𝒯hin,  i=1,2, may vary between time steps. For a fixed n denote by 𝒯hin the final triangulation obtained by the adaptive algorithm 𝒜. For the new time step n+1, we use 𝒯hin as the initial triangulation. If the algorithm 𝒜 suggests that 𝒯hin+1 is modified, we need to project ρ~i,hin to the new grid so it can be used as initial data for the new step, but ρ~i,hin differs from their interpolation only in the elements where the triangulation is coarsened. Another way to deal with different triangulations 𝒯hin+1, 𝒯hin is to compute the a posteriori error estimator in a common refinement of the two triangulations 𝒯hin+1, 𝒯hin [29, 42], but this will not be pursued further.Implementation
The decision to use two grids implies that we have to interpolate between the finite element spacesVh1n and Vh2n, see (4.1). Suppose that {ζi,j}j,  i=1,2, is a basis for Vhin. Then the operator Π in (4.1) is represented by the matrix 𝕀h1,h2 whose components are given by 𝕀h1,h2(j,k)=(ζ1,j,ζ2,k).
Consider now and compare the computational cost associated with one-grid and two-grid approaches. To solve (4.2) with one-grid Vh1n we need to assemble the stiffness matrix S1,h1(j,k)=(k~1∇ζ1,j,∇ζ1,k), the mass matrix M1,h1(j,k)=(ζ1,j,ζ1,k), and the block matrix
(4.15)Bone-grid=[-ΔtS1,h1+(ϕ~1+cΔt)M1,h1-cΔtM1,h1-cΔtM1,h1(ϕ~2+cΔt)M1,h1].
At each time step we solve the linear system with Bone-grid; this is referred to as 𝕊one-grid.
To solve (4.2) with two grids and two spaces Vh1n,Vh2n, we need to assemble the stiffness matrix S1,h1(j,k)=(k~1∇ζ1,j,∇ζ1,k), the mass matrices M1,h1(j,k)=(ζ1,j,ζ1,k), M2,h2(j,k)=(ζ2,j,ζ2,k), the interpolation matrix 𝕀h1,h2(j,k)=(ζ1,j,ζ2,k), and the block matrix
(4.16)Btwo-grid=[-ΔtS1,h1+(ϕ~1+cΔt)M1,h1-cΔtIh1,h2T-cΔtIh1,h2(ϕ~2+cΔt)M2,h2].
To finish we solve the linear system which we refer to as 𝕊two-grid.Suppose now we solve (4.2) and wish to maintain Sn≤δ for some tolerance δ over 1,000 time steps. With data as in Section 3.1.2 we run the simulations using one grid and find that we need h1=h2=1/80. Using the two-grid approach we only need h1=1/80,h2=1/20.The computational time needed to assemble the matrices and solve the systems in our MATLAB implementation is displayed in Table8. Clearly the assembly process takes more time for two-grid case than for one-grid discretization. However, the economy in the cost of solving the linear system makes up for that cost.Table 8
Cost of one-grid versus two-grid discretization in a simulation with 1,000 time steps.
One-grid withh1=h2=1/80
Matrix
Cost (seconds)
S
1
,
h
1 + M1,h1
3.6547
B
one-grid
0.0123
Solve𝕊one-grid 1,000 times
0.2761× 1,000
Total
279.7670
Two-grid withh1=1/80,h2=1/20
Matrix
Cost (seconds)
S
1
,
h
1 + M1,h1
3.6547
M
2
,
h
2
0.0438
𝕀
h
1
,
h
2
0.2121
B
two-grid
0.0063
Solve𝕊two-grid 1,000 times
0.1030× 1,000
Total
106.9169


## 5. Conclusions

We compared the solutions to the micro-model to those of the macro-model and have shown that the latter is a good approximation to the former. This remains true also when a very coarse computational grid is used. Moreover, we established a linear convergence rate in function of the periodicity parameterε in a certain quantity of interest which further shows that the double-diffusion model is an excellent approximation to the micro-model, at least for the considered scenarios.To make the double-diffusion model computationally efficient, we proposed an algorithm for local grid adaptivity which allows each component to live on its own grid. The grid refinement is guided by an a posteriori error estimator for which we proved theoretical results.



---

*Source: 938727-2012-08-08.xml*