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.

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 Ξ©: πœ™πœ•πœŒπœ•π‘‘βˆ’βˆ‡β‹…(π‘˜βˆ‡πœŒ)=𝑓,inΞ©,𝑑>0,(1.1) where the coefficients πœ™, π‘˜ 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πœ•ΜƒπœŒ1ξ€·Μƒπ‘˜πœ•π‘‘βˆ’βˆ‡β‹…1βˆ‡ΜƒπœŒ1ξ€Έξ€·+π‘ΜƒπœŒ1βˆ’ΜƒπœŒ2ξ€Έ=ξ‚ξ‚πœ™π‘“,inΞ©,𝑑>0,(1.2a)2πœ•ΜƒπœŒ2ξ€·Μƒπ‘˜πœ•π‘‘βˆ’βˆ‡β‹…2βˆ‡ΜƒπœŒ2ξ€Έξ€·+π‘ΜƒπœŒ2βˆ’ΜƒπœŒ1ξ€Έ=0,inΞ©,𝑑>0,(1.2b)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 term 𝑐(ΜƒπœŒ2βˆ’ΜƒπœŒ1), and the coefficients Μƒπ‘˜π‘–, ξ‚πœ™π‘– 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 Μƒπ‘˜2=0 known as the Warren-Root model [6]ξ‚πœ™1πœ•ΜƒπœŒ1ξ€·Μƒπ‘˜πœ•π‘‘βˆ’βˆ‡β‹…1βˆ‡ΜƒπœŒ1ξ€Έξ€·+π‘ΜƒπœŒ1βˆ’ΜƒπœŒ2ξ€Έ=ξ‚ξ‚πœ™π‘“,inΞ©,𝑑>0,(1.3a)2πœ•ΜƒπœŒ2ξ€·πœ•π‘‘+π‘ΜƒπœŒ2βˆ’ΜƒπœŒ1ξ€Έ=0,inΞ©,𝑑>0,(1.3b)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 Μƒπ‘˜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≀2β‰ͺΜƒπ‘˜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 Μƒπ‘˜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 Section 2 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 Ξ©βŠ‚β„π‘‘, 𝑑=2,3, be a bounded polygonal region with piecewise smooth boundary πœ•Ξ©=Γ𝐷βˆͺΓ𝑁, where Γ𝐷, Γ𝑁 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 𝐿2(πœ”) and Sobolev spaces π»π‘˜(πœ”), π‘˜βˆˆβ„•. These are equipped with the usual seminorms |β‹…|π‘˜,πœ”, norms β€–β‹…β€–π‘˜,πœ”βˆΆ=β€–β‹…β€–π»π‘˜(πœ”) and the usual scalar product (s) (𝑓,𝑔)πœ”βˆΆ=(𝑓,𝑔)𝐿2(πœ”). We denote β€–π‘“β€–πœ”βˆΆ=(𝑓,𝑓)πœ”1/2. If πœ”=Ξ©, the subscript πœ” will be omitted. We also let 𝐻1𝐷(Ξ©)∢={π‘£βˆˆπ»1(Ξ©)βˆΆπ‘£|Γ𝐷=0} under the standard norm in 𝐻1(Ξ©), ‖𝑣‖1∢=(‖𝑣‖2+β€–βˆ‡π‘£β€–2)1/2.

For any 𝑇>0 and any vector space 𝑆, the space 𝐿2(0,𝑇;𝑆) consists of all square-integrable functions with values in 𝑆 such that ‖𝑣‖𝐿2(0,𝑇;𝑆)∫∢=(𝑇0‖𝑣‖2𝑆)1/2 is finite. The space 𝐢([0,𝑇];𝑆) 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 πœ™=πœ™(π‘₯) and π‘˜=π‘˜(π‘₯),π‘₯∈Ω denote the nonnegative coefficients of porosity and conductivities, respectively, and 𝑓 accounts for external sources. We assume for simplicity that π‘˜ 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 𝜌=𝜌0,inΩ×{0},𝜌=𝜌𝐷,onΓ𝐷,π‘˜βˆ‡πœŒβ‹…β†’π‘›=0,onΓ𝑁.(2.1) For simplicity we assume that πœŒπ·β‰‘0 and that the flow is driven by 𝑓. However, in simulation examples we use 𝑓≑0, πœŒπ·β‰’0.

The variational formulation of problem (1.1) and (2.1) reads

Find 𝜌∈𝐢([0,𝑇];𝐿2(Ξ©))∩𝐿2(0,𝑇;𝐻1𝐷(Ξ©))  so thatξ‚΅πœ™πœ•πœŒξ‚Άπœ•π‘‘,πœ“+(π‘˜βˆ‡πœŒ,πœ“)=(𝑓,πœ“),βˆ€πœ“βˆˆπ»1𝐷(Ξ©).(2.2) For any 𝜌0∈𝐿2(Ξ©), π‘“βˆˆπΆ([0,𝑇];𝐿2(Ξ©)), there is a unique solution 𝜌∈𝐢([0,𝑇];𝐿2(Ξ©))∩𝐿2(0,𝑇;𝐻1𝐷(Ξ©)) 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, Μƒπ‘˜1, Μƒπ‘˜2, and π‘βˆˆπΏβˆž(Ξ©) 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 Μƒπ‘˜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 Μƒπ‘˜2>0, we assume the following initial and boundary conditions defined for 𝑖=1,2: ΜƒπœŒπ‘–=ΜƒπœŒπ‘–,0,inΩ×{0},ΜƒπœŒπ‘–=𝜌𝐷,onΓ𝐷,Μƒπ‘˜π‘–βˆ‡ΜƒπœŒπ‘–β‹…β†’π‘›=0,onΓ𝑁.(2.3) As before, we assume that πœŒπ·β‰‘0 for now.

Let ΜƒπœŒπ‘–(β‹…,0)=ΜƒπœŒπ‘–,0∈𝐿2(Ξ©), π‘“βˆˆπΆ([0,𝑇],𝐿2(Ξ©)), 𝑒𝐷≑0. The variational form of (1.2a), (1.2b), and (2.3) is as follows:

Find ΜƒπœŒπ‘–βˆˆπΆ([0,𝑇];𝐿2(Ξ©))∩𝐿2(0,𝑇;𝐻1𝐷(Ξ©)),  𝑖=1,2, such thatξ‚΅ξ‚πœ™1πœ•ΜƒπœŒ1ξ‚Ά+ξ€·Μƒπ‘˜πœ•π‘‘,πœ“1βˆ‡ΜƒπœŒ1ξ€Έ+𝑐,βˆ‡πœ“ΜƒπœŒ1βˆ’ΜƒπœŒ2ξ€Έξ€Έ,πœ“=(𝑓,πœ“),βˆ€πœ“βˆˆπ»1π·ξ‚΅ξ‚πœ™(Ξ©),2πœ•ΜƒπœŒ2ξ‚Ά+ξ€·Μƒπ‘˜πœ•π‘‘,πœ‰2βˆ‡ΜƒπœŒ2ξ€Έ+𝑐,βˆ‡πœ‰ΜƒπœŒ2βˆ’ΜƒπœŒ1ξ€Έξ€Έ,πœ‰=0,βˆ€πœ‰βˆˆπ»1𝐷(Ξ©).(2.4) The well posedness of this system (2.4) follows immediately by showing that the operator β„¬βˆΆπ’ŸβŠ‚π»1𝐷(Ξ©)×𝐻1𝐷(Ξ©)→𝐿2(Ξ©)×𝐿2(Ξ©)βŽ‘βŽ’βŽ’βŽ£Μƒπ‘˜β„¬=βˆ’πœ•1Μƒπ‘˜πœ•+π‘πΌβˆ’π‘πΌβˆ’π‘πΌβˆ’πœ•2⎀βŽ₯βŽ₯βŽ¦πœ•+𝑐𝐼,where𝐼istheidentityoperator,(2.5) is m-accretive. The accretiveness is straightforward by showing that (𝐡𝑒,𝑒)β‰₯0 in the product space. The maximal property follows from that of the operators Μƒπ‘˜βˆ’πœ•π‘–πœ•, 𝐼.

Now, if Μƒπ‘˜2≑0, we do not impose boundary conditions on ΜƒπœŒ2. The variational formulation of (3) and (2.3), reads

Find ΜƒπœŒ1∈𝐢([0,𝑇];𝐿2(Ξ©))∩𝐿2(0,𝑇;𝐻1𝐷(Ξ©)),β€‰β€‰ΜƒπœŒ2∈𝐢([0,𝑇];𝐿2(Ξ©))∩𝐿2(0,𝑇;𝐿2(Ξ©))   such thatξ‚΅ξ‚πœ™1πœ•ΜƒπœŒ1ξ‚Ά+ξ€·Μƒπ‘˜πœ•π‘‘,πœ“1βˆ‡ΜƒπœŒ1ξ€Έ+𝑐,βˆ‡πœ“ΜƒπœŒ1βˆ’ΜƒπœŒ2ξ€Έξ€Έ,πœ“=(𝑓,πœ“),βˆ€πœ“βˆˆπ»1π·ξ‚΅ξ‚πœ™(Ξ©),2πœ•ΜƒπœŒ2ξ‚Ά+ξ€·π‘ξ€·πœ•π‘‘,πœ‰ΜƒπœŒ2βˆ’ΜƒπœŒ1ξ€Έξ€Έ,πœ‰=0,βˆ€πœ‰βˆˆπΏ2(Ξ©).(2.6) 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,𝑇] into subintervals [π‘‘π‘›βˆ’1,𝑑𝑛], 1≀𝑛≀𝑁, such that 0=𝑑0<𝑑1<β‹―<𝑑𝑁=𝑇. We denote Δ𝑑𝑛=π‘‘π‘›βˆ’π‘‘π‘›βˆ’1 and 𝑧(𝑑𝑛)=𝑧𝑛 for any function 𝑧.

For spatial discretization, we adopt standard finite element nomenclature that can be found in textbooks such as [23–25]. At each time step 𝑛, we denote by π’―π‘›β„Ž, β„Ž>0, a family of admissible and shape-regular partitions of Ξ© into a finite number of elements. For any element π‘‡βˆˆπ’―π‘›β„Ž, we let β„Žπ‘‡ be the diameter of 𝑇 and denote β„Ž=maxπ‘‡βˆˆπ’―β„Žβ„Žπ‘‡. Also, β„°π‘›β„Ž is the set of all edges of π’―π‘›β„Ž. For any edge πΈβˆˆβ„°π‘›β„Ž, we let β„ŽπΈ denote the diameter of the edge 𝐸. Finally, we denote by π’«π‘˜(𝑇) the space of polynomials of degree π‘˜ in ℝ𝑑.

Numerical Approximation of Micro-Model
At each time step 𝑛, we define π‘‰π‘›β„Ž=ξ‚†π‘£β„Žξ‚€βˆˆπΆΞ©ξ‚βˆΆβˆ€π‘‡βˆˆπ’―π‘›β„Ž,π‘£β„Žβˆ£π‘‡βˆˆπ’«π‘˜(𝑇),π‘£βˆ£Ξ“π·ξ‚‡=0,(2.7) and seek an approximation πœŒπ‘›β„Žβˆˆπ‘‰π‘›β„Ž such that for all πœ“β„Žβˆˆπ‘‰π‘›β„Žξƒ©πœ™πœŒπ‘›β„Žβˆ’πœŒβ„Žπ‘›βˆ’1Δ𝑑𝑛,πœ“β„Žξƒͺ+ξ€·π‘˜βˆ‡πœŒπ‘›β„Ž,βˆ‡πœ“β„Žξ€Έ=𝑓,πœ“β„Žξ€ΈπœŒ,(2.8a)0β„Ž=β„β„ŽπœŒ0inΞ©.(2.8b)Here β„β„Ž denotes an interpolation or projection operator into π‘‰π‘›β„Ž.
It is standard [20] that there is a unique solution πœŒπ‘›β„Žβˆˆπ‘‰π‘›β„Ž for (2.8a) at each 𝑛, 1≀𝑛≀𝑁. 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 ΜƒπœŒπ‘– can be associated with its own grid π’―π‘›β„Žπ‘–; 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 πœ™, π‘˜ are highly varying coefficients. In particular, consider a porous domain Ξ© composed of three disjoint regions Ω𝑖, 𝑖=1,…,3, see Figure 1, with positive constants (πœ™π‘–,π‘˜π‘–)3𝑖=1 so that πœ™(π‘₯)=πœ™1πœ’1(π‘₯)+πœ™2πœ’2(π‘₯)+πœ™3πœ’3π‘˜(π‘₯),(π‘₯)=π‘˜1πœ’1(π‘₯)+π‘˜2πœ’2(π‘₯)+π‘˜3πœ’3(π‘₯).(3.1) Here πœ’π‘– denotes the characteristic function of Ω𝑖.

An accurate numerical approximation to (1.1) with (3.1) requires that the grid π’―β„Ž is very fine so that it lines up with the interfaces πœ•Ξ©π‘–βˆ©πœ•Ξ©π‘˜. 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 ξ‚πœ™π‘–,Μƒπ‘˜π‘– 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 π’―β„Ž.

3.1. Derivation of Macro-Model

One successful avenue to reduce the complexity of a highly heterogeneous model is to consider a homogenized 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 π‘˜1β‰«π‘˜2 and/or π‘˜1β‰«π‘˜3.

The double-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 Ω𝑖 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 𝑐 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 𝑑=3. In 𝑑=2 only one of these regions, say Ξ©1, can be connected; see Figure 1(a). As we point out below, this implies Μƒπ‘˜2=0, and the macro-model (1.2a) and (1.2b) in 𝑑=2 becomes the Warren-Root relaxation model (3) as in [6].

In this paper we are interested in the double-diffusion models with a constant 𝑐 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, Μƒπ‘˜1, Μƒπ‘˜2, and 𝑐 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 πœ™π‘–, π‘˜π‘– 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 π‘Œ=(0,1)π‘‘βŠ‚β„π‘‘ be divided into three distinct regions π‘Œ1, π‘Œ2, and π‘Œ3, such that πœ•π‘Œ1βˆ©πœ•π‘Œ2=βˆ… and π‘Œ=βˆͺ3𝑖=1π‘Œπ‘– as illustrated by Figure 1(b). Denote by 𝐢𝑖 the characteristic function of π‘Œπ‘–, 𝑖=1,2,3, extended π‘Œ-periodically to all ℝ𝑑. We assume that the domains {π‘₯βˆˆβ„π‘‘βˆΆπΆπ‘–(π‘₯)=1}, 𝑖=1,2,3, have a smooth boundary. Also, we define the πœ€-periodic characteristic functions πΆπœ€π‘–(π‘₯)=𝐢𝑖(π‘₯/πœ€), 𝑖=1,2,3,  π‘₯∈Ω. Thus, for a fixed πœ€, Ξ© is subdivided into three distinct regions Ξ©πœ€π‘–={π‘₯βˆˆΞ©βˆΆπΆπœ€π‘–(π‘₯)=1}, 𝑖=1,2,3. Another subdivision of Ξ© is into 𝑂(πœ€βˆ’π‘‘)-periodic cells.

Macroporosities ξ‚πœ™π‘–
From [11, equation 16, page 206] we have ξ‚πœ™π‘–=βˆ«π‘Œπ‘–πœ™π‘–(𝑦)𝑑𝑦.

Macroconductivities Μƒπ‘˜π‘–
From [11, equations 16, 14, page 205-206] we have Μƒπ‘˜π‘ξ€œ(𝑖,𝑗)=π‘Œπ‘π‘˜π‘ξ‚΅(𝑦)→𝑒𝑖+βˆ‡π‘Šπ‘π‘–ξ‚Άβ‹…ξ‚΅(𝑦)→𝑒𝑗+βˆ‡π‘Šπ‘π‘—ξ‚Ά(𝑦)𝑑𝑦,𝑝=1,2;𝑖,𝑗=1,…𝑑,(3.2) where π‘Šπ‘π‘–, for 𝑝=1,2;𝑖=1,…,𝑑, is the solution to an auxiliary PDEξ‚Έπ‘˜βˆ‡β‹…π‘ξ‚΅β†’π‘’π‘–+βˆ‡π‘Šπ‘π‘–(𝑦)ξ‚Άξ‚Ή=0,inπ‘Œπ‘π‘˜,(3.3a)𝑝→𝑒𝑖+βˆ‡π‘Šπ‘π‘–ξ‚Άβ‹…(𝑦)→𝑛𝑝=0,onΓ𝑝,3π‘Š,(3.3b)𝑝𝑖,π‘˜π‘βˆ‡π‘Šπ‘π‘–β‹…β†’π‘›π‘areperiodic.(3.3c)One can prove [1, 3, 11] that Μƒπ‘˜π‘ is a symmetric matrix.
Now consider 𝑑=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 →𝑒𝑖+βˆ‡π‘Š2𝑖(𝑦)≑0 for 𝑖=1,2. Consequently, Μƒπ‘˜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 Parameter 𝑐
In general, the flow between Ξ©1 and Ξ©2 across Ξ©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 𝑐(ΜƒπœŒ1βˆ’ΜƒπœŒ2). To compute 𝑐, we consider the solution π‘ˆ ofξ€Ίπ‘˜βˆ’βˆ‡β‹…3ξ€»βˆ‡π‘ˆ=0,inπ‘Œ3,(3.4a)π‘ˆ=1,onΞ“1,3,(3.4b)π‘ˆ=0,onΞ“2,3,(3.4c)π‘ˆ,π‘˜3βˆ‡π‘ˆβ‹…β†’π‘›3areperiodiconΞ“3,3.(3.4d)The parameter 𝑐 see [11, equation 19, page 209] is given by ξ€œπ‘=Ξ“1,3π‘˜3βˆ‡π‘¦π‘ˆβ‹…β†’π‘›3𝑑𝑠.(3.5) We note that the general non-quasistatic case is when πœ™3(πœ•π‘ˆ/πœ•π‘‘) 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 β„Žπ‘˜βˆ’11=40 and β„Žπ‘βˆ’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 parameters Μƒπ‘˜π‘–, 𝑐 and provide examples of their calculations for various values of π‘˜3 and choices of geometry of π‘Œ.

In all examples 𝑑=2 hence Μƒπ‘˜2=ξ€Ί0000ξ€». Also, in all examples π‘˜1≑const and the geometry of π‘Œ is axisymmetric, and thus Μƒπ‘˜1 is a scalar constant.

Example 3.1. Here π‘Œ2=(0.3,0.7)2,β€‰β€‰π‘Œ3=(0.2,0.8)2β§΅π‘Œ2,β€‰β€‰πœ™=1πœ’1+10βˆ’4πœ’2,  and π‘˜=1πœ’1+10βˆ’4πœ’2+10βˆ’7πœ’3. We obtain ξ‚πœ™1ξ‚πœ™=0.64,2=1.6Γ—10βˆ’5,Μƒπ‘˜1=0.4519,𝑐=1.8634Γ—10βˆ’6.(3.6)

Example 3.2. Let π‘Œ2=(0.3,0.7)2,β€‰β€‰π‘Œ3=(0.2,0.8)2β§΅π‘Œ2. πœ™=1πœ’1+10βˆ’4πœ’2, and β€‰π‘˜=1πœ’1+10βˆ’1πœ’2+10βˆ’4πœ’3. We obtain ξ‚πœ™1ξ‚πœ™=0.64,2=1.6Γ—10βˆ’5,Μƒπ‘˜1=0.4519,𝑐=1.8634Γ—10βˆ’3.(3.7)

Example 3.3. Here π‘Œ2=(0.3,0.7)2,β€‰β€‰π‘Œ3=(0.2,0.8)2β§΅π‘Œ2,β€‰β€‰πœ™=1πœ’1+10βˆ’4πœ’2,  and π‘˜=1πœ’1+10βˆ’1πœ’2+10βˆ’2πœ’3. We obtain ξ‚πœ™1ξ‚πœ™=0.64,2=1.6Γ—10βˆ’5,Μƒπ‘˜1=0.4519,𝑐=1.8634Γ—10βˆ’1.(3.8)

Example 3.4. Let π‘Œ2=(0.3,0.7)2,β€‰β€‰π‘Œ3=(0.1,0.9)2β§΅π‘Œ2,β€‰β€‰πœ™=1πœ’1+10βˆ’4πœ’2, and π‘˜=1πœ’1+10βˆ’1πœ’2+10βˆ’4πœ’3. We get ξ‚πœ™1ξ‚πœ™=0.36,2=1.6Γ—10βˆ’5,Μƒπ‘˜1=0.2130,𝑐=1.0412Γ—10βˆ’3.(3.9)

Example 3.5. π‘Œ2=(0.25,0.75)2,β€‰β€‰π‘Œ3=(0.2,0.8)2β§΅π‘Œ2,β€‰β€‰πœ™=1πœ’1+10βˆ’4πœ’2, and π‘˜=1πœ’1+10βˆ’1πœ’2+10βˆ’4πœ’3. We get ξ‚πœ™1ξ‚πœ™=0.64,2=2.5Γ—10βˆ’5,Μƒπ‘˜1=0.4519,𝑐=4.3000Γ—10βˆ’3.(3.10)

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 ΜƒπœŒπ‘– to πœ’π‘–πœŒπœ€. However, this notion of convergence does not give information about the rate at which ΜƒπœŒπ‘–βˆ’πœ’π‘–πœŒπœ€ 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 π‘˜3 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 set 𝑓≑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 Γ𝐷={0}Γ—[0,1]βˆͺ{1}Γ—[0,1]. On Γ𝐷 we define 𝜌𝐷(0,𝑦)=1,𝜌𝐷(1,𝑦)=0,π‘¦βˆˆ[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βˆ’π‘₯,1βˆ’π‘₯⟩. Due to the incompatibility between initial and boundary data, we have high gradients of the solution close to π‘₯=1 for small 𝑑.

We fix the final time 𝑇=0.05 and use uniform time stepping with Δ𝑑=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 𝑑𝑛 by πœŒπ‘›πœ€,β„Ž where β„Ž denotes the grid parameter for the micro-model.

For the macro-model, we use ΜƒπœŒπ‘–,0=𝜌0≑1. Also, we use ΜƒπœŒπ·,𝑖=𝜌𝐷 for 𝑖=1,2 when Μƒπ‘˜2β‰’0. If Μƒπ‘˜2≑0, then the boundary condition for ΜƒπœŒ2 is not prescribed.

3.2.2. Qualitative Comparison

In our first comparison we use parameters from Example 3.1. We solve the macro-model (4.2) for {βŸ¨ΜƒπœŒπ‘›1,β„Ž1,ΜƒπœŒπ‘›2,β„Ž2⟩}𝑁𝑛=1 and the micro-model (2.8a) and (2.8b) for πœŒπ‘›πœ€,β„Ž. The plots of πœŒπ‘›β„Ž and ΜƒπœŒπ‘›π‘–,β„Ž are in Figure 2, and a different view is shown in Figure 1(d).

The heterogeneous structure of Ξ©πœ€ is well visible from the behavior of the micro-model solution πœŒπ‘›β„Ž. Large gradients of solution are visible on cell boundaries due to the large difference between π‘˜1, π‘˜2, and π‘˜3.

However, the behavior in the fast region is well approximated by ΜƒπœŒπ‘›1,β„Ž1 which envelopes πœŒπ‘›β„Žπœ’1, while ΜƒπœŒπ‘›2,β„Ž2 envelopes πœŒπ‘›β„Žπœ’2 well. This is very well seen in a side view in Figure 2(d).

Quite interesting is the behavior of 𝜌2|π‘₯=1. Since in our examples Μƒπ‘˜2≑0, no boundary condition at π‘₯=1 is imposed on 𝜌2. Thus, for small 𝑑, 𝜌2(1,𝑑) 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 π‘₯=1. In particular, with constant ξ‚πœ™2, 𝑐, we have ΜƒπœŒ2(π‘₯,𝑦,𝑑)=ΜƒπœŒ2,0(π‘₯,𝑦)π‘’βˆ’(𝑐/ξ‚πœ™2)𝑑+ξ€œπ‘‘0ΜƒπœŒ1(π‘₯,𝑦,𝑠)π‘’βˆ’(𝑐/ξ‚πœ™2)(π‘‘βˆ’π‘ )𝑑𝑠.(3.11) Thus, for small 𝑑 and small ξ‚πœ™π‘/2, ΜƒπœŒ2(1,𝑦,𝑑) is away from 0, but, as time increases, its magnitude decreases proportionally to π‘’βˆ’(𝑐/ξ‚πœ™2)𝑑.

Next, we use parameters from Example 3.2 and plot them in Figure 3. In this example the conductivity π‘˜1 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.

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/(2𝑙+1),  𝑙β‰₯1. We use uniform mesh with β„Ž=β„Ž1=β„Ž2=πœ€/20 for now. See Section 3.2.4 for β„Žπ‘–β‰«β„Ž 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 πœ’π‘– are involved. While the results in [11] use two-scale convergence with πœ€, that is, rely on special periodic test functions, we compare πœŒπ‘›πœ€,β„Ž and ΜƒπœŒπ‘›π‘–,β„Ž directly. We see that it is easy to do so only in the connected region Ξ©1.

We define a quantity π‘’π‘›πœ€=(πœŒπ‘›πœ€,β„Žβˆ’ΜƒπœŒπ‘›1,β„Ž1)πœ’1 to be used in comparison. We also define, for a fixed π‘‘βˆˆ(0,𝑇] and some function 𝑧, ‖‖𝑧(π‘₯,𝑦;𝑑)βˆ—,𝐿1ξ€œβˆΆ=10||||‖𝑧(π‘₯,0.5;𝑑)𝑑π‘₯,𝑧(π‘₯,𝑦;𝑑)β€–βˆ—,𝐿2ξ‚»ξ€œβˆΆ=10||||𝑧(π‘₯,0.5;𝑑)2𝑑π‘₯1/2,(3.12) which will be applied to π‘’π‘›πœ€. Clearly β€–β‹…β€–βˆ—,π‘Š is not a norm but is a useful quantity of interest.

In Tables 1, 2, 3, and 4 we present β€–π‘’π‘›πœ€β€–βˆ—,π‘Š for different πœ€ and π‘Š=𝐿2,𝐿1. For each case we consider convergence of β€–π‘’π‘›πœ€β€–βˆ—,π‘Š with πœ€β†“0.

From homogenization theory it is not clear what order of convergence should be expected. We observe that β€–π‘’π‘›πœ€β€–βˆ—,π‘Š decreases linearly with πœ€ for both βˆ—=𝐿1,𝐿2 -like quantities in all cases.

Next we discuss the dependence of the results on conductivity and other data. Consider Example 3.2 as a reference example with its corresponding Table 1. Compare it now with that from Table 2 where π‘˜3 is larger than that in the reference case to see how this influences the errors. As expected, and as shown in Example 3.3, 𝑐 changes by the same factor, and this means that the coupling between ΜƒπœŒ1 and ΜƒπœŒ2 is stronger which corresponds to faster flow across π‘Œ3 in the micro-model. However, the influence on the convergence of the micro-macro difference, from Table 2 is rather weak.

Next consider Example 3.4 which uses a different geometry than that in the reference case, with thicker region π‘Œ3. This produces a slightly smaller 𝑐 since the gradient of π‘ˆ is smaller. Also, ξ‚πœ™1 and Μƒπ‘˜1 are predictably different. The error quantity is influenced insignificantly.

The opposite effect is seen in Table 4 when π‘Œ3 is selected as thinner than in the reference Example 3.2.

3.2.4. Coarse Macrogrid and Parameter Grid

In the examples presented in Section 3.2.3 we used β„Ž=β„Ž1=β„Ž2. 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 β„Žβ‰ͺβ„Ž1 in the macro-model.

Example 3.6. We use data from Example 3.2 with πœ€=0.2. Consider a fixed microgrid with β„Žβˆ’1=100 and a few cases of macrogrid with β„Ž1βˆ’1=β„Ž2βˆ’1=10,50,100. Now we compare the micro- and macrosolutions πœŒπ‘›πœ€,β„Ž and 𝜌1,β„Ž1, see Table 5. We also compare the run time of the macro-model to that of the reference micro-model which is 120 seconds.

We see that, as predicted above, the macro-model run time at β„Žπ‘–=β„Ž 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 β„Žπ‘–β‰«β„Ž, 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 β„Žπ‘– is by a factor of 10 coarser than β„Ž 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 with β„Ž1β‰ β„Ž2 are implemented, see Section 4.

Influence of Parameter Grid
Last we address the effect of the computed parameters 𝑐, Μƒπ‘˜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 β„Žβˆ’1=60. Use β„Ž1βˆ’1=β„Ž2βˆ’1=60. Compute numerical parameters π‘β„Žπ‘, Μƒπ‘˜1,β„Žπ‘˜1 by approximating solutions to (3.4a), (3.4b), (3.4c), and (3.4d) and (3.3a), (3.3b), and (3.3c), respectively.
Example 3.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.

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 that Μƒπ‘˜1β‰«Μƒπ‘˜2, that is, that the first component varies in space faster than the second. As a special case, this includes the case Μƒπ‘˜2=0. We consider two triangulations π’―β„Žπ‘–, 𝑖=1,2, which will be used for ΜƒπœŒβ„Žπ‘–, respectively. In principle, these can be chosen independently. In fact, we allow the triangulations to vary in time and thus consider π’―π‘›β„Žπ‘– and the associated spaces π‘‰π‘›β„Žπ‘– as in (2.7). To avoid the loss of accuracy due to excessive interpolation and intergrid projections and because Μƒπ‘˜1β‰«Μƒπ‘˜2, we assume that π’―π‘›β„Ž1 is a refinement of the partition π’―π‘›β„Ž2.

We need intergrid operators to handle two components that live on separate grids. Let Ξ βˆΆπ‘‰π‘›β„Ž2β†’π‘‰π‘›β„Ž1 be the interpolation operator and Ξ β€²βˆΆπ‘‰π‘›β„Ž1β†’π‘‰π‘›β„Ž2 the 𝐿2 projection defined by ξ‚€Ξ ξ…žπœ™β„Ž1,πœ“β„Ž2ξ‚ξ€·πœ™βˆΆ=β„Ž1,Ξ πœ“β„Ž2ξ€Έ,βˆ€πœ“β„Ž2βˆˆπ‘‰β„Ž2.(4.1)

Now we define the discrete solutions. At each time step 1≀𝑛≀𝑁, we seek ΜƒπœŒπ‘›π‘–,β„Žπ‘–βˆˆπ‘‰π‘›β„Žπ‘–, 𝑖=1,2, which satisfy the discrete problem, that is, (2.4) restricted to the finite dimensional subspaces π‘‰π‘›β„Žπ‘– so that, for all πœ“β„Ž1βˆˆπ‘‰π‘›β„Ž1, for all πœ‰β„Ž2βˆˆπ‘‰π‘›β„Ž2ξƒ©ξ‚πœ™1ΜƒπœŒπ‘›1,β„Ž1βˆ’ΜƒπœŒπ‘›βˆ’11,β„Ž1Δ𝑑𝑛,πœ“β„Ž1ξƒͺ+ξ‚€Μƒπ‘˜1βˆ‡ΜƒπœŒπ‘›1,β„Ž1,βˆ‡πœ“β„Ž1+ξ‚€π‘ξ‚€ΜƒπœŒπ‘›1,β„Ž1βˆ’Ξ ΜƒπœŒπ‘›2,β„Ž2,πœ“β„Ž1=𝑓,πœ“β„Ž1ξ€Έ,ξƒ©ξ‚πœ™2ΜƒπœŒπ‘›2,β„Ž2βˆ’ΜƒπœŒπ‘›βˆ’12,β„Ž2Δ𝑑𝑛,πœ‰β„Ž2ξƒͺ+ξ‚€Μƒπ‘˜2βˆ‡ΜƒπœŒπ‘›2,β„Ž2,βˆ‡πœ‰β„Ž2+ξ‚€π‘ξ‚€ΜƒπœŒπ‘›2,β„Ž2βˆ’Ξ ξ…žΜƒπœŒπ‘›1,β„Ž1,πœ‰β„Ž2=0,ΜƒπœŒ01,β„Ž1,ΜƒπœŒ02,β„Ž2ξ‚­=ξ«β„β„Ž1ΜƒπœŒ1,0,β„β„Ž2ΜƒπœŒ2,0inΞ©.(4.2) If Μƒπ‘˜2=0, the system (4.2) is modified appropriately so that ΜƒπœŒ2(β‹…,𝑑)∈𝐿2(Ξ©) instead of 𝐻1𝐷(Ξ©), and we do not impose boundary conditions on this component. In fact, instead of (2.7) we define π‘‰π‘›β„Ž2=ξ‚†π‘£β„Žξ‚€βˆˆπΆΞ©ξ‚βˆΆβˆ€π‘‡βˆˆπ’―π‘›β„Ž,π‘£β„Žβˆ£π‘‡βˆˆπ’«π‘˜ξ‚‡(𝑇).(4.3) (It holds that π‘£β„ŽβˆˆπΆ(Ξ©) 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 follows 1≀𝑛≀𝑁 and 𝑖=1,2. We denote 𝑒𝑖,β„Žπ‘–(𝑑𝑛)=ΜƒπœŒπ‘›π‘–,β„Žπ‘–βˆ’ΜƒπœŒπ‘›π‘– and define the energy norm for the product space ||||β€–βŸ¨π‘§,π‘€βŸ©β€–2ξ€·π‘‘π‘›ξ€Έξ‚πœ™βˆΆ=1‖‖𝑧𝑑𝑛‖‖2+ξ‚πœ™2‖‖𝑀𝑑𝑛‖‖2+π‘›ξ“π‘š=1Ξ”π‘‘π‘šξ‚€Μƒπ‘˜1||π‘§ξ€·π‘‘π‘šξ€Έ||21+Μƒπ‘˜2||π‘€ξ€·π‘‘π‘šξ€Έ||21.(4.4)

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,β„Ž1,ΜƒπœŒπ‘›2,β„Ž2⟩}𝑁𝑛=1 be the smooth solutions of (2.4) and (4.2), respectively. Then one has ||‖‖𝑒1,β„Ž1,𝑒2,β„Ž2‖‖||ξ€·π‘‘π‘›ξ€Έβ‰€Μƒπ‘˜11/2𝐢1β„Ž1+Μƒπ‘˜21/2𝐢2β„Ž2+𝐢3maxπ‘š=1βˆΆπ‘›Ξ”π‘‘π‘š+higherorderterms(h.o.t.),(4.5) where 𝐢1, 𝐢2, and 𝐢3 are independent of β„Ž.
Furthermore, if Μƒπ‘˜2=0, then ||‖‖𝑒1,β„Ž1,𝑒2,β„Ž2‖‖||ξ€·π‘‘π‘›ξ€Έβ‰€Μƒπ‘˜11/2𝐢1β„Ž1+𝐢2β„Ž22+𝐢3maxπ‘š=1βˆΆπ‘›Ξ”π‘‘π‘š+higherorderterms(β„Ž.π‘œ.𝑑.).(4.6)

Now we notice that (4.5) suggests the following choice of grid parameters β„Ž1, β„Ž2. If Μƒπ‘˜1β‰«Μƒπ‘˜2, then to balance the components of the error one can use β„Ž1β‰ͺβ„Ž2. Similarly, if Μƒπ‘˜2=0, then (4.6) implies that β„Ž2 can be chosen to be of the order of β„Ž11/2. The use of coarse grid β„Ž2β‰«β„Ž1 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 πœ‚π‘›βˆ‘=(π‘›π‘š=1[Ξ”π‘‘π‘šπ‘†2π‘š+𝑇2π‘š])1/2 for (4.2) is composed of the temporal part 𝑇𝑛 which adapts the time stepping and the spatial part 𝑆𝑛 which guides the spatial grid adaptivity. We define 𝑇2π‘›βˆΆ=Δ𝑑𝑛3ξ‚΅β€–β€–Μƒπ‘˜11/2ξ‚€ΜƒπœŒπ‘›1,β„Ž1βˆ’ΜƒπœŒπ‘›βˆ’11,β„Ž1‖‖2+β€–β€–Μƒπ‘˜21/2ξ‚€ΜƒπœŒπ‘›2,β„Ž2βˆ’ΜƒπœŒπ‘›βˆ’12,β„Ž2‖‖2ξ‚Ά+Δ𝑑𝑛3‖‖𝑐1/2ξ‚€ΜƒπœŒπ‘›1,β„Ž1βˆ’ΜƒπœŒπ‘›βˆ’11,β„Ž1ξ‚βˆ’π‘1/2ξ‚€ΜƒπœŒπ‘›2,β„Ž2βˆ’ΜƒπœŒπ‘›βˆ’12,β„Ž2‖‖2,(4.7) and 𝑆2𝑛=𝑆1,𝑛+𝑆2,𝑛, where 𝑆𝑖,𝑛=βˆ‘π‘‡π‘–βˆˆπ’―β„Žπ‘–π‘†π‘›π‘–,𝑇𝑖, and where, as usual (see [23]), one defines 𝑆𝑛𝑖,π‘‡π‘–βˆΆ=πœƒ2𝑖,𝑇𝑖‖‖𝑅𝑛𝑖,𝑇𝑖‖‖2+12ξ“πΈβˆˆβ„°π‘‡π‘–π›Ύ2𝑖,𝐸‖‖𝑅𝑛𝑖,𝐸‖‖2.(4.8) The element and edge residuals in 𝑆𝑛𝑖,𝑇𝑖,𝑖=1,2, are given by 𝑅𝑛𝑖,𝑇𝑖=π‘“π‘–βˆ’ξ‚πœ™π‘–ΜƒπœŒπ‘›π‘–,β„Žπ‘–βˆ’ΜƒπœŒπ‘›βˆ’1𝑖,β„Žπ‘–Ξ”π‘‘π‘›ξ‚€Μƒπ‘˜+βˆ‡β‹…π‘–βˆ‡ΜƒπœŒπ‘›π‘–,β„Žπ‘–ξ‚+(βˆ’1)π‘–π‘ξ‚€ΜƒπœŒπ‘›1,β„Ž1βˆ’ΜƒπœŒπ‘›2,β„Ž2,𝑅(4.9)𝑛𝑖,𝐸=⎧βŽͺ⎨βŽͺβŽ©ξ‚ƒΜƒπ‘˜π‘–πœ•πœˆΜƒπœŒπ‘›π‘–,β„Žπ‘–ξ‚„πΈβˆ’Μƒπ‘˜ifπΈβŠ„πœ•Ξ©,π‘–πœ•πœˆΜƒπœŒπ‘›π‘–,β„Žπ‘–ifπΈβŠ‚Ξ“π‘,0,otherwise.(4.10) Here [𝑧]𝐸 denotes the jump of the flux 𝑧 across the edge 𝐸 of an element. Usually 𝑓1≑𝑓, 𝑓2≑0. Also, if Μƒπ‘˜2=0, then 𝑅𝑛2,𝐸≑0.

The scaling constants πœƒπ‘–,𝑇,𝛾𝑖,𝐸 take into account the contribution of the exchange term 𝑐(ΜƒπœŒ1βˆ’ΜƒπœŒ2) and are defined so that the estimators remain robust when the coefficients of the problem change substantially: πœƒ1,π‘‡ξ‚†β„Ž=min𝑇1Μƒπ‘˜1βˆ’1/2𝑐,maxβˆ’1/2,β„Žπ‘‡1Μƒπ‘˜2βˆ’1/2,𝑇1βˆˆπ’―π‘›β„Ž1βˆͺβ„°π‘›β„Ž1πœƒ,(4.11)2,π‘‡ξ‚†β„Ž=min𝑇2Μƒπ‘˜2βˆ’1/2𝑐,maxβˆ’1/2,β„Žπ‘‡2Μƒπ‘˜1βˆ’1/2,𝑇2βˆˆπ’―π‘›β„Ž2βˆͺβ„°π‘›β„Ž2𝛾,(4.12)𝑖,𝐸=2β„ŽπΈβˆ’1/2πœƒπ‘–,𝐸,𝑖=1,2.(4.13) This a posteriori estimator works well with two-grid discretizations as shown in the examples to follow. It is well formulated also for Μƒπ‘˜2≑0 when ΜƒπœŒ2∈𝐿2(Ξ©). We have the following result on reliability of the estimator proven in the appendix.

Theorem 4.2. Let βŸ¨ΜƒπœŒ1,ΜƒπœŒ2⟩, {βŸ¨ΜƒπœŒπ‘›1,β„Ž1,ΜƒπœŒπ‘›2,β„Ž2⟩}𝑁𝑛=1 be the solution of (1.2a) and (1.2b), (2.3) and (4.2), respectively. Then ||‖‖𝑒1,β„Ž1,𝑒2,β„Ž2‖‖||ξ€·π‘‘π‘›ξ€Έβ‰€πœ‚π‘›+h.o.t.;βˆ€π‘›β‰₯1(4.14)

4.3. Adaptive Two-Grid Discretization and Implementation

Consider a fixed time step 𝑛=π‘›π‘œ for which some triangulation is chosen, the solution is found, and the spatial error indicators 𝑆𝑖,𝑛 are computed. To apply the two-grid spatial adaptivity we use the following refinement algorithm.

Adaptive Two-Grid Algorithm 𝓐
(1)Select the triangulation(s) to be refined: If 𝑆1,π‘›π‘œ>3𝑆2,π‘›π‘œ, refine only π’―π‘›π‘œβ„Ž1. If 𝑆2,π‘›π‘œ>3𝑆1,π‘›π‘œ, refine only π’―π‘›π‘œβ„Ž2. Otherwise, refine both. (2)In π’―β„Žπ‘–, 𝑖=1,2 selected in Step 1, refine any π‘‡βˆˆπ’―β„Ž,𝑖 for which π‘†π‘›π‘œπ‘–,𝑇β‰₯0.5maxπ‘‡βˆˆπ’―β„Žπ‘–π‘†π‘›π‘œπ‘–,𝑇. (3)Enforce the requirement that π’―π‘›β„Ž1 is a refinement of π’―π‘›β„Ž2 by adding extra elements as needed.
The steps in π’œ should be repeated a certain number π‘π’œ 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 Example 3.2, with 𝑇=0.02 and uniform time stepping with Δ𝑑=10βˆ’4. Other data is as in Section 3.2.1.
Consider a fixed 𝑛0=200 with a uniform triangulation π’―π‘›π‘œβ„Ž1=π’―π‘›π‘œβ„Ž2 with β„Ž1=β„Ž2=0.1. Then apply 𝑁𝐴=6 times the steps in π’œ, see the solution and the refined meshes in Figure 4. The process works as expected: since at the time 𝑑𝑛0 the solution has a high gradient near π‘₯=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 Μƒπ‘˜2=0. This can be compared to the a-priori estimate in Theorem 4.1 from which we know the convergence is 𝑂(β„Ž1+β„Ž22+Δ𝑑𝑛).
In Table 7 we provide the details on π’―π‘›π‘œβ„Žπ‘– and compare the effectiveness of the local two-grid refinement by algorithm π’œ with uniform grid refinement. With adaptive refinement we get π‘†π‘›π‘œ=0.35834 with 435+121=556 unknowns, while the uniform refinement needs 441+441=882 unknowns to achieve comparable π‘†π‘›π‘œ=0.39221.

Last, we describe the process with which π’―π‘›β„Žπ‘–,  𝑖=1,2, may vary between time steps. For a fixed 𝑛 denote by π’―π‘›β„Žπ‘– the final triangulation obtained by the adaptive algorithm π’œ. For the new time step 𝑛+1, we use π’―π‘›β„Žπ‘– as the initial triangulation. If the algorithm π’œ suggests that π’―β„Žπ‘›+1𝑖 is modified, we need to project ΜƒπœŒπ‘›π‘–,β„Žπ‘– to the new grid so it can be used as initial data for the new step, but ΜƒπœŒπ‘›π‘–,β„Žπ‘– differs from their interpolation only in the elements where the triangulation is coarsened. Another way to deal with different triangulations π’―β„Žπ‘›+1𝑖, π’―π‘›β„Žπ‘– is to compute the a posteriori error estimator in a common refinement of the two triangulations π’―β„Žπ‘›+1𝑖, π’―π‘›β„Žπ‘– [29, 42], but this will not be pursued further.

The decision to use two grids implies that we have to interpolate between the finite element spaces π‘‰π‘›β„Ž1 and π‘‰π‘›β„Ž2, see (4.1). Suppose that {πœπ‘–,𝑗}𝑗,  𝑖=1,2, is a basis for π‘‰π‘›β„Žπ‘–. Then the operator Ξ  in (4.1) is represented by the matrix π•€β„Ž1,β„Ž2 whose components are given by π•€β„Ž1,β„Ž2(𝑗,π‘˜)=(𝜁1,𝑗,𝜁2,π‘˜).
Consider now and compare the computational cost associated with one-grid and two-grid approaches. To solve (4.2) with one-grid π‘‰π‘›β„Ž1 we need to assemble the stiffness matrix 𝑆1,β„Ž1Μƒπ‘˜(𝑗,π‘˜)=(1βˆ‡πœ1,𝑗,βˆ‡πœ1,π‘˜), the mass matrix 𝑀1,β„Ž1(𝑗,π‘˜)=(𝜁1,𝑗,𝜁1,π‘˜), and the block matrix 𝐡one-grid=βŽ‘βŽ’βŽ’βŽ£βˆ’Ξ”π‘‘π‘†1,β„Ž1+ξ‚€ξ‚πœ™1𝑀+𝑐Δ𝑑1,β„Ž1βˆ’π‘Ξ”π‘‘π‘€1,β„Ž1βˆ’π‘Ξ”π‘‘π‘€1,β„Ž1ξ‚€ξ‚πœ™2𝑀+𝑐Δ𝑑1,β„Ž1⎀βŽ₯βŽ₯⎦.(4.15) At each time step we solve the linear system with 𝐡oneβˆ’grid; this is referred to as π•Šone-grid.
To solve (4.2) with two grids and two spaces π‘‰π‘›β„Ž1,π‘‰π‘›β„Ž2, we need to assemble the stiffness matrix 𝑆1,β„Ž1Μƒπ‘˜(𝑗,π‘˜)=(1βˆ‡πœ1,𝑗,βˆ‡πœ1,π‘˜), the mass matrices 𝑀1,β„Ž1(𝑗,π‘˜)=(𝜁1,𝑗,𝜁1,π‘˜), 𝑀2,β„Ž2(𝑗,π‘˜)=(𝜁2,𝑗,𝜁2,π‘˜), the interpolation matrix π•€β„Ž1,β„Ž2(𝑗,π‘˜)=(𝜁1,𝑗,𝜁2,π‘˜), and the block matrix 𝐡two-grid=βŽ‘βŽ’βŽ’βŽ£βˆ’Ξ”π‘‘π‘†1,β„Ž1+ξ‚€ξ‚πœ™1𝑀+𝑐Δ𝑑1,β„Ž1βˆ’π‘Ξ”π‘‘π•€π‘‡β„Ž1,β„Ž2βˆ’π‘Ξ”π‘‘π•€β„Ž1,β„Ž2ξ‚€ξ‚πœ™2𝑀+𝑐Δ𝑑2,β„Ž2⎀βŽ₯βŽ₯⎦.(4.16) To finish we solve the linear system which we refer to as π•Štwo-grid.

Suppose now we solve (4.2) and wish to maintain 𝑆𝑛≀𝛿 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 β„Ž1=β„Ž2=1/80. Using the two-grid approach we only need β„Ž1=1/80,β„Ž2=1/20.

The computational time needed to assemble the matrices and solve the systems in our MATLAB implementation is displayed in Table 8. 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.

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.


A. Proof of Theorem 4.1

Proof. The proof follows the techniques presented in [25, 26] for a scalar PDE. We present an outline of the proof for the double-diffusion system. By (4.1), if π’―π‘›β„Ž1 is a refinement of π’―π‘›β„Ž2, we can rewrite (4.2) as ξƒ©ξ‚πœ™1ΜƒπœŒπ‘›1,β„Ž1βˆ’ΜƒπœŒπ‘›βˆ’11,β„Ž1Δ𝑑𝑛,πœ“β„Ž1ξƒͺ+ξ‚€Μƒπ‘˜1βˆ‡ΜƒπœŒπ‘›1,β„Ž1,βˆ‡πœ“β„Ž1+ξ‚€π‘ξ‚€ΜƒπœŒπ‘›1,β„Ž1βˆ’ΜƒπœŒπ‘›2,β„Ž2,πœ“β„Ž1=𝑓,πœ“β„Ž1ξ€Έ,ξƒ©ξ‚πœ™2ΜƒπœŒπ‘›2,β„Ž2βˆ’ΜƒπœŒπ‘›βˆ’12,β„Ž2Δ𝑑𝑛,πœ‰β„Ž2ξƒͺ+ξ‚€Μƒπ‘˜2βˆ‡ΜƒπœŒπ‘›2,β„Ž2,βˆ‡πœ‰β„Ž2+ξ‚€π‘ξ‚€ΜƒπœŒπ‘›2,β„Ž2βˆ’ΜƒπœŒπ‘›1,β„Ž1,πœ‰β„Ž2=0.(A.1) Next, for Μƒπ‘˜2β‰’0, as in [26] we define the elliptic projection βŸ¨π‘ƒΜƒπœŒ1,π‘ƒΜƒπœŒ2⟩ of βŸ¨ΜƒπœŒ1,ΜƒπœŒ2⟩ into π‘‰β„Ž1Γ—π‘‰β„Ž2 via ξ€·βˆ‡ξ€·π‘ƒΜƒπœŒπ‘–βˆ’ΜƒπœŒπ‘–ξ€Έξ€Έ,βˆ‡πœ’=0,βˆ€πœ’βˆˆπ‘‰β„Žπ‘–,𝑖=1,2.(A.2) Applying (A.2) for 𝑖=1,2 into the weak formulation (2.4), with arbitrary πœ“,πœ‰βˆˆπ»1𝐷(Ξ©) we arrive at ξ‚΅ξ‚πœ™1πœ•π‘ƒΜƒπœŒ1ξ‚Ά+ξ€·Μƒπ‘˜πœ•π‘‘,πœ“1βˆ‡π‘ƒΜƒπœŒ1ξ€Έ+𝑐,βˆ‡πœ“ΜƒπœŒ1βˆ’ΜƒπœŒ2ξ€Έξ€Έξƒ©ξ‚πœ™,πœ“=(𝑓,πœ“)+1πœ•ξ€·π‘ƒΜƒπœŒ1βˆ’ΜƒπœŒ1ξ€Έξƒͺ,ξ‚΅ξ‚πœ™πœ•π‘‘,πœ“2πœ•π‘ƒΜƒπœŒ2ξ‚Ά+ξ€·Μƒπ‘˜πœ•π‘‘,πœ‰2βˆ‡π‘ƒΜƒπœŒ2ξ€Έ+𝑐,βˆ‡πœ‰ΜƒπœŒ2βˆ’ΜƒπœŒ1ξ€Έξ€Έ=ξƒ©ξ‚πœ™,πœ‰2πœ•ξ€·π‘ƒΜƒπœŒ2βˆ’ΜƒπœŒ2ξ€Έξƒͺ.πœ•π‘‘,πœ‰(A.3)
Next we define πœ•π‘›π‘€=(𝑀(𝑑𝑛)βˆ’π‘€(π‘‘π‘›βˆ’1))/Δ𝑑𝑛, Δ𝑛(𝑀)=πœ•π‘›π‘€βˆ’(πœ•π‘€/πœ•π‘‘)(𝑑𝑛), subtract (A.3) from (A.1), add the resulting equations, and use as test functions πœ“=ΜƒπœŒπ‘›1,β„Ž1βˆ’π‘ƒΜƒπœŒπ‘›1 and πœ‰=ΜƒπœŒπ‘›2,β„Ž2βˆ’π‘ƒΜƒπœŒπ‘›2 to get ξ‚πœ™1β€–β€–ΜƒπœŒπ‘›1,β„Ž1βˆ’π‘ƒΜƒπœŒπ‘›1β€–β€–2+ξ‚πœ™2β€–β€–ΜƒπœŒπ‘›2,β„Ž2βˆ’π‘ƒΜƒπœŒπ‘›2β€–β€–2+Μƒπ‘˜1Δ𝑑𝑛|||ΜƒπœŒπ‘›1,β„Ž1βˆ’π‘ƒΜƒπœŒπ‘›1|||21+Μƒπ‘˜2Δ𝑑𝑛|||ΜƒπœŒπ‘›2,β„Ž2βˆ’π‘ƒΜƒπœŒπ‘›2|||21β‰€ξ‚πœ™1ξ‚€ΜƒπœŒπ‘›βˆ’11,β„Ž1βˆ’π‘ƒΜƒπœŒ1π‘›βˆ’1,ΜƒπœŒπ‘›1,β„Ž1βˆ’π‘ƒΜƒπœŒπ‘›1+ξ‚πœ™2ξ‚€ΜƒπœŒπ‘›βˆ’12,β„Ž2βˆ’π‘ƒΜƒπœŒ2π‘›βˆ’1,ΜƒπœŒπ‘›2,β„Ž2βˆ’π‘ƒΜƒπœŒπ‘›2+Ξ”π‘‘π‘›ξƒ©ξ‚πœ™1πœ•ξ€·π‘ƒΜƒπœŒ1βˆ’ΜƒπœŒ1ξ€Έπœ•π‘‘,ΜƒπœŒπ‘›1,β„Ž1βˆ’π‘ƒΜƒπœŒπ‘›1ξƒͺ+Ξ”π‘‘π‘›ξƒ©ξ‚πœ™2πœ•ξ€·π‘ƒΜƒπœŒ2βˆ’ΜƒπœŒ2ξ€Έπœ•π‘‘,ΜƒπœŒπ‘›2,β„Ž2βˆ’π‘ƒΜƒπœŒπ‘›2ξƒͺβˆ’Ξ”π‘‘π‘›ξ‚€ξ‚πœ™1Ξ”π‘›ξ€·π‘ƒΜƒπœŒ1ξ€Έ,ΜƒπœŒπ‘›1,β„Ž1βˆ’π‘ƒΜƒπœŒπ‘›1ξ‚βˆ’Ξ”π‘‘π‘›ξ‚€ξ‚πœ™2Ξ”π‘›ξ€·π‘ƒΜƒπœŒ2ξ€Έ,ΜƒπœŒπ‘›2,β„Ž2βˆ’π‘ƒΜƒπœŒπ‘›2.(A.4) With the help of standard inequalities, summing the equations from 𝑛=1…,𝑁, and applying the discrete Gronwall's lemma with ξ‚πœ™π›½=2max{1,ξ‚πœ™2}, we get |||β€–β€–ξ‚¬ΜƒπœŒπ‘›1,β„Ž1βˆ’π‘ƒΜƒπœŒπ‘›1,ΜƒπœŒπ‘›2,β„Ž2βˆ’π‘ƒΜƒπœŒπ‘›2ξ‚­β€–β€–|||2β‰€π‘’π›½π‘‘π‘ξƒ―ξ‚πœ™1β€–β€–ΜƒπœŒ01,β„Ž1βˆ’π‘ƒΜƒπœŒ01β€–β€–2+ξ‚πœ™2β€–β€–ΜƒπœŒ02,β„Ž2βˆ’π‘ƒΜƒπœŒ02β€–β€–2+𝑁𝑛=1Ξ”π‘‘π‘›ξ‚πœ™1ξƒ©β€–β€–β€–β€–πœ•ξ€·π‘ƒΜƒπœŒ1βˆ’ΜƒπœŒ1ξ€Έβ€–β€–β€–β€–πœ•π‘‘2+β€–β€–Ξ”π‘›ξ€·π‘ƒΜƒπœŒ1ξ€Έβ€–β€–2ξƒͺ+𝑁𝑛=1Ξ”π‘‘π‘›ξ‚πœ™2ξƒ©β€–β€–β€–πœ•(π‘ƒΜƒπœŒ2βˆ’ΜƒπœŒ2)β€–β€–β€–πœ•π‘‘2+‖‖Δ𝑛(π‘ƒΜƒπœŒ2)β€–β€–2.ξƒͺξƒ°(A.5) Next step is to estimate the quantities β€–πœ•(π‘ƒΜƒπœŒπ‘–βˆ’ΜƒπœŒπ‘–)/πœ•π‘‘β€–,‖Δ𝑛(π‘ƒΜƒπœŒπ‘–)β€–, 𝑖=1,2. This is done using standard properties of elliptic projections as in [25, Lemmas  13.2 and 13.4, pages 233 and 241] assuming ΜƒπœŒπ‘–,πœ•ΜƒπœŒπ‘–/πœ•π‘‘ are smooth enough. Combining these we arrive at the desired inequality (4.5).
For the case Μƒπ‘˜2≑0, we define π‘ƒβˆ—ΜƒπœŒ2 as the 𝐿2-projection of ΜƒπœŒ2 over π‘‰π‘›β„Ž2 via ξ€·π‘ƒβˆ—ΜƒπœŒ2βˆ’ΜƒπœŒ2ξ€Έ,πœ’=0βˆ€πœ’βˆˆπ‘‰π‘›β„Ž2.(A.6) The proof follows along the same lines as for Μƒπ‘˜2β‰’0 except that we do not (and cannot) use any terms with βˆ‡ΜƒπœŒ2.

B. Proof of Theorem 4.2

Proof. The proof is a combination of techniques in [13–15] which we extend to a coupled system and propose an estimator robust with respect to the five parameters of the problem.
First we define the semidiscrete problem: find {βŸ¨ΜƒπœŒπ‘›1,ΜƒπœŒπ‘›2⟩}𝑁𝑛=1∈(𝐻1𝐷(Ξ©))2 so that for all βŸ¨πœ“,πœ‰βŸ©βˆˆ(𝐻1𝐷(Ξ©))2ξ‚€ξ‚πœ™1ΜƒπœŒπ‘›1+ξ€·Μƒπ‘˜,πœ“1βˆ‡ΜƒπœŒπ‘›1ξ€Έ+𝑐,βˆ‡πœ“ΜƒπœŒπ‘›1βˆ’ΜƒπœŒπ‘›2ξ€Έξ€Έξ‚€ξ‚πœ™,πœ“=(𝑓,πœ“)+1ΜƒπœŒ1π‘›βˆ’1,ξ‚€ξ‚πœ™,πœ“2ΜƒπœŒπ‘›2+ξ€·Μƒπ‘˜,πœ‰2βˆ‡ΜƒπœŒπ‘›2ξ€Έ+𝑐,βˆ‡πœ‰ΜƒπœŒπ‘›2βˆ’ΜƒπœŒπ‘›1ξ€Έξ€Έ=ξ‚€ξ‚πœ™,πœ‰2ΜƒπœŒ2π‘›βˆ’1.,πœ‰(B.1) We estimate the error between the solution of (2.4) and (B.1) and this is how the time estimator 𝑇𝑛 (4.7) arises; see [15] for details.
Now the semidiscrete system (B.1) is a stationary coupled reaction-diffusion model of a type we considered in [13] and for which we proposed a spatial a posteriori error estimator accounting for the coupling terms in the system. That estimator is robust, that is, the efficiency ratio remains essentially constant when the coefficients change by orders of magnitude; this is relevant for our problem since we may have small Μƒπ‘˜2. The robustness in [13] is achieved by an appropriate scaling of constants in the estimator; the scaling in (4.11)–(4.13) serves the same purpose.
Recall the standard set-up first. Fix 𝑖=1,2, and for any π‘‡βˆˆπ’―π‘›β„Žπ‘– denote by ξ‚πœ”π‘‡ the set of all elements in π’―π‘›β„Žπ‘– that share at least one vertex with 𝑇. For the quasi-interpolator π‘„β„Ž defined in [14, Lemma  3.1, page 482], we have that for any π‘£βˆˆπ»π‘˜(𝑀𝑇), 0β‰€π‘˜β‰€1, β€–β€–βˆ‡π‘™ξ€·π‘£βˆ’π‘„β„Žπ‘£ξ€Έβ€–β€–π‘‡β‰€πΆβ„Žπ‘‡π‘˜βˆ’π‘™β€–β€–βˆ‡π‘˜π‘£β€–β€–ξ‚πœ”π‘‡0β‰€π‘™β‰€π‘˜β‰€1,(B.2) where the constant 𝐢 is independent of β„Ž, 𝑣, and where βˆ‡0𝑣=𝑣. The plan is to extend (B.2) for 𝑙=0 to an inequality involving the form ℬ𝑇=ξ‚΅ξ€œξ‚πœ”π‘‡π‘(πœ‰βˆ’πœ“)2+Μƒπ‘˜2(βˆ‡πœ“)2+Μƒπ‘˜1(βˆ‡πœ‰)2ξ‚Ά1/2,(B.3) which involves the coupling term as well as other coefficients of the problem and thus leads to robustness.
If Μƒπ‘˜2>0, let πœ‰,πœ“βˆˆπ»1(Ξ©) and recall that π‘‰π‘›β„Ž2βŠ†π‘‰π‘›β„Ž1. If Μƒπ‘˜2≑0, we consider πœ‰βˆˆπ»1𝐷(Ξ©), πœ“βˆˆπΏ2(Ξ©), and away from the boundary on some ∘Ω⊊Ω we have π‘‰π‘›β„Ž2(∘Ω)βŠ†π‘‰π‘›β„Ž1(∘Ω). For a fixed 𝑛, 𝑖, let πœ‰β„Žπ‘–, πœ“β„Žπ‘– be the quasi-interpolators of πœ™, πœ“ in π‘‰π‘›β„Žπ‘–, respectively. Note that πœ“β„Ž2 is defined correctly when π‘˜2=0 and πœ“βˆˆπΏ2(Ξ©). To get the desired estimate, we apply (B.2) with π‘‡βˆˆπ’―π‘›β„Ž1, and 𝑙=0, π‘˜=1, to get β€–πœ‰βˆ’πœ‰β„Ž1β€–π‘‡β‰€β„Žπ‘‡β€–βˆ‡πœ‰β€–ξ‚πœ”π‘‡πœ‰β€–ξ‚πœ”π‘‡β‰€(β„Žπ‘‡/Μƒπ‘˜11/2)ℬ𝑇. To extend, we add and subtract πœ“βˆ’πœ“β„Ž1 and use the triangle inequality followed by (B.2) with π‘˜=𝑙=0 and π‘˜=1,𝑙=0 to get β€–β€–πœ‰βˆ’πœ‰β„Ž1β€–β€–π‘‡β‰€β€–β€–ξ€·πœ‰βˆ’πœ‰β„Ž1ξ€Έβˆ’ξ€·πœ“βˆ’πœ“β„Ž1‖‖𝑇+β€–β€–πœ“βˆ’πœ“β„Ž1β€–β€–π‘‡β‰€β€–πœ‰βˆ’πœ“β€–ξ‚πœ”π‘‡+β„Žπ‘‡β€–βˆ‡πœ“β€–ξ‚πœ”π‘‡β‰€βˆšξƒ―12max𝑐1/2,β„Žπ‘‡Μƒπ‘˜21/2ℬ𝑇.(B.4) Thus we have by (4.11) β€–β€–πœ‰βˆ’πœ‰β„Ž1β€–β€–π‘‡β‰€πœƒ1,𝑇ℬ𝑇,(B.5) which is the desired element estimate.
To prove the edge interpolation estimate for an edge 𝐸 of the element 𝑇, we use the following trace inequality [43, Lemma  3.1, page 645] for π‘£βˆˆπ»1(𝑇): β€–π‘£β€–πΈξ€·β„Žβ‰€πΆπ‘‡βˆ’1/2‖𝑣‖𝑇+‖𝑣‖𝑇1/2β€–βˆ‡π‘£β€–π‘‡1/2ξ€Έ,(B.6) where 𝐢 is a constant independent of 𝑣, β„Žπ‘‡. Let 𝑣=πœ‰βˆ’πœ‰β„Ž1. By (B.2) β€–β€–βˆ‡ξ€·πœ‰βˆ’πœ‰β„Ž1ξ€Έβ€–β€–π‘‡β‰€β€–βˆ‡πœ‰β€–ξ‚πœ”π‘‡β‰€β„Žπ‘‡βˆ’1πœƒ1,𝑇ℬ𝑇.(B.7) Applying (B.6) followed by (B.5) and (B.7) we get β€–β€–πœ‰βˆ’πœ‰β„Ž1β€–β€–πΈβ‰€β„Žπ‘‡βˆ’1/2β€–β€–πœ‰βˆ’πœ‰β„Ž1‖‖𝑇+β€–β€–πœ‰βˆ’πœ‰β„Ž1‖‖𝑇1/2β€–β€–βˆ‡πœ‰βˆ’πœ‰β„Ž1β€–β€–1/2β‰€β„Žπ‘‡βˆ’1/2πœƒ1,𝑇ℬ𝑇+πœƒ1/21,𝑇ℬ𝑇1/2β„Žπ‘‡βˆ’1/2πœƒ1/21,𝑇ℬ𝑇1/2.(B.8) Thus we obtain β€–πœ‰βˆ’πœ‰β„Ž1‖𝐸≀𝛾1,𝐸ℬ𝑇, and the estimate for the second component follows similarly if Μƒπ‘˜2β‰’0.
When Μƒπ‘˜2≑0, we have πœ‰βˆˆπ»1𝐷(Ξ©),πœ“βˆˆπΏ2(Ξ©), and (4.12) is reduced to πœƒ2,𝑇=max{1/𝑐1/2,β„Žπ‘‡/Μƒπ‘˜11/2}. Also, 𝛾2,𝑇 does not need to be defined since 𝑅𝑛2,𝐸≑0. Now we add and subtract πœ“βˆ’πœ“β„Ž1 and use the triangle inequality followed by (B.2) with π‘˜=𝑙=0 and π‘˜=1, 𝑙=0 to get β€–β€–πœ“βˆ’πœ“β„Ž2β€–β€–π‘‡β‰€β€–β€–ξ€·πœ‰βˆ’πœ‰β„Ž1ξ€Έβˆ’ξ€·πœ“βˆ’πœ“β„Ž1‖‖𝑇+β€–β€–πœ‰βˆ’πœ‰β„Ž1β€–β€–π‘‡β‰€β€–πœ‰βˆ’πœ“β€–ξ‚πœ”π‘‡+β„Žπ‘‡β€–βˆ‡πœ‰β€–ξ‚πœ”π‘‡β‰€βˆšξƒ―12max𝑐1/2,β„Žπ‘‡Μƒπ‘˜11/2ℬ𝑇.(B.9) Thus, by (4.12) we have β€–πœ“βˆ’πœ“β„Ž1β€–π‘‡β‰€πœƒ2,𝑇ℬ𝑇, which is the desired edge estimate.
Combining these interpolation estimates with the usual technique of upper bounds for residual a posteriori error estimators extended to systems in [13] completes the proof the theorem.


The authors would like to thank Professor Ralph Showalter from Oregon State University for numerous discussions on the double-diffusion models. they would like also to acknowledge the package 𝑖FEM [44] which was used to refine and coarsen the grid in our adaptive implementation. they are also grateful to the anonymous referee whose remarks helped to improve the paper. This work was partially supported by the Grants NSF DMS-0511190, DMS-1115827, DOE 98089, and by FAPESP/Brazil.