Multiscale model reduction for fluid infiltration simulation through dual-continuum porous media with localized uncertainties

Accepted Manuscript Multiscale model reduction for fluid infiltration simulation through dual-continuum porous media with localized uncertainties Qiuqi Li, Yuhe Wang, Maria Vasilyeva

PII: DOI: Reference:

S0377-0427(18)30008-6 https://doi.org/10.1016/j.cam.2017.12.040 CAM 11459

To appear in:

Journal of Computational and Applied Mathematics

Received date : 10 August 2017 Revised date : 20 November 2017 Please cite this article as: Q. Li, Y. Wang, M. Vasilyeva, Multiscale model reduction for fluid infiltration simulation through dual-continuum porous media with localized uncertainties, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2017.12.040 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Manuscript Click here to view linked References

Multiscale Model Reduction for Fluid Infiltration Simulation through Dual-Continuum Porous Media with Localized Uncertainties Qiuqi Li∗ Yuhe Wang† Maria Vasilyeva‡

ABSTRACT Here, we present some Reduced Basis (RB) methods for fluid infiltration problems through certain porous media modeled as dual-continuum with localized uncertainties. We apply dimension reduction techniques to construct a reduced order model. In the RB methods, to perform the offline-online computation decomposition, the model inputs need to be affinely dependent on the uncertainties. We develop a Proper Orthogonal Decomposition and Greedy (POD-Greedy) RB method for stochastic dual-continuum models. In the POD-Greedy RB framework, for heterogeneous porous media, we need to solve the stochastic dual-continuum models many times using very fine grid to construct a set of snapshots for building optimal reduced basis. This offline computation may be very expensive. To improve the offline computational efficiency, we further develop a local global RB method, which integrates the coupled multiscale and multicontinuum approach using Generalized Multiscale Finite Element Method (GMsFEM) to the POD-Greedy RB method. To illustrate the efficiency of the proposed methods, we present two numerical examples for stochastic dual-continuum models. Our numerical results show that both the POD-Greedy RB method and the local global RB method greatly improve the computation efficiency with high approximation accuracy. keywords: Reduced basis methods, GMsFEM, Local global RB method, Greedy algorithm, Proper orthogonal decomposition, Stochastic Dual-continuum model

1

Introduction

Fluid infiltration in porous media has a wide range of scientific and engineering applications [18, 22, 6, 23]. Typical examples include the fluid flow through natural strata (e.g. groundwater and hydrocarbon) and composite materials (e.g. membrane and catalysis) [43, 38, 12, 27, 26]. Very often, there are some degree of fissuring developed in porous media, participially in natural strata presented in the format of fractures [3, 28]. For those fractured porous media, the dual-continuum model after Barenblatt et. al [4] has been the standard approach, based on a very insightful concept and reasonable simplification. The concept of dual-continuum model comes into place by considering dual interconnected parallel systems that distribute all over the domain at any point ∗

College of Mathematics and Econometrics, Hunan University, Changsha 410082, China. Email: [email protected]. † Department of Petroleum Engineering, Texas A&M University at Qatar, Doha, Qatar. Email: [email protected] ‡ Institute for Scientific Computation, Texas A&M University, College Station, TX 77843-3368 & Department of Computational Technologies, North-Eastern Federal University, Yakutsk, Russia. Email: [email protected]

1

in the flow field. In a standard dual-continuum model, mass transfers are coupled between two different continua of distinct scales. One continuum is normally called fracture continuum with apparently high permeability but small volume and the other is the matrix continuum with low permeability but large volume. The coupling between fracture and matrix continua are realized using a transfer function quantifying the mass transfer in between. A terms called shape factor is introduced as a characteristic scale representing the existence and significance of fracture or matrix. Although this dual-continuum approach presents an efficient conceptual way of modeling fluid infiltration through fractured porous media, the shape factors are generally uncertain and may have a big impact on simulation results [32, 33, 2, 42]. Following the introduction of dual-continuum, subsequent extensions and developments have been carried out in the past half century. Wu and Pruess extended it using multiple interacting continua to capture the transient characteristics of inter-continuum flow [48]. Recently, Yan et al. [50] generalized this approach by proposing a generalized multi-continuum model, which can handle any number of continua with arbitrary connections among them. On the other hand, multiscale approaches have been proposed to couple the multiple model equations for multiple continua [9, 1, 10, 44, 45] using Generalized Multiscale Finite Element Method (GMsFEM)[15]. This coupling brings many advantages. The most significant one is its ability to systematically and robustly consider the multiscale heterogeneity of each continuum, such that one can incorporate the heterogenous permeability field quite easily. Though one can use the multiscale approach to consider the model parameter heterogeneity, the ways of treating shape factors, which are generally uncertain, are still unimproved. Generally speaking, uncertain parameters are usually parameterized by a few random variables. In this way, one can describe the models using Parameterized Partial Differential Equations (PPDEs) [34, 49]. It is thus crucially important to solve those PPDEs in an efficient and accurate fashion. Thus, some type of model reduction are required [19]. In this work, we focus on developing the reduced basis approximations based on the coupled GMsFEM method for dual-continuum models with localized uncertainties. Those localized uncertainties refer to the localized uncertain shape factors, which will be detailed in Section 2. Recent works have been devoted to the propagation of the localized uncertainties for stochastic partial differential equation [13, 35]. Here we use a different method to deal with the localized uncertainties. Method of snapshots [24] is used to globalize the localized uncertainties. Reduced Basis (RB) approach has been applied to solve PPDEs in a lower dimensional manifold [36, 41, 7, 40, 21, 5, 39, 20]. RB method constructs a basis functions of much smaller dimension using a set of snapshots, which are the solutions of those PPDEs for certain parameters chosen by some sampling methods. Then RB projects the solution space onto a low-dimensional space, spanned by the set of basis functions. The efficiency is gained by applying a offline-online computation decomposition. We compute the snapshots and compute the reduced basis functions during the offline phase. Then, we solve the reduced model for many cases of parameters in the online phase to evaluate the uncertainty. To obtain the snapshots, one can solve the fine model in a fine grid using some conventional numerical methods such as finite volume method, finite element method and their variations [46, 29]. When model parameters are highly heterogeneous with multiscale features such as porous media with highly heterogeneous permeability field and fractures, one have to use very fine grid to resolve the necessary heterogeneity. This could be very computational demanding, since one has to run the fine model many times to generate the snapshots. To overcome this limitation, we apply the Generalized Multiscale Finite Element Method (GMsFEM) to compute the snapshots and develop the local-global reduced basis methods. In particular, we develop a Local Global RB method for the fluid infiltration simulation through dual-continuum porous media with localized uncertainties. We integrate the coupled GMsFEM [10] and multicontinuum approach to the POD-Greedy RB 2

method. This paper is organized as follows. In Section 2, we present the dual-continuum model, including the variational formation and multiple patches with independent and localized uncertainties about exchange term coefficients. In Section 3, the Reduced Basis (RB) method for dual-continuum problem is described. In Section 4, we introduce the POD-Greedy algorithm to construct the reduced basis functions for the stochastic dual-continuum model. In Section 5, we describe the coupled GMsFEM and Local-global reduced basis method. The method of snapshots is given in Section 6 with numerical examples provided in Section 7. Finally we summarize and conclude in Section 8.

2

Dual-continuum model with localized uncertainty

We first describe briefly the dual-continuum model for fractured porous media. We refer to Barenblatt et. al [4] for the details of this conceptual model. Using the concept of dual-continuum, we assume the porous media can be decomposed into two continua, namely fractures and matrix. This basic model introduces two fluid pressures at the same spatial point, representing the fluid pressures in fracture and matrix. In addition, the model takes into consideration the transfer of fluid between the fracture and matrix media using an exchange term. Let domain D be divided into the fracture domain Df and the matrix domain Dm , where f and m represent the fracture and the matrix regions in the following text. Generally, the matrix domain has a large volume but with small permeability, while the fracture domain has very small volume but with high permeability. We denote the permeability of the fracture region by kf , and denote the permeability of the matrix region by km . Similarly, we use ϕf for the porosity of the fracture domain and ϕm for the porosity of the matrix domain. We consider a dual-continuum model with localized uncertainties in the exchange term coefficient. Thus, the mass balance equations describing this coupled dual-continuum system read kf ∂pf = div( ∇pf ) + Q(σ)(pm − pf ), ∂t µ ∂pm km ϕm c = div( ∇pm ) + Q(σ)(pf − pm ) ∂t µ ϕf c

(2.1)

with boundary conditions on ∂D = ∂DD ∪ ∂DN given by pf = gf , pm = gm on ∂DD , kf ∂pf km ∂pm = hf , = hm , on ∂DN , µ ∂⃗n µ ∂⃗n where ⃗n is the outward-pointing normal to the boundary. Here σ = (σ1 , σ2 , · · · , σr ) denotes localized uncertainties in the exchange term coefficient, Q(σ)(pm , pf ) is the exchange term, c is the total compressibility, pf is the pressure in the fracture region, and pm is the pressure in the matrix region, µ is the fluid viscosity.

2.1

Variational formulation of the problem

System (2.1) leads to the variational problem as follows: for any test functions v, ∂pf , v) + af (pf , v) − q(pm − pf , v) = sf (v) ∂t ∂pm mm ( , v) + am (pm , v) − q(pf − pm , v) = sm (v), ∂t mf (

3

(2.2)

Figure 2.1: Patches Λ = Λ1 ∪ · · · Λ4 ⊆ D containing localized uncertainties. where

∫ ∂pf ∂pf , v) = vdx, mf ( ϕf c ∂t ∂t D ∫ kf ∇pf ∇vdx, af (pf , v) = D µ ∫ q(pm − pf , v) = Q(pm − pf )vdx,

∫ ∂pm ∂pm mm ( , v) = vdx, ϕm c ∂t ∂t D ∫ km am (pm , v) = ∇pm ∇vdx, D µ ∫ q(pf − pm , v) = Q(pf − pm )vdx,

D

and

sf (v) =

∫

∂DN

hf vdx −

∫

D

D

kf ∇Gf ∇vdx, sm (v) = µ

∫

∂DN

hm vdx −

∫

D

km ∇Gm ∇vdx. µ

Here Gf and Gm are continuous piecewise polynomials, and have the same values as gf and gm , respectively, at the boundary nodes belonging to ∂DN .

2.2

Multiple patches with independent localized uncertainties

The exchange term coefficient is related to the underlying heterogeneous reservoir properties, such as matrix permeability, fracture conductivity and fracture density. In practice, we often divide the domain into a number of zones with different attributes, often according to the geological facies of the reservoir. Thus, based on this concept, we consider the Q in equation (2.1) is uncertain only on some parts of the domain Λ = Λ1 ∪ · · · Λr . Each Λi (1 ≤ i ≤ r) is disjoint called as a patch, such that the domain D can be considered as (see figure 2.1) D = (D \ Λ) ⊔ Λ. We then consider Q such that Q(x, σ) =

{

Q0 (x), for x ∈ D \ Λ,

Q(x, σi ), for x ∈ Λi , 1 ≤ i ≤ r,

(2.3)

where {σ1 , · · · , σr } are mutually independent. Here Q0 (x) is a function defined on D \ Λ, and is independent on the uncertainties σ. 4

To fulfill an offline-online computation decomposition, we need to use the affine assumption, i.e., ∑nq i q(w, v; σ) ≈ i=1 ι (σ)q i (w, v), ∀σ ∈ Ω, (2.4)

for the bilinear form q(·, ·; σ) with localized uncertainties. We apply method of snapshots (MOS) [24] to construct a Karhunen-Love expansion for Q to obtain the bilinear form q(·, ·; σ). The MOS will be discussed in Section 6, which is based on Karhunen-Love expansion and can extract most important information from a set of snapshots.

3

Reduced basis method (RBM)

In this section, we follow the framework [25, 40, 41, 30, 31] to present the reduced basis method for the dual-continuum model. Let X be a given finite element (FE) space in a fine grid and dim(X) = NF . Then the FE method of problem (2.2) reads: for any σ ∈ Ω, we find pf (σ) ∈ X and pm (σ) ∈ X such that 1 mf (pnf − pn−1 , v) + af (pnf , v) − q(pnm − pnf , v) = sf (v) f △t 1 n n n mm (pnm − pn−1 m , v) + am (pm , v) − q(pf − pm , v) = sm (v). △t

(3.5)

The backward Euler temporal discretization is used with time step △t. Here n is the number of Nf time steps. Let {Xfn }n=1 be sequence of finite dimensional reduced spaces satisfying Nf

Xf1 ⊂ Xf2 ⊂ · · · ⊂ Xf

⊂ X,

n }Nm be sequence of finite dimensional reduced spaces satisfying and {Xm n=1 1 2 Nm Xm ⊂ Xm ⊂ · · · ⊂ Xm ⊂ X. N

Nm , respectively. RBM is devoted to We denote by Nf and Nm the dimensions of Xf f and Xm approximating the solution pf (σ) and pm (σ) of the problem (2.1) by a few pre-computed solutions pf (σ i ) and pm (σ i ) of (2.2) for some selected values σ i , i = 1, · · ·, N . In a particular case, we take Nf = Nm = N .

3.1

Galerkin projection and offline-online computation

Now we consider the reduced order model. The approximation of RBM is then defined as the Galerkin projection onto these low-dimensional subspaces: for any given σ ∈ Ω, we find pNf (σ) ∈ Nf

Xf

Nm such that and pNm (σ) ∈ Xm

1 Nf n n n mf (pnNf − pn−1 Nf , v) + af (pNf , v) − q(pNm − pNf , v) = sf (v), v ∈ Xf , △t 1 n−1 Nm , v) + am (pnNm , v) − q(pnNf − pnNm , v) = sm (v), v ∈ Xm . mm (pnNm − pN m △t Nf

N

f m and {ψim }N Let {ψif }i=1 i=1 be the set of basis functions of Xf n n solutions pNf and pNm can be represented by

pnNf

=

Nf ∑

f pi,n Nf (σ)ψi ,

pnNm

=

Nm ∑ r=1

i=1

5

(3.6)

Nm , respectively. The and Xm

m pr,n Nm (σ)ψr .

(3.7)

By plugging v = ψj into (3.6), we have Nf

Nf

∑ 1 ∑ i,n−1 mf (ψif , ψjf )(pi,n − p ) + af (ψif , ψjf )pi,n Nf Nf Nf − △t i=1

i=1

Nm ∑

Nf

q(ψrm , ψjf )pr,n Nm +

r=1

1 △t

Nm ∑ r=1

∑

f q(ψif , ψjf )pi,n Nf = sf (ψj ),

i=1

r,n−1 mm (ψrm , ψlm )(pr,n Nm − pNm ) + Nf

∑

q(ψif , ψlm )pi,n Nf +

Nm ∑

Nm ∑ r=1

am (ψrm , ψlm )pr,n Nm −

(3.8)

m q(ψrm , ψlm )pr,n Nm = sm (ψl ),

r=1

i=1

1 ≤ i, j ≤ Nf ,

1 ≤ r, l ≤ Nm ,

Equation (3.8) implies a linear algebraic system with Nf + Nm unknowns. It involves the Nf m computation of inner products with entities {ψif }i=1 and {ψim }N i=1 , each of which is represented by fine grid finite element basis functions of X. This will lead to substantial computation for the inputoutput evaluation σ −→ {pnNf (σ), pnNm (σ)}. With the assumption (2.4) of affine decomposition, equation (3.8) can be rewritten by Nf

Nf

∑ 1 ∑ i,n−1 mf (ψif , ψjf )(pi,n − p ) + af (ψif , ψjf )pi,n Nf Nf Nf − △t i=1

Nm ∑ r=1

nq

(∑ q=1

nq

i=1

q=1

∑(∑

nq

i=1

q=1

) ) ∑(∑ q f ι (σ)q q (ψif , ψjf )pi,n ιq (σ)q q (ψrm , ψjf )pr,n Nf = sf (ψj ), Nm + 1 △t

Nf

i=1

Nf

Nm ∑

mm (ψrm , ψlm )(pr,n Nm

r=1

) ιq (σ)q q (ψif , ψlm )pi,n Nf +

Nm ∑ r=1

−

nq

(∑ q=1

pr,n−1 Nm )

+

Nm ∑ r=1

am (ψrm , ψlm )pr,n Nm −

(3.9)

) m ιq (σ)q q (ψrm , ψlm )pr,n Nm = sm (ψl ), 1 ≤ i, j ≤ Nf ,

1 ≤ r, l ≤ Nm ,

This gives rise to the matrix form  nq nq ∑ ∑ − ιq (σ)QqNmf MNf + ANf + ιq (σ)QqNf  q=1 q=1  nq nq  ∑ ∑ − ιq (σ)QqNf m MNm + ANm + ιq (σ)QqNm q=1

q=1

where (MNf )ij = (MNm )rl =

1 mf (ψif , ψjf ), △t

1 mm (ψrm , ψlm ), △t

(QqNmf )rj = q q (ψrm , ψjf ),



[ ] [ ]  pnN MNf pn−1 Nf + sNf f  = ,  pn MNm pn−1 Nm Nm + sNm

(ANf )ij = af (ψif , ψjf ), (ANm )rl = am (ψrm , ψlm ), (QqNf m )il = q q (ψif , ψlm ), 6

(QqNf )ij = q q (ψif , ψjf ), (QqNm )rl = q q (ψrm , ψlm ), (sNf )j = sf (ψjf ),

(3.10)

and (sNm )l = sm (ψlm ), N

1 ≤ i, j ≤ Nf ,

1 ≤ r, l ≤ Nm , N

1 ≤ q ≤ nq .

f f m Let {ξk }k=1 be the FE basis of X. Because basis functions {ψif }i=1 and {ψim }N i=1 belong to the FE space X, they can be written as

ψif =

Nf ∑

f Zki ξk ,

k=1

1 ≤ i ≤ Nf .

ψrm =

Nm ∑

m Zkr ξk ,

k=1

1 ≤ r ≤ Nm .

f m , 1 ≤ r ≤ N . Thus we can get Let (Z f )ki = Zki , 1 ≤ i ≤ Nf , (Z m )kr = Zkr m

MNf = MNm = and

1 (Z f )T Mf Z f , △t

1 (Z m )T Mm Z m , △t

QqNmf = (Z m )T Qq Z f ,

ANf = (Z f )T Af Z f , ANm = (Z m )T Am Z m ,

QqNf m = (Z f )T Qq Z m ,

QqNf = (Z f )T Qq Z f , QqNm = (Z m )T Qq Z m ,

sNf = (Z f )T Ff ,

sNm = (Z m )T Fm .

where (Mf )ij = mf (ξj , ξi ), (Mm )ij = mm (ξj , ξi ), (Af )ij = af (ξj , ξi ), (Am )ij = am (ξj , ξi ), (Qq )ij = q q (ξj , ξi ), and (Ff )i = sf (ξi ), (Fm )i = sm (ξi ), 1 ≤ i, j ≤ NF . The matrixes MNf , ANf , QqNf , MNm , ANm , QqNm , QqNmf , QqNf m , and the vectors sNm , sNf are independent of localized uncertainties σ. As a result, their computation is once and could be done in offline phase. The online computation is then to solve equation (3.10) for any σ ∈ Ω. Because the online computation only involves Nf + Nm unknowns (Nf + Nm ≪ 2NF ), one can obtain significant efficiency by this setting.

4

POD-greedy algorithm for sampling reduced basis

To find a few optimal reduced basis and to assure the fidelity of the reduced model to approximate the original model, we combine the POD in t with the greedy procedure in σ (POD-greedy algorithm) to obtain reduce basis functions for the dual-continuum stochastic models. To this end, we need to make the posteriori error bounds for the dual-continuum stochastic models. We first divide the time interval I into J = T /△t subintervals of equal length △t and define tn = n△t, 0 ≤ n ≤ J. We consider the residual error for dual-continuum stochastic models. Let enm (σ) := pnm (σ) − pnNm (σ).

enf (σ) := pnf (σ) − pnNf (σ), By equation (3.5), we get: ∀v ∈ X,

1 n mf (enf (σ) − en−1 (σ), v) + af (enf (σ), v) − q(enm (σ) − enf (σ), v) = rN (v; σ), f f △t 1 n−1 n mm (enm (σ) − em (σ), v) + am (enm (σ), v) − q(enf (σ) − enm (σ), v) = rN (v; σ). m △t where 1 Nf Nf Nm mf (pnNf − pn−1 Nf , v) − af (pf , v) + q(pm − pf , v), △t 1 Nf n Nm Nm rN (v; σ) = sm (v) − mm (pnNm − pn−1 Nm , v) − am (pm , v) + q(pf − pm , v), m △t

n rN (v; σ) = sf (v) − f

7

(4.11)

n (v; σ) and r n (v; σ) are the residuals associated with the reduced basis approximation. and rN Nm f By Riesz representation theory, there exists eˆnNf (σ) ∈ X and eˆnNm (σ) ∈ X such that: ∀v ∈ X,

(

(

eˆnNf (σ), v

)

eˆnNm (σ), v

X

n = rN (v; σ) f

X

n = rN (v; σ). m

)

(4.12)

n (v; σ) and r n (v; σ) can be evaluated through the Consequently, the dual norms of the residual rN Nm f Riesz representation, n ∥rN (v; σ)∥X ∗ := sup f n ∥rN (v; σ)∥X ∗ m

n (v; σ) rN f

= ∥ˆ enNf (σ)∥X , ∥v∥ X v∈X n (v; σ) rN := sup m = ∥ˆ enNm (σ)∥X . ∥v∥ X v∈X

(4.13)

The computation of the residual is crucial to perform the offline-online computation decomposition. We define an error estimator [37, 41] for the solution as ∆N (σ) :=

n √ ∑ ekNm (σ)∥2X △t. ∥ˆ ekNf (σ)∥2X + ∥ˆ

(4.14)

k=1

n (v; σ) and r n (v; σ), respectively, We note that the eˆnNf (σ) ∈ X and eˆnNm (σ) ∈ X are related to rN Nm f through the equation (4.13). By (4.11), (3.7) and (2.4), the residual can be expressed as n rN (v; σ) =sf (v) − f

1 Nf Nf Nm mf (pnNf − pn−1 Nf , v) − af (pf , v) + q(pm − pf , v), △t Nf

=sf (v) − Nm ∑ r=1

r=1

i=1

i=1

Nf

q(ψrm , v)pr,n Nm − Nf

=sf (v) − Nm ∑

Nf

∑ 1 ∑ i,n−1 mf (ψif , v)(pi,n af (ψif , v)pi,n Nf − pNf ) − Nf + △t

nq

∑

q(ψif , v)pi,n Nf

(4.15)

i=1

Nf

∑ 1 ∑ i,n−1 mf (ψif , v)(pi,n af (ψif , v)pi,n Nf − pNf ) − Nf + △t

(∑ q=1

i=1

i=1

Nf

nq

i=1

q=1

) ∑(∑ q ) ιq (σ)q q (ψrm , v)pr,n ι (σ)q q (ψif , v)pi,n Nm − Nf ,

8

1 Nf Nm Nm mm (pnNm − pn−1 Nm , v) − am (pm , v) + q(pf − pm , v), △t Nm Nm ∑ 1 ∑ r,n r,n−1 m =sm (v) − mm (ψr , v)(pNm − pNm ) − am (ψrm , v)pr,n Nm + △t

n rN (v; σ) =sm (v) − m

r=1

Nf

∑ i=1

q(ψif , v)pi,n Nf − N

r=1

Nm ∑

q(ψrm , v)pr,n Nm

(4.16)

r=1

N

m m ∑ 1 ∑ r,n−1 mm (ψrm , v)(pr,n − p ) − am (ψrm , v)pr,n Nm Nm Nm + △t

=sm (v) −

r=1

Nf

nq

i=1

q=1

∑(∑

r=1

) ιq (σ)q q (ψif , v)pi,n Nf −

Nm ∑ r=1

nq

(∑ q=1

) ιq (σ)q q (ψrm , v)pr,n Nm .

By (4.15), (4.16) and (4.12), we have (

eˆnNf (σ), v

)

X

Nm ∑ r=1

(

eˆnNm (σ), v

)

Nf

Nf

=sf (v) −

X

∑ 1 ∑ i,n−1 mf (ψif , v)(pi,n af (ψif , v)pi,n Nf − pNf ) − Nf + △t i=1

i=1

nq

(∑ q=1

=sm (v) − Nf

nq

i=1

q=1

nq

i=1

q=1

) ) ∑(∑ q ιq (σ)q q (ψrm , v)pr,n ι (σ)q q (ψif , v)pi,n Nm − Nf ,

1 △t

∑(∑

Nf

Nm ∑ r=1

r,n−1 mm (ψrm , v)(pr,n Nm − pNm ) −

) ιq (σ)q q (ψif , v)pi,n Nf −

Nm ∑ r=1

nq

(∑ q=1

Nm ∑

am (ψrm , v)pr,n Nm +

r=1

) ιq (σ)q q (ψrm , v)pr,n Nm .

This implies that Nf

eˆnNf (σ)

Nf

∑ i,n 1 ∑ i,n i =Cf + (pNf − pi,n−1 )(σ)X + pNf (σ)Lif + f Nf △t i=1

nq

Nm ∑ ∑ ( r=1

q=1

i=1

Nf

nq

i=1

q=1

) ∑(∑ q i,q ) r,q ι (σ)pi,n ιq (σ)pr,n (σ)I m − Nm Nf (σ)If , N

N

m m ∑ 1 ∑ r,n−1 r r (pr,n − p )(σ)X + pr,n eˆnNm (σ) =Cm + m Nm Nm Nm (σ)Lm + △t

r=1

Nf

nq

i=1

q=1

∑(∑

(4.17)

r=1

i,q ) ιq (σ)pi,n − Nf (σ)If

Nm ∑ r=1

nq

(∑ q=1

) r,q ιq (σ)pr,n Nm (σ)Im .

where Cf is the Riesz representation of sf v, i.e., (Cf , v)X = sf (v) ; Cm is the Riesz representation of sm v, i.e., (Cm , v)X = sm (v) ; Xfi is the Riesz representation of mf (ψif , v), i.e.,(Xfi , v)X =

−mf (ψif , v); Lif is the Riesz representation of af (ψif , v), i.e., (Lif , v)X = −af (ψif , v); Ifi,q is

r is the Riesz representation the Riesz representation of q q (ψif , v), i.e., (Ifi,q , v)X = q q (ψif , v); Xm r m r m of mm (ψr , v), i.e.,(Xm , v)X = −mm (ψr , v); Lm is the Riesz representation of am (ψrm , v), i.e.,

9

r,q r,q (Lrm , v)X = −am (ψrm , v); Im is the Riesz representation of q q (ψrm , v), i.e., (Im , v)X = q q (ψrm , v) for any v ∈ X, where 1 ≤ i ≤ Nf , 1 ≤ r ≤ Nm and 1 ≤ q ≤ nq . The equation (4.17) gives rise to Nf

∥ˆ enNf (σ)∥2X

Nf

∑ i,n 2 ∑ i,n i =(Cf , Cf )X + (pNf − pi,n−1 pNf (σ)(Cf , Lif )X + Nf )(σ)(Cf , Xf )X + 2 △t i=1

2

nq Nm ( ∑ ∑ r=1

q

ι

q=1

i=1

r,q (σ)pr,n Nm (σ)(Cf , Im )X

)

nq ∑(∑ Nf

−2

q=1

i=1

q

ι

i,q (σ)pi,n Nf (σ)(Cf , If )X

) +

Nf Nf 1 ∑ ∑ i,n ′ i,n−1 i′ ,n i′ ,n−1 (p − p )(σ)(p − p )(σ)(Xfi , Xfi )X + N N N Nf f f f △t2 ′ i=1 i =1 Nf Nf

′ ′ 2 ∑ ∑ i,n ′ pNf (σ)(piN,n − piN,n−1 )(σ)(Xfi , Lif )X + f f △t ′

i=1 i =1

2 △t

Nf nq Nm ∑ ∑ ∑ r=1 i′ =1 q=1 Nf Nf

′

′

′

i r,q ιq (σ)(piN,n − piN,n−1 )(σ)pr,n Nm (σ)(Xf , Im )X − f f

nq

′ ′ 2 ∑∑∑ q i,q i′ ι (σ)(piN,n − piN,n−1 )(σ)pi,n Nf (σ)(Xf , If )X + f f △t ′

(4.18)

i=1 i =1 q=1

Nf Nf ∑ ∑

i=1 i′ =1

2

′

Nf Nf nq ∑ ∑∑

′

Nf nq Nm ∑ ∑ ∑ r=1 i′ =1 q=1

′

∑∑∑∑

′

r=1 r′ =1 q=1 q ′ =1

i=1 r′ =1 q=1 q ′ =1

i=1 i′ =1 q=1 q ′ =1

′

′

r ,n q r ,q r,q ιq (σ)pr,n Nm (σ)ι (σ)pNm (σ)(Im , Im )X −

Nf Nm nq nq ∑ ∑∑∑

Nf Nf nq nq ∑ ∑∑∑

′

′

′

′

′

′

r ,n i,q q r ,q ιq (σ)pi,n Nf (σ)ι (σ)pNm (σ)(Im , If )X + ′

′

′

′

i ,n i ,q i,q q ιq (σ)pi,n Nf (σ)ι (σ)pNf (σ)(If , If )X ,

10

′

i r,q ιq (σ)piN,n (σ)pr,n Nm (σ)(Lf , Im )X − f

i,q i ιq (σ)piN,n (σ)pi,n Nf (σ)(Lf , If )X + f

i=1 i′ =1 q=1 Nm Nm nq nq

2

′

i ,n i i pi,n Nf (σ)pNf (σ)(Lf , Lf )X + 2

and N

∥ˆ enNm (σ)∥2X

N

m m ∑ 2 ∑ r,n r,n−1 r r (pNm − pNm )(σ)(Cm , Xm )X + 2 pr,n =(Cm , Cm )X + Nm (σ)(Cm , Lm )X + △t

r=1

nq ∑(∑

r=1

Nf

2

i=1

q

ι

q=1

i,q (σ)pi,n Nf (σ)(Cm , If )X

)

−2

nq Nm ( ∑ ∑ r=1

q

ι

q=1

r,q (σ)pr,n Nm (σ)(Cm , Im )X

) +

Nm ∑ Nm 1 ∑ r′ ,n r′ ,n−1 r,n−1 r,n r′ r )X + , Xm )(σ)(p − p )(σ)(Xm − p (p Nm N N N m m m △t2 ′

2 △t

r=1 r =1 N Nm m ∑ ∑ r=1 r′ =1 Nf N

′

′

′

r ,n r ,n−1 r pr,n )(σ)(Xm , Lrm )X + Nm (σ)(pNm − pNm nq

m ∑ ′ ,n ′ ,n−1 2 ∑∑ i,q r′ ιq (σ)(prNm − prNm )(σ)pi,n Nf (σ)(Xm , If )X − △t ′

i=1 r =1 q=1 Nm Nm nq

′ ,n ′ ,n−1 2 ∑∑∑ q r′ r,q ι (σ)(prNm − prNm )(σ)pr,n Nm (σ)(Xm , Im )X + △t ′

(4.19)

r=1 r =1 q=1

Nm ∑ Nm ∑

r=1 r′ =1

2

′

nq Nm ∑ Nm ∑ ∑ r=1 r′ =1 q=1

′

Nf Nm nq ∑ ∑∑ i=1 r′ =1 q=1

′

′

,n i,q r ιq (σ)prNm (σ)pi,n Nf (σ)(Lm , If )X −

′

,n r r,q ιq (σ)prNm (σ)pr,n Nm (σ)(Lm , Im )X +

Nf Nf nq nq ∑ ∑∑∑

′

′

′

′

i ,n i ,q i,q q ιq (σ)pi,n Nf (σ)ι (σ)pNf (σ)(If , If )X −

i=1 i′ =1 q=1 q ′ =1

2

′

r ,n r r pr,n Nm (σ)pNm (σ)(Lm , Lm )X + 2

Nf nq nq Nm ∑ ∑ ∑∑

r=1 i′ =1 q=1 q ′ =1

nq nq Nm ∑ Nm ∑ ∑ ∑

r=1 r′ =1 q=1 q ′ =1

′

′

′

′

i ,n i ,q q r,q ιq (σ)pr,n Nm (σ)ι (σ)pNf (σ)(If , Im )X + ′

′

′

′

r ,n q r ,q r,q ιq (σ)pr,n Nm (σ)ι (σ)pNm (σ)(Im , Im )X ,

To efficiently compute ∥ˆ enNm (σ)∥X , ∥ˆ enNf (σ)∥X and ∆N (σ), we apply an offline-online procedure. In the offline stage we compute and store the quantities independent of the uncertainties.In particr,q r , and Lr , and store (C , C ) , (C , X i ) , (C , Li ) , ular, we compute Cf , Cm , Xfi , Lif , Im , Ifi,q , Xm f f X f f m f X f X ′

′

′

′

′

r,q r,q r,q )X , )X , (Xfi , Ifi,q )X , (Lif , Lif )X , (Lif , Im (Cf , Im )X , (Cf , Ifi,q )X , (Xfi , Xfi )X , (Xfi , Lif )X , (Xfi , Im ′

′

′

′

′

′

′

r ,q r,q r ,q r ) , (C , Lr ) , (C , I i,q ) , (Lif , Ifi,q )X , (Im , Im )X , (Im , Ifi,q )X , (Ifi ,q , Ifi,q )X , (Cm , Cm )X , (Cm , Xm m f X m X m X ′

′

′

′

′

′

r,q r , X r ) , (X r , Lr ) , (X r , I i,q ) , (X r , I r,q ) , (Lr , Lr ) , (Lr , I i,q ) , and (Lr , I r,q ) , (Cm , Im )X , (Xm m X m f X m X m m X m f X m X m m m m X ′ ′ ′ where 1 ≤ i, i ≤ Nf , 1 ≤ r, r ≤ Nm , 1 ≤ q, q ≤ nq . In the online stage, for any σ, we compute pr,n Nm , i,n n pNf (1 ≤ i ≤ Nf , 1 ≤ r ≤ Nm and 1 ≤ n ≤ J) and use (4.18)and (4.19) to compute ∥ˆ eNm (σ)∥X , ∥ˆ enNf (σ)∥X , and use (4.14) to compute ∆N (σ). Next we describe the POD-greedy algorithm to construct the reduced basis functions for the dual-continuum stochastic model (2.2). Firstly, we summarize the basic idea of POD, which evaluates the eigenfunctions and eigenvalues of a covariance operator to determine an orthonormal basis

11

and to estimate the approximation quality of the corresponding subspace. Based upon a set of snapshots, the POD is often formulated as an optimization problem as follows: given M elements wj ∈ X, 1 ≤ j ≤ M , the algorithm POD(w1 , w2 , · · · , wM , N ) returns N (≤ M ) orthonormal functions with regard to (·, ·)X . Thus the optimal space of dimension N , i.e., XN = span{ψi , 1 ≤ i ≤ N } is obtained by the following optimization problem XN =

arg inf Y ⊂span{wi ,1≤i≤N }

(

)1/2 M 1 ∑ 2 inf ∥wi − v∥X . v∈Y M i=1

Let Ξtrain be a given training set. Then we describe the POD-greedy algorithm in Algorithm 1. Algorithm 1 POD-greedy algorithm constructing the reduced basis functions for the dualcontinuum stochastic model. Input: A training set Ξtrain ⊂ Γ, σ1 ∈ Ξtrain , Ntf , Ntm , Mtf and Mtm N f Output:The RB approximation spaces Xf f = span{ψ1f , ψ2f , · · · , ψN } f Nm = span{ψ m , ψ m , · · · , ψ m } and Xm 1 2 Nm 1: Initialize N = 1, Xf = Ø, Xm = Ø, Nf = 0, Nm = 0 and σN = σ1 ; 2: Construct the snapshots {pf (tn , σN )}Jn=1 and {pm (tn , σN )}Jn=1 by solving the equation (2.2) with FEM and the backward Euler temporal discretization; 3: Compute the optimal basis functions with regard to different time levels {ϕfi , 1 ≤ i ≤ Mtf } = POD({pf (tj , σN ), 1 ≤ i ≤ J}, Mtf ) m j m and {ϕm i , 1 ≤ i ≤ Mt } = POD({pm (t , σN ), 1 ≤ i ≤ J}, Mt ); f f m 4: Update Xf ← {Xf , {ϕi , 1 ≤ i ≤ Mt }}, Xm ← {Xm , {ϕi , 1 ≤ i ≤ Mtm }}, Nf ← Nf + Ntf and Nm ← Nm + Ntm ; 5: Construct the reduced basis approximation spaces via POD algorithm again, i.e., N f {ψif , 1 ≤ i ≤ Nf } = POD(Xf , Nf ), Xf f = span{ψ1f , ψ2f , · · · , ψN } f

6: 7: 8: 9:

5

Nm = span{ψ m , ψ m , · · · , ψ m }; and {ψim , 1 ≤ i ≤ Nm } = POD(Xm , Nm ), Xm 1 2 Nm Update Ξtrain with Ξtrain ← Ξtrain = Ξtrain \ σN , for each σ ∈ Ξtrain , evaluate the error estimator ∆N (σ) by (4.14). Choose σN +1 = arg maxµ∈Ξtrain ∆N (σ), and set εN = maxµ∈Ξtrain ∆N (σ); N ← N + 1; Return to Step 2 if εN ≤ εN −1 , otherwise terminate.

Local-global reduced basis method

In the reduced basis method, we need to compute a set of snapshots {pf (tn , σN )}Jn=1 and {pm (tn , σN )}Jn=1 to construct reduced basis. These snapshots can be computed by traditional finite element methods in a fine grid and the backward Euler temporal discretization, as discussed in Section 4. If the snapshots have strong multiscale features, then we have to use a very fine mesh to resolve the features in all scales. This computation of computing snapshots may be very expensive. To overcome the difficulty and improve the offline computation, we can use the coupled Generalized Multiscale Finite Element Method (coupled GMsFEM) to compute the snapshots. In this section, we briefly present coupled GMsFEM and Local-global reduced basis method.

12

Figure 5.2: Illustration of a coarse neighborhood and coarse element.

5.1

Coupled GMsFEM

Let KH be a conforming coarse partition for the computational domain D, where H is the coarse mesh size. Each coarse-grid block is further partitioned into a connected union of fine-grid blocks, c we get the fine grid partition Kh . Let εH := {xi }N i=1 (Nc is the number of coarse nodes) be the H vertices of the coarse mesh K . The coarse neighborhood wi corresponding to the coarse edge xi is defined by ∪ wi = {Kj ∈ KH ; xi ∈ ∂Kj }. m Li i See Figure 5.2 for illustration. Let {φfj }L j=1 and {φj }j=1 be the sets of the coupled GMsFE basis functions corresponding to coarse node xi , which are supported in wi . We define the coupled GMsFE spaces as the linear span of all local basis functions, i.e.,

XfH =

⊕ f Lf i , {φj }j=1

H Xm =

εH

⊕ Lm i {φm j }j=1 . εH

H H Thus the coupled GMsFEM is to find {pfH (σ), pm H (σ)} ∈ Xf × Xm such that

1 f,n−1 m,n f,n H mf (pf,n , v) + af (pf,n H − pH H , v) − q(pH − pH , v) = sf (v), v ∈ Xf , △t 1 m,n−1 f,n m,n H mm (pm,n , v) + am (pm,n H − pH H , v) − q(pH − pH , v) = sm (v), v ∈ Xm . △t

(5.20)

We follow GMsFEM [17, 15] and present the construction of the coupled GMsFE spaces XfH H . We first generate the snapshot space by solving local problems subject to some boundary and Xm conditions, and then use local spectral decomposition [8, 15, 14, 16, 11] to obtain a lower-dimensional offline space. Let xi ∈ εH be a coarse node, and Jh (wi ) denotes the fine-grid boundary nodes on ∂wi . We define a piecewise constant function δjl on ∂wi as δjl

=

{

1, j = l, 0, j ∈ Jh (wi ) \ l. 13

(5.21)

We define ϕsnap = (ϕsnap,f , ϕsnap,m ) by solving the following problems on the coarse neighbori,l i,l i,l hood wi corresponding to the edge xi , { 2 a(ϕsnap ¯(ϕsnap i,l , v) − q i,l , v) = 0 ∀v = (v1, v2) ∈ [X(wi )] , (5.22) = δl on ∂wi , ϕsnap i,l where

q¯(ϕsnap i,l , v)

:=

∫

wi

∫

kf km ∇ϕsnap,f ∇v + ∇ϕsnap,m ∇vdx, i,l i,l µ wi µ ∫ snap,f snap,m ¯ ¯ snap,m − ϕsnap,f )v2 dx Q(ϕi,l − ϕi,l )v1 dx + Q(ϕ i,l i,l

ai (ϕsnap i,l , v)

=

wi

¯ being the mean of Q with regard to the localized uncertainties. In the above definition, the with Q discrete delta function δl is defined as δl = (δj,1 , δj,2 ), and each δj,k is the discrete delta function defined as (5.21). Then we define the local snapshot space Vsnap (wi ) by Vsnap (wi ) = span{ϕsnap i,l (σ) : ∀l}, After the construction of the snapshot space, the offline space is constructed by performing some local spectral problem on the snapshot space. More precisely, for each coarse region wi , we consider the following local eigenvalue problem: find the l−th eigenpair {λil , φil } (φil ∈ Vsnap (wi )) such that ai (φil , v) − q¯(φil , v) = λil si (φil , v)

∀ v ∈ Vsnap (wi ),

(5.23)

where the bilinear form si is defined by ∫ kf km i |∇χi |ϕsnap,f v+ |∇χi |ϕsnap,m vdx, s (v, w) = i,l i,l µ wi µ where {χi } is a set of partition of unity functions for the coarse-grid partition of the domain D. Suppose that the eigenvalues of (5.23) are arranged in increasing order, the corresponding i corresponding to the first li eigenvectors are denoted by φil . We take the eigenfunctions {φil }ll=1 eigenvalues for each local eigenvalue problem to get the GMsFEM basis functions, i.e., {ψli := φil χi : 1 ≤ l ≤ li , 1 ≤ i ≤ Nc }. To simplify notation, we use the following single-index notation ∑ c f m {ψj : 1 ≤ j ≤ NH }, where NH = N i=1 li . We note that ψj can be represented as ψj = (ψj , ψj ), H are defined as follows then the coupled GMsFE spaces XfH and Xm XfH = span{φfj : 1 ≤ j ≤ NH } and H Xm = span{φm j : 1 ≤ j ≤ NH }.

We define

T m RfT = [φf1 , · · · , φfNH ], Rm = [φm 1 , · · · , φNH ],

where φfj and φm j are nodal values of each basis function defined on the fine grid. Thus we can rewrite equation (5.20) in a discrete form as  nq nq ∑ ∑ q H q + + A ι (σ)Q − MH ιq (σ)QqH,mf f H,f  f q=1 q=1  nq nq  ∑ ∑ H − ιq (σ)QqH,f m MH ιq (σ)QqH,m m+ Am+ q=1

q=1

14



] [ H n−1 [ ]  pnH,f Mf pH,f + sH,f  = , n−1  pnH,m MH m pH,m + sH,m

(5.24)

where MH f = MH m = and

1 (Rf )T MRf , △t

1 (Rm )T MRm , △t

QqH,mf = (Rm )T Qq Rf ,

5.2

T AH f = (Rf ) Af Rf , T AH m = (Rm ) Am Rm ,

QqH,f m = (Rf )T Qq Rm ,

QqH,f = (Rf )T Qq Rf , QqH,m = (Rm )T Qq Rm ,

sH,f = (Rf )T Ff ,

sH,m = (Rm )T Fm .

Local-global reduced basis method

For the dual-continuum stochastic model (2.2) in heterogeneous random porous media, we use the Local-global reduced basis method to improve the computation efficiency. We outlines the Local-global reduced basis method in Algorithm 2. Algorithm 2 Local-global reduced basis method for the dual-continuum stochastic model Input: The dual-continuum stochastic model (2.2) Output: The local-global reduced basis model under the form (3.10) H and get the coupled GMsFE model ◃:Construct the coupled GMsFE spaces XfH and Xm (5.24); ◃:We use Algorithm 1 with Step 2 being replaced by “Construct the snapshots {pf (tn , σN )}Jn=1 and {pm (tn , σN )}Jn=1 by solving the equation (2.2) with coupled GMsFE and backward Euler temporal discretization, i.e., equation (5.24);” to construct the local-global N f Nm = span{ψ m , ψ m , · · · , ψ m }; RB approximation spaces Xf f = span{ψ1f , ψ2f , · · · , ψN } and Xm 1 2 Nm f ◃: Construct the local-global reduced basis model (3.10); ◃: For random input σ ∈ Ω, we get the solution of (2.2) by solving local-global reduced basis model (3.10), then we can estimate the influence of the uncertainties quickly by (3.10). Remark 5.1. We note that the Local-global reduced basis method is proposed for the dual-continuum stochastic model in highly developed fractured media without long fractures.

6

Method of snapshots

To achieve the offline-online computation decomposition, we need to use the affine assumption (2.4) for the bilinear form q(·, ·; σ), where σ are localized uncertainties defined as in equation (2.3). If q(·, ·; σ) is not affine with respect to σ, then we can use method of snapshots (MOS) to get an affine decomposition. MOS is based on Karhunen-Love expansion and can extract most important information from a set of snapshots. When it is difficult to get the knowledge of a covariance function for the models inputs Q(x, σ), we can use method of snapshots to construct a KL expansion for the models inputs Q(x, σ). Let Ξt be a collection of a finite number of samples in Ω and the cardinality |Ξt | = nt . For ∀ σ ∈ Ξt , we can split G(x, σ) into two parts, i.e., ¯ ˜ σ), Q(x, σ) = Q(x) + Q(x, ∑ t ˜ σ) = Q(x, σ) − Q(x) ¯ ¯ Q(x, σi ) is the mean, and Q(x, is a where Q(x) := E[Q(x, ·)] = n1t ni=1 nt ˜ ˜ random fluctuating part. To obtain Q(x, σ), we take a set of snapshots {Q(x, σi )}i=1 and compute 15

a covariance matrixes C, whose entries can be defined by ( ) 1 ˜ ˜ Cn,m := Q(x, σn ), Q(x, σm ) . nt X Let {λk , ek } be the eigen-pairs (normalized) of C, 1 ≤ k ≤ nt . Set (ek )j = ejk , we define the functions nt 1 ∑ ˜ σj ). gk (x) := √ ej Q(x, λk nt j=1 k

(6.25)

It is easy to get (gk , gl )X = δk,l , 1 ≤ k, l ≤ nt . Then it holds that ˜ σ) ≈ Q(x,

nt √ ∑

ˆ i ιi (σ)gi (x), λ

i=1

t where {ιi (σ)}ni=1 are given by

Thus we get the decomposition

) 1 (˜ ιi (σ) := √ Q(·, σ), gi X . ˆi λ

¯ Q(x, σ) ≈ Q(x) +

7

MQ √ ∑

ˆ i ιi (σ)gi (x). λ

(6.26)

(6.27)

i=1

Numerical results

In this section, we present two numerical examples to illustrate the applicability of the proposed reduced basis methods for solving the dual-continuum stochastic models. In Section 7.1, we consider an example to illustrate performance of the POD-greedy reduced basis method for solving the dualcontinuum stochastic model with 8 dimensional localized uncertainties. In Section 7.2, we use the local global reduced basis method to a dual-continuum stochastic model in highly heterogeneous porous media with 8 dimensional localized uncertainties. Before presenting the individual examples, we describe the computational domain, that is, spatial domain D = (0, 100m)2 , and the other model parameters used are as follows: ◃ ϕf = 0.01, ϕm = 0.1, c = 1.45 × 10−8 , µ = 0.001 [pa · see], ◃ For example in Section 7.1 kf = 9.869233 × 10−13 [m2 ], km = 9.869233 × 10−16 [m2 ]; for example in Section 7.2 kf is shown in Figure 7.4, km = kf /1000, ◃ Q = σkµm , where σ = (σ1 , · · · , σ8 ) are localized uncertainties, that is,  km   , for x ∈ D \ Λ,  µ Q(σ) = σi km    , for x ∈ Λi , 1 ≤ i ≤ 8, µ

where D = (D \ Λ) ⊔ Λ is depicted in Figure 7.3. ◃ Initial conditions p = 3000 × 6894.76 [Pa], ◃ Boundary conditions pf = pm = 2000 × 6894.76 [Pa]. (Here we consider Dirichlet boundary 16

Figure 7.3: Patches Λ = Λ1 ∪ · · · Λ8 ⊆ D containing localized uncertainties. condition) ◃ The final simulation time T = 1500 day. Let (pf , pm ) be the solutions of the proposed reduced basis methods, and (pf,h , pm,h ) be the reference solutions, which are solved by FEM on a fine grid. Then the relative mean errors in the weighted L2 norm are defined by εf =

N 1 ∑ ∥pf (x, σ (i) ) − pf,h (x, σ (i) )∥2L2 , N ∥pf,h (x, σ (i) )∥2L2

(7.28)

i=1

εm

N 1 ∑ ∥pm (x, σ (i) ) − pm,h (x, σ (i) )∥2L2 = , N ∥pm,h (x, σ (i) )∥2L2

(7.29)

i=1

∫ where ∥v∥2L2 = D v 2 dx, and N is the number of samples. The reference solution is computed on a 200 × 200 fine grid. Thus the step for x is h = 0.5m, and the time step is △t = 5 up to the final simulation time T = 1500. We note that each σi obey power law distribution with the minimum value being 1/(3h2 ) and the exponent being 2.5, which make the mean of σi be 1/h2 . To fulfill offline-online computation decomposition, we use MOS for Q(x, σ) and get an affine expression (6.27) with MQ = 8, and the relative mean error by MOS is 3.52 × 10−5 .

7.1

Numerical results for POD-greedy reduced basis method

In this numerical example, we take kf = 9.869233 × 10−13 [m2 ], km = 9.869233 × 10−16 [m2 ], and apply the POD-greedy algorithm described in Section 4 to get the reduced basis model for the dual-continuum stochastic models. We randomly choose nt = |Ξtrain | = 200 samples for localized uncertainties, take Ntf = 4, Ntm = 4, Mtf = 8 and Mtm = 8, and use the Algorithm 1 to construct the reduced basis (RB) approximation spaces. We randomly choose 10000 samples and compute the average relative errors defined as (7.28), (7.29). Figure 7.5 shows the average relative error for the solutions versus different time levels with the number of the reduced basis functions being Nf = 2, 5 and 8, Nm = 2, 5 and 8. By the 17

100 −12 90 −12.5 80 −13 70 −13.5 60 −14 50 −14.5 40 −15 30 −15.5 20 −16 10

−16.5

10

20

30

40

50

60

70

80

90

100

(a) log10 (kf )

Figure 7.4: High-contrast permeability in the numerical experiments.

p

10

2 global basis functions 4 global basis functions 8 global basis functions

0

10

−2

10

−4

10

m

2

2 global basis functions 4 global basis functions 8 global basis functions

0

−2

−4

2

10

10

10

10

10

10

Relative L

Relative L

2

error

10

p

f

2

error

10

−6

−8

10

10

−10

10

−12

0

250

500

750

1000

1250

1500

Time

10

−6

−8

−10

−12

0

250

500

750

1000

1250

1500

Time

Figure 7.5: The average relative error versus different time levels with the number of the reduced basis functions being Nf = 2, 4 and 8, Nm = 2, 4 and 8, the numbers in the legend of Figures are just for pf or pm .

18

figure, we find that the average relative error becomes smaller when the number of the reduced basis functions Nf and Nm increase. To visualize the individual errors of 100 samples, we plot the relative errors for the number of the reduced basis functions being Nf = 2, 4 and 8, Nm = 2, 4 and 8 in Figure 7.6, which shows that (1) the relative error for each sample becomes smaller when the number of the reduced basis functions Nf and Nm increase. (2) POD-greedy RBM gives the robust and accurate approximation with the number of reduced basis functions Nf ≥ 4 and Nm ≥ 4. pf

2

pm

2

10

10

2 global basis functions 4 global basis functions 8 global basis functions

0

10

2 global basis functions 4 global basis functions 8 global basis functions

0

10

−2 −2

Relative L error

−4

10

2

Relative L2 error

10

10

−6

10

−4

10

−6

10

−8

10

−8

10 −10

10

−10

0

20

40

60

80

10

100

Index of the random samples

pf

2

20

40

60

80

100

pm

2

10

10

2 global basis functions 4 global basis functions 8 global basis functions

0

10

2 global basis functions 4 global basis functions 8 global basis functions

0

10

−2

−2

10

Relative L error

10

−4

10

−4

10

2

Relative L2 error

0

Index of the random samples

−6

10

−8

−8

10

10

−10

−10

10

10

−12

10

−6

10

−12

0

20

40

60

80

100

Index of the random samples

10

0

20

40

60

80

100

Index of the random samples

Figure 7.6: The relative error for 100 samples versus the number of the reduced basis functions being Nf = 2, 4 and 8, Nm = 2, 4 and 8 at time levels t = 250 (First line) and t = 500(Second line). Table 1 lists the average online CPU and the average relative error at different time levels for the two approaches: FEM and the POD-greedy reduced basis method. We find that: (1) the POD-greedy RBM significantly improves computation efficiency. (2) from the forth and fifth rows 19

Table 1: Comparison of the average online CPU time and the average relative error for the solution at different time levels for the two approaches: FEM and the POD-greedy reduced basis method (RBM), the number of the reduced basis functions is Nf = 8 and Nm = 8. CPU time Relative error

Strategies FEM RBM RBM εf RBM εm

t = 500 76.05 s 0.0682 s 3.91 × 10−10 5.28 × 10−10

t = 1000 149.37 s 0.1227 s 1.29 × 10−11 1.99 × 10−11

t = 1500 221.72 s 0.1587 s 1.30 × 10−11 2.01 × 10−11

of the table, we can see that the POD-greedy RBM gives a good accuracy with the number of the reduced basis functions being Nf = 8 and Nm = 8.

7.2

Numerical results for local-global reduced basis method

We consider the dual-continuum stochastic model in highly heterogeneous porous media in this numerical example. Here consider kf depicted in Figure 7.4, and km = kf /1000. The coupled GMsFEM is performed on a 10 × 10 coarse grid. we use L to denote the number of basis functions per coarse element. Let DOF be the number of degrees of freedom. In Table 2, we present the Table 2: The average relative errors (%) of the coupled GMsFEM at different time levels. (10 × 10 coarse grid) L

DOF

L=2 L=5 L=8

484 1210 1936

L=2 L=5 L=8

484 1210 1936

εf t = 100 0.74 8.10 × 10−2 2.94 × 10−2 t = 1000 0.96 7.71 × 10−2 2.81 × 10−2

εm 0.76 9.73 × 10−2 4.12 × 10−2 1.02 9.32 × 10−2 3.91 × 10−2

average relative errors for the coupled GMsFEM. From this table, we observe that when choosing L = 5 basis functions per coarse node, we can obtain an excellent result using the coupled GMsFEM. To show the influence of the proportion between coarse grid and fine grid to the coupled GMsFEM, we use 20 × 20 coarse grid and 10 × 10 coarse grid, and list the average relative errors for the coupled GMsFEM in Table 3. For this example, we still use 10 × 10 coarse grid, and choose L = 5 Table 3: The average relative errors (%) of the coupled GMsFEM at different time levels. (20 × 20 coarse grid) L

DOF

L=2 L=5 L=8

1764 4410 7056

εf t = 100 0.15 1.25 × 10−2 2.51 × 10−3

εm 0.17 1.86 × 10−2 5.73 × 10−3

basis functions per coarse node, and use the coupled GMsFEM to solve the equation (2.2). To get local-global reduced basis model for this dual-continuum stochastic model, we randomly choose nt = |Ξtrain | = 200 samples for localized uncertainties, take Ntf = 4, Ntm = 4, Mtf = 8 and 20

Mtm = 8, and use the Algorithm 2 to construct the local global reduced basis (RB) approximation spaces. We randomly choose 10000 samples and compute the average relative errors defined as (7.28), (7.29). Figure 7.7 shows the average relative error for the solutions versus different time levels with the number of the reduced basis functions being Nf = 2, 4 and 8, Nm = 2, 4 and 8. By the figure, we find that: (1) the average relative error becomes smaller when the number of local global reduced basis functions Nf and Nm increase. (2) From the pictures at the second line, we can see when choosing Nf = 8 and Nm = 8 local global reduced basis functions, the average relative errors for local global RBM is almost identical to the errors for the coupled GMsFEM at all simulation time levels. p

p

f

2

m

2

10

10

2 local global basis functions 4 local global basis functions 8 local global basis functions

1

10

2 local global basis functions 4 local global basis functions 8 local global basis functions

1

10

0 0

10

−1

Relative L error

10

−2

10

−1

10

2

Relative L2 error

10

−3

10

−4

10

−2

10

−3

10

−4

−5

10

−6

10

10

−5

10

−7

10

−6

0

500

1000

10

1500

0

500

Time

pf

2

1500

pm

2

10

10

2 local global basis functions 4 local global basis functions 8 local global basis functions coupled GMsFEM with L=5

1

10

2 local global basis functions 4 local global basis functions 8 local global basis functions coupled GMsFEM with L=5

1

10

0

Relative L error

10

0

10

2

Relative L2 error

1000

Time

−1

10

−2

10

−3

−2

10

−3

10

10

−4

10

−1

10

−4

0

500

1000

1500

Time

10

0

500

1000

1500

Time

Figure 7.7: The average relative error versus different time levels with the number of the reduced basis functions being Nf = 2, 4 and 8, Nm = 2, 4 and 8. The reference solution of the first (second) line is solved by FEM (coupled GMsFEM).

21

To visualize the individual errors of 100 samples, we plot the relative errors for the number of the reduced basis functions being Nf = 2, 4 and 8, Nm = 2, 4 and 8 in Figure 7.8, which shows that: (1) the relative error for each sample becomes smaller when the number of the reduced basis functions Nf and Nm increase. (2) POD-greedy RBM gives the robust and accurate approximation with the number of local global reduced basis functions Nf ≥ 4 and Nm ≥ 4. p

p

f

2

10

2 local global basis functions 4 local global basis functions 8 local global basis functions

1

10

0

0

Relative L2 error

10

−1

10

2

Relative L error

2 local global basis functions 4 local global basis functions 8 local global basis functions

1

10

10

−2

10

−3

10

−4

−1

10

−2

10

−3

10

−4

10

10

−5

−5

10

10

−6

10

m

2

10

−6

0

20

40

60

80

10

100

0

Index of the random samples

pf

1

20

60

80

100

pm

1

10

10

2 local global basis functions 4 local global basis functions 8 local global basis functions

2 local global basis functions 4 local global basis functions 8 local global basis functions 0

0

10

2

2

Relative L error

10

Relative L error

40

Index of the random samples

−1

10

−2

−2

10

10

−3

10

−1

10

−3

0

20

40

60

80

100

Index of the random samples

10

0

20

40

60

80

100

Index of the random samples

Figure 7.8: The relative error for 100 samples versus the number of the reduced basis functions being Nf = 2, 4 and 8, Nm = 2, 5 and 8 at time level t = 250. The reference solution of the first (second) line is solved by FEM (coupled GMsFEM). Table 4 lists the average online CPU and the average relative error at different time levels for the three approaches: FEM, coupled GMsFEM and local global RBM. We find that: (1) the approximation by the two approaches, i.e., the coupled GMsFEM and local global RBM is almost identical, which means that the coupled GMsFEM model can be accurately expressed by the local global RBM model. (2) the average CPU time per sample by local global RBM is the smallest 22

Table 4: Comparison of the average online CPU time and the average relative error for the solution at different time levels for the three approaches: FEM, coupled GMsFEM and the local global reduced basis method (RBM), the number of the local global reduced basis functions is Nf = 8 and Nm = 8. CPU time

Relative error

Strategies FEM coupled GMsFEM LG-RBM coupled GMsFEM εf coupled GMsFEM εm LG-RBM εf LG-RBM εm

t = 500 78.24 s 4.4238 s 0.0638 s 1.10 × 10−3 1.30 × 10−3 1.10 × 10−3 1.31 × 10−3

t = 1000 146.26 s 12.4505 s 0.1208 s 1.10 × 10−3 1.31 × 10−3 1.10 × 10−3 1.31 × 10−3

t = 1500 223.83 s 19.888 s 0.1860 s 1.10 × 10−3 1.31 × 10−3 1.10 × 10−3 1.31 × 10−3

among all of the approaches. The numerical example shows that the local global RBM achieves a good trade-off in both approximation accuracy and computation efficiency. Table 5 lists the average Table 5: Comparison of the average online CPU time and the average relative error for the velocity at different time levels for the three approaches: FEM, coupled GMsFEM and the local global reduced basis method (RBM). The number of the basis functions per coarse node for the coupled GMsFEM is L = 10, and the number of the local global reduced basis functions is Nf = 8 and Nm = 8. CPU time

Relative error

Strategies FEM coupled GMsFEM LG-RBM coupled GMsFEM εf coupled GMsFEM εm LG-RBM εf LG-RBM εm

t = 500 78.24 s 13.271438 s 0.0669 s 2.66 × 10−2 2.68 × 10−2 2.66 × 10−2 2.68 × 10−2

t = 1000 146.26 s 35.3515 s 0.1305 s 2.94 × 10−2 2.96 × 10−2 2.94 × 10−2 2.96 × 10−2

t = 1500 223.83 s 53.6640 s 0.1981 s 2.95 × 10−2 2.96 × 10−2 2.95 × 10−2 2.96 × 10−2

online CPU time and the average relative error for the velocity at different time levels for the three approaches: FEM, coupled GMsFEM and local global RBM. We note that the number of the basis functions per coarse node for the coupled GMsFEM is L = 10.

8

Conclusions

In the paper, we have presented reduced basis methods for solving dual-continuum stochastic models. It will be a great challenge when we solve the equations directly on the fine grid for a many-query situation. To significantly improve computation efficiency, we apply the reduced basis methods to construct a reduced order model for the dual-continuum stochastic model. To this end, we choose a few snapshots scattered in the manifold of the solutions. Then the proper orthogonal decomposition and greedy algorithm are used to select the reduced basis functions and get a reduced order model. For the reduced basis methods, an offline-online computational decomposition is achieved in the reduced basis methods to significantly improve computation efficiency. Although the offline computation may be expensive, the online computation is very efficient, which is very desirable for predicting the model’s outputs for various stochastic influence. If dual-continuum stochastic models have multiscale structures in inputs or is in heterogeneous random porous media, the offline computation will be very expensive. Because we have to use a very fine mesh to resolve 23

the features in all scales. To overcome the difficulty and improve the offline computation, we use the coupled Generalized Multiscale Finite Element Method (GMsFEM) to computer the snapshots in the offline stage of the POD-greedy RB method, then construct the reduced basis functions based on the coupled GMsFEM. The combination of the coupled GMsFEM and the POD-greedy RB method is called local global RB method in this work. The numerical results showed that the POD-greedy RB method performs well in both approximation accuracy and computation efficiency; Moreover, from the numerical results we can find that local global RB method has good performance in both approximation accuracy and computation efficiency for the dual-continuum stochastic models in heterogeneous random porous media.

Acknowledgments MV’s work is supported by the mega-grant of the Russian Federation Government (N 14.Y26.31.0013).

References [1] IY. Akkutlu, Y. Efendiev, M. Vasilyeva and Y. Wang, Multiscale model reduction for shale gas transport in a coupled discrete fracture and dual-continuum porous media, Journal of Natural Gas Science and Engineering, 2017. [2] P. Ackerer, N. Trottier and F. Delay, Flow in double-porosity aquifers: Parameter estimation using an adaptive multiscale method, Advances in Water Resources, 73 (2014), pp. 108–122. [3] J. Bear, Dynamics of fluids in porous media, Courier Corporation, 2013. [4] GI. Barenblatt, IP. Zheltov, and IN. Kochina, Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata], Journal of applied mathematics and mechanics, 24 (1960), pp. 1286–1303. [5] A. Buffa, Y. Maday, A.T. Patera, C. Prud’homme, and G. Turinici, A priori convergence of the greedy algorithm for the parametrized reduced basis method, ESAIM: Mathematical Modelling and Numerical Analysis, 46 (2012), pp. 595-603. [6] MBGM. Biermans, KM. Dijkema, and DA. De Vries, Water movement in porous media towards an ice front, Journal of Hydrology, 37 (1978), pp. 137–148. [7] Y. Chen, J.S.Hesthaven, Y. Maday, and J.Rodriguez, Certified reduced basis methods and output bounds for the harmonic Maxwells equations, SIAM Journal on Scientific Computing, 32 (2010), pp. 970-996. [8] E. T. Chung and Y. Efendiev, Adaptive multiscale model reduction with generalized multiscale finite element methods, Journal of Computational Physics, 320 (2016), pp. 69C95. [9] E. T. Chung Y. Efendiev W. T. Leung, Y Wang and M Vasilyeva, Non-local Multicontinua Upscaling for Flows in Heterogeneous Fractured Media 2017. [10] E. T. Chung, Y. Efendiev, W.T. Leung, and M. Vasilyeva, Coupling of multiscale and multi-continuum approaches, GEM-International Journal on Geomathematics, 8 (2017), pp. 9–41. 24

[11] V. Calo, Y. Efendiev, J. Galvis and G. Li, Randomized oversampling for generalized multiscale finite element methods, Multiscale Modeling and Simulation, 14 (2016), pp. 482– 501. [12] Z. Chen, G. Huan, and Y. Ma, Computational methods for multiphase flows in porous media, SIAM, 2006. [13] M.Chevreuil, A. Nouy, and E.Safatly, A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties, Comput. Methods Appl. Mech. Engrg. 255(2013), pp. 255C274 [14] Y. Efendiev, J. Galvis and E. Gildin, Local-global multiscale model reduction for flows in highly heterogeneous media, Journal of Computational Physivs, 231 (2012), pp.8100C8113. [15] Y. Efendiev, G. Galvis, and TY. Hou, Generalized multiscale finite element methods (GMsFEM),Journal of Computational Physics, 251 (2013), pp. 116–135. [16] Y. Efendiev, J. Galvis and T. Hou, Generalized multiscale finite element methods. oversampling strategies, International Journal for Multiscale Computational Engineering, accepted, 2013. [17] Y. Efendiev, T. Hou, and V. Ginting, Multiscale finite element methods for nonlinear problems and their applications, Comm. Math. Sci., 2 (2004), pp. 553C589. [18] P. Eisenklam, Flow Through Porous Media, Nature, 180(1957), pp. 1008–1009. [19] HC. Elman, and V. Forstall, Numerical solution of the parameterized steady-state Navier– Stokes equations using empirical interpolation methods, Computer Methods in Applied Mechanics and Engineering, 317 (2017), pp. 380–399. [20] Q. Guan, M. Gunzburger, CG. Webster, and G. Zhang, Reduced basis methods for nonlocal diffusion problems with random input data, Computer Methods in Applied Mechanics and Engineering, 317 (2017), pp. 746–770. [21] M.A. Grepl and M. Karcher, Reduced basis a posteriori error bounds for parametrized linear quadratic elliptic optimal control problems, Comptes Rendus Mathematique, 349 (2011), pp. 873-877. [22] C. Harris, Sedimentation rates, fluidization and flow through porous media, Nature, 184(1959), pp. 716. [23] L. Howle, R. Behringer, and J. Georgiadis, Visualization of convective fluid flow in a porous medium, Nature, 362(1993), pp. 230–232. [24] R. Handler, K. Housiadas, and A. Beris, Karhunen–loeve representations of turbulent channel flows using the method of snapshots, International journal for numerical methods in fluids, 52 (2006), pp. 1339–1360. [25] D.B.P. Huynh, G. Rozza, S. Sen, and A.T. Patera, A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants, Comptes Rendus Mathematique, 345 (2007), pp. 473-478. [26] T. Hayat, Z. Hussain, A. Alsaedi, and B. Ahmad, Numerical study for slip flow of carbon-water nanofluids, Computer Methods in Applied Mechanics and Engineering, 2017. 25

[27] J. Han, M. Wang, P. Bai, FR. Fikile, and MZ. Bazant, Dendrite suppression by shock electrodeposition in charged porous media, Scientific Reports, 6 (2016). [28] U. Hornung, Homogenization and porous media, Springer Science & Business Media, 6 (2012). [29] TJR. Hughes, The finite element method: linear static and dynamic finite element analysis, Courier Corporation, 2012. [30] L. Jiang and Q. Li, Reduced multiscale finite element basis methods for elliptic pdes with parameterized inputs, Journal of Computational and Applied Mathematics, 301 (2016), pp. 101–120. [31] L. Jiang and Q. Li, Model sparse representation based on reduced mixed GMsFE basis methods, Journal of Computational Physics, 338 (2017), pp. 285-312. [32] KT. Lim, and K. Aziz, Matrix-fracture transfer shape factors for dual-porosity simulators, Journal of Petroleum Science and Engineering, 13 (1995), pp. 169–178. [33] AF. Moench, Double-porosity models for a fissured groundwater reservoir with fracture skin, Water Resources Research, 20 (1984), pp. 831–846. [34] D. Modesto, S. Zlotnik, and A. Huerta, Proper generalized decomposition for parameterized helmholtz problems in heterogeneous and unbounded domains: application to harbor agitation, Computer Methods in Applied Mechanics and Engineering, 295 (2015), pp. 127–149. [35] A. Nouy, and F. Pled, A multiscale method for semi-linear elliptic equations with localized uncertainties and non-linearities, Mathematical Modelling and Numerical Analysis, 2017. [36] N. C. Nguyen, A multiscale reduced-basis method for parametrized elliptic partial differential equations with multiple scales, J. Comput. Phys. 227 (2008), pp. 9807-9822. [37] C. Prud’homme, D. Rovas, K. Veroy, Y. Maday, A. Patera, and G. Turinici, Reliable real-time solution of parametrized partial differential equations: reduced-basis output bounds methods, Journal of Fluids Engineering, 124 (2002), pp. 70-80. [38] DM. Pedroso, A solution to transient seepage in unsaturated porous media, Computer Methods in Applied Mechanics and Engineering, 285 (2015), pp. 791–816. [39] P. Pacciarini, and G. Rozza, Stabilized reduced basis method for parametrized advection– diffusion PDEs, Computer Methods in Applied Mechanics and Engineering, 274 (2014), pp. 1–18. [40] A. Quarteroni, G. Rozza, and A. Manzoni, Certified reduced basis approximation for parametrized partial differential equations and applications, Journal of Mathematics in Industry, 1 (2011), pp. 1–49. [41] G. Rozza, D.B.P. Huynh, and A.T. Patera, Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations, Arch Comput Methods Eng. 15 (2008), pp. 229-275. [42] E. Ranjbar, H. Hassanzadeh and Z. Chen, Effect of fracture pressure depletion regimes on the dual-porosity shape factor for flow of compressible fluids in fractured porous media, Advances in Water Resources, 34 (2011), pp. 1681–1693. 26

[43] BR. Scanlon, RE. Mace, ME. Barrett, and B. Smith, Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer, USA, Journal of Hydrology, 276 (2003), pp. 137–158. [44] M. Tene, MS. Al. Kobaisi and H Hajibeygi, Algebraic multiscale method for flow in heterogeneous porous media with embedded discrete fractures (F-AMS), Journal of Computational Physics, 321 (2016), pp. 819-845. [45] M. Tene, SBM. Bosma, MS. Al. Kobaisi and H. Hajibeygi, Projection-based Embedded Discrete Fracture Model (pEDFM), Advances in Water Resources, 105 (2017), pp. 205-216. [46] HK. Versteeg, and W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method, Pearson Education, 2007. [47] S. Whitaker, Flow in porous media I: A theoretical derivation of Darcy’s law, Transport in porous media, Springer, 1(1986), pp. 3–25. [48] YS. Wu, and K. Pruess, A multiple-porosity method for simulation of naturally fractured petroleum reservoirs, SPE Reservoir Engineering, 3 (1988), pp. 327–336. [49] D. Xiao, F. Fang, CC. Pain, and IM. Navon, A parameterized non-intrusive reduced order model and error analysis for general time-dependent nonlinear partial differential equations and its applications, Computer Methods in Applied Mechanics and Engineering, 317 (2017), pp. 868–889. [50] B. Yan, M. Alfi, C. An, Y. Cao, Y. Wang, and JE. Killough, General Multi-Porosity simulation for fractured reservoir modeling, Journal of Natural Gas Science and Engineering, 33 (2016), pp. 777–791.

27