The multiscale perturbation method for second order elliptic equations

ARTICLE IN PRESS

JID: AMC

[m3Gsc;January 24, 2020;10:27]

Applied Mathematics and Computation xxx (xxxx) xxx

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

The multiscale perturbation method for second order elliptic equations Alsadig Ali a, Het Mankad a, Felipe Pereira a,∗, Fabrício S. Sousa b a

Department of Mathematical Sciences, The University of Texas at Dallas, 800 W. Campbell Road, Richardson, Texas 75080-3021, USA Instituto de Ciências Matemáticas e de Computação, Universidade de São Paulo, Av. Trabalhador São-carlense, 400, São Carlos 13566-590, SP, Brazil

b

a r t i c l e

i n f o

Article history: Available online xxx Keywords: Porous media Domain decomposition Multiscale basis functions Robin boundary conditions Multiphase flows

a b s t r a c t In the numerical solution of elliptic equations, multiscale methods typically involve two steps: the solution of families of local solutions or multiscale basis functions (an embarrassingly parallel task) associated with subdomains of a domain decomposition of the original domain, followed by the solution of a global problem. In the solution of multiphase flow problems approximated by an operator splitting method one has to solve an elliptic equation every time step of a simulation, that would require that all multiscale basis functions be recomputed. In this work, we focus on the development of a novel method that replaces a full update of local solutions by reusing multiscale basis functions that are computed at an earlier time of a simulation. The procedure is based on classical perturbation theory. It can take advantage of both an offline stage (where multiscale basis functions are computed at the initial time of a simulation) as well as of a good initial guess for velocity and pressure. The formulation of the method is carefully explained and several numerical studies are presented and discussed. They provide an indication that the proposed procedure can be of value in speeding-up the solution of multiphase flow problems by multiscale methods. © 2020 Elsevier Inc. All rights reserved.

1. Introduction Problems of interest to the oil industry, such as the petroleum reservoirs in Brazil’s pre-salt layer, may require spatial discretizations involving billions of cells such that numerical solutions of multiphase flows can accurately reflect the heterogeneous nature of the underlying subsurface formations. Frequently used reservoir simulators cannot handle such large computational tasks. Multiscale methods are a class of numerical procedures that have good potential to be used in producing approximations for such challenging problems. These methods have been developed both in the context of finite volume and finite element approximations [1,2]. See [3,4] for references for more recent work. In the context of mixed finite elements we mention [5–7] that have been generalized by the recently introduced Multiscale Robin Coupled Method [3,4]. The aforementioned multiscale mixed methods share a similar algorithmic structure. Initially a set of local solutions (or mixed multiscale basis functions, MMBF) of boundary value problems (BVPs, the methods differ primarily in the choice of boundary conditions for these local problems) is solved for non-overlapping subdomains of the domain of interest. Then, a ∗

Corresponding author. E-mail addresses: [email protected] [email protected] (F.S. Sousa).

(A.

Ali),

[email protected]

(H.

Mankad),

[email protected]

(F.

Pereira),

https://doi.org/10.1016/j.amc.2019.125023 0 096-30 03/© 2020 Elsevier Inc. All rights reserved.

Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC 2

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

Fig. 1. The operator splitting technique. Note that the procedure allows for the use of distinct time steps: δ t for velocity and pressure and δ th (typically smaller because of a CFL restriction) for the water saturation.

global problem is defined, that will couple the multiscale basis functions to produce the final approximate solution for the global problem at hand. We are concerned with the development of the Multiscale Perturbation Method (MPM), a novel procedure to speed-up the parallel (multi-core) numerical solution of multiphase flow problems by multiscale methods that share the structure just described. In order to motivate the discussion of the proposed procedure for second order elliptic equations, we consider the numerical approximation of two-phase (water-oil) flows in porous media by an operator splitting method. Operator splitting methods [8–11] are well known strategies to handle the coupled systems of partial differential equation by decomposing them into smaller problems. The problems that result from this decomposition describe simpler physics and numerically they are solved sequentially. The governing system of equations for two phase flows is given by (see [9,12])

−λ(sw )κ∇ p = u,

∇ · u = q, φ

(1) (2)





∂ sw λw (sw ) +∇ · u = 0, ∂t λ ( sw )

(3)

where the unknowns are the velocity u(x, t) and pressure p(x, t) fields, and the water saturation sw (x, t), as a function of space x ∈  ⊂ Rd , d = 2, 3 and time t ≥ 0. For simplicity in the presentation additional assumptions considered in this model are the absence of capillary pressure and gravity effects, along with a fully saturated media; (sw + so = 1). The procedure introduced here could incorporate both gravity effects and capillary pressure, in that one would still need to solve a second order elliptic equation for the global pressure (see [9,12]). The assumption of a fully saturated media is typical of oil-water flows in petroleum reservoirs. In the above equations, q = q(x, t ) is a volumetric source term for wells, φ = φ (x ) is the porosity, κ(x) is the permeability tensor (we assume it is a scalar). The coefficient λ(sw ) is the total phase mobility, given by

λ ( sw ) =

κrw (sw ) κro (sw ) + , μw μo

(4)

where κrα (sw ), (α = w, o) are the relative permeabilities and μα is the viscosity of the α phase. We define the α phase mobility by

λα (sw ) =

κrα (sw ) , μα

(5)

such that λ(sw ) = λo (sw ) + λw (sw ). We define K (x, sw ) = λ(sw )κ (x ). Initial conditions for one of the saturations and boundary conditions have to be provided. For computational efficiency an operator splitting method is frequently used for the numerical approximation of the above system. This splitting is performed by “freezing” the saturation in time for the solution of the second order elliptic equation, allowing for the segregated solution of both equations. In one time step of the overall method, say from time tn to time tn+1 , saturation is kept at time tn for the solution of Eqs. (1) and (2), to obtain pn+1 and un+1 both at time tn+1 . These newly computed pressure and velocity are then used in Eq. (3) to evolve saturation to the next time step tn+1 . The main advantage is that now each equation can be solved by well known techniques that can benefit from different time steps. A larger time step is used for the update of velocity and pressure and smaller time steps, that satisfy a CFL condition [13], are used in the solution of a hyperbolic conservation law for the water saturation. The use of distinct time steps in illustrated in Fig. 1. The solution of two phase flows can naturally take advantage of multiscale methods that have been developed for second order elliptic equations by replacing a global solver by a multiscale one in the update of velocity and pressure. The proposed procedure aims at finding the solution for velocity and pressure by replacing the construction of a set of multiscale basis Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

3

functions or local solutions by reusing the multiscale basis functions computed for another elliptic equation. Thus, the proposed method is formulated in terms of two elliptic boundary value problems (let us refer to them as the initial time and target problems) that are approximated in terms of a common decomposition of the domain. For the initial time problem a set of multiscale basis functions is computed (say, in the first elliptic solution of a time dependent problem). Within the MPM the solution of the target elliptic equation is obtained in two steps. In the first step an interface problem is solved: perturbation theory is applied to predict the normal component of the fluxes on the skeleton of the domain decomposition. Through perturbation theory we derive a hierarchy of BVPs that share the same coefficients with the time zero problem. Thus, these BVPs can be solved by taking advantage of the pre-computed, time zero MMBFs. The first step is followed by a task that fits well in multi-core devices: for each subdomain a BVP is solved using the boundary conditions computed in the first step, giving the final approximation for the velocity field. The proposed method has some distinctive properties that can be summarized as follows: • The method works quite well for high-contrast formations; • A considerable fraction of the computational work can be performed in an embarrassingly parallel offline stage (all solutions of local problems can be naturally solved independently from each other); • If a full set of multiscale basis functions (one multiscale basis function for each edge of the boundary of a subdomain) is computed for the time zero problem the target solution has continuous fluxes (thus, a downscaling step can be avoided); • The method can incorporate a good initial guess (that would be relevant in the solution of time dependent problems); • The cost of the procedure is considerably smaller than a typical multiscale method; • The procedure in non-intrusive: it can incorporate various multiscale methods in their original formulation and numerical implementation; • Numerical studies indicate that the accuracy of the procedure is comparable (or better) to that of the multiscale method used in its formulation. This work is organized as follows. In Section 2 we recall the Multiscale Mixed Method [7] that will be the underlying multiscale method for the description of the proposed procedure. Then, the new method is derived in detail in Section 3. Initially we discuss a mixed (for the flux and pressure variables) perturbation method for a BVP. Then we derive an auxiliary residual problem that allows one to produce improved solutions by taking advantage of an initial guess (that is frequently available in the solution of time dependent methods). Section 4 is dedicated to a discussion of numerical solutions of a model problem. Our conclusions and the prospects of future work appear in Section 5. The dimensionless form of the Poisson problem appears in an appendix. 2. A multiscale mixed method for the elliptic equation As pointed out earlier in the discussion of the operator splitting method the second order elliptic equation and the hyperbolic equation are solved sequentially. In this work we are concerned with the solution of the elliptic problem that is usually the most expensive one, and how we can accelerate it within the operator splitting framework. Let us consider the following model elliptic problem

u = −K ∇ p x ∈ ,

(6)

∇ · u = q x ∈ ,

(7)

which holds in  ⊂ d = 2, 3, a bounded domain with Lipschitz boundary, and with suitable boundary conditions p = g p ˇ = gu in x ∈ ∂ u (n ˇ is the exterior unit normal vector), with ∂  = ∂  p ∪ ∂ u and ∂  p ∩ ∂ u = ∅. for x ∈ ∂ p and u · n The Multiscale Perturbation Method that we propose will be described using the Multiscale Mixed Method (MuMM) [7] to approximate the elliptic problems (6) and (7). The MuMM is a non-overlapping multiscale domain decomposition method that greatly accelerates the parallel solution of the elliptic equation described in (6) and (7). Here flux and pressure conservation is ensured by enforcing suitable Robin type boundary conditions between subdomains. Let Rd ,

H (div ; ) = {v ∈ (L2 ())d : ∇ · v ∈ L2 ()}, ˇ = g on Vg () = {v ∈ H (div ; ) : v · n

∂ u } ,

be suitable Sobolev spaces, and inner products (·, · ) in L2 () and ( · , · )∂  in L2 (∂ ). Consider i , i = 1, . . . , M, a subdivision of the domain  such that  = ∪M i and let be the skeleton of this decomposition, i.e., the union of all i i j = ji = ∂ i ∩ ∂  j . Moreover, let us define ip = ∂  p ∩ ∂ i as well as ui = ∂ u ∩ ∂ i to simplify notation. In the original MuMM, local solutions are solved for each i , i = 1, . . . , M, with Robin boundary conditions enforced on , that are

ˇ i + pi = β j u j · n ˇ j + pj, −βi ui · n for each ij and ji , that are equivalent to the imposition of the continuity conditions

pi | i j = p j | i j

and

ˇi + uj · n ˇ j = 0, ui · n

on i j .

The variational formulation of these local problems reads: Find { pi , ui } ∈ L2 (i ) × Vgu (i ) such that Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC 4

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

(∇ · ui , w )i = (q, w )i ,

(8)

(K −1 ui , v )i = ( pi , ∇ · v )i − (g p , v · nˇ i ) ip − (βi ui · nˇ i + r, v · nˇ i )∂ i \∂  ,

(9)

for all w ∈ is

L2 (i )

and v ∈ V0 (i ). Note that in the above formulation the general Robin boundary condition that is imposed

ˇ i + pi = r, −βi ui · n assuming that r is given. The MuMM method is derived from this variational formulation by a suitable choice of approximation spaces. Let Th be a partition of  into d-dimensional rectangles of size h, which is referred to as the fine mesh. Without loss of generality, let the subdomains i , i = 1, . . . , M be defined in terms of this fine mesh, as d-dimensional rectangles of size H h such that H/h is an integer number. The authors in [7] introduce a third length scale H bounded below (above) by h (H), such that h ≤ H ≤ H and H/H is also an integer number. While the domain decomposition is performed at the H scale, the enforcing of conservation through the Robin boundary conditions between subdomains will be performed at the H scale. 

Note that ij can be further divided in H/H segments i j such that



H/H

i j =



i j .

=1

Consider the lowest order Raviart–Thomas finite element space (RT0 ), by selecting suitable finite dimensional spaces i ⊂ V ( ) and Q i ⊂ L2 ( ). The MuMM method is an iterative procedure (for α = 1, 2 . . . ), that for each subdomain i = Vh,g g i i h i 1, . . . , M seeks {uih , pih } ∈ Vh,g × Qhi , with u

uih

=

lim ui,(α ) α →∞ h

and

pih = lim pi,h (α ) , α →∞

such that

(∇ · ui,h (α ) , w )i = (q, w )i ,

(10) 

(K −1 ui,h (α ) , v )i = ( pi,h(α ) , ∇ · v )i − (g p , v · nˇ i ) ip −

 ij⊂

∂ i

(βi ui,h (α ) · nˇ i + ri,j ( ) , v · nˇ i )  ,

(11)

ij

hold for all w ∈ and v ∈ The Robin boundary condition values are taken from the nearest neighbors j (to i ), from a previous known iteration. In the MuMM [7] a red-black ordering is introduced by classifying subdomains as red or black in a checkerboard fashion, so that one can take advantage of the most recent data available. Hence Qhi

Vh,i 0 .



ri,j ( )

=

ri,j (α −1) if

∂ i is a red subdomain,

ri,j (α )

∂ i is a black subdomain,

where

ri,j (α ) =

1



H

if

  ij

β j uhj,(α ) · nˇ j + phj,(α ) dl,

(12)

(13)

is the L2 -projection of the fine scale Robin boundary condition in the neighbour j to the H scale. This is well defined, provided that all red domains are processed before the black ones. The iteration sequence is as follows: Eqs. (10) and (11) are first solved for all red subdomains, taking ri,j (α −1 ) from the previous iteration from the neighboring black subdomains, and then the same is performed for all black subdomains, taking the most recent ri,j (α ) from the red subdomains.

Remark 1. As in [7], Eqs. (10) and (11) are solved using a hybridized mixed finite element approximation, so that the j, (α ) pressure unknowns evaluated at the edges appearing in Eq. (13) are in fact Lagrange multipliers, say λh , which are computed by the hybridized method. Remark 2. The red-black iteration procedure is in fact a Gauss–Seidel-like iteration that is incorporated in the parallel procedure of the original MuMM method. 2.1. Multiscale basis functions As pointed out in [7], although the iterative method presented in the previous section converges to an approximate solution of the original elliptic equation, there is no need to fully solve local problems given by Eqs. (10) and (11) for every Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

5

iteration step. This can be avoided by precomputing a set of MMBFs ij = {ψij , φij } and i = {χi , ϕi } in such a way that the numerical solution is written as



{uih , pih } =

 ij⊂

Ai j ij + i .

(14)

∂ i

The basis functions ij = {ψij , φij } are computed through the solution of the uncoupled local problems

(∇ · ψij , w )i = 0,

(15)

(K −1 ψij , v )i = (φij , ∇ · v )i − (βi ψij · nˇ i + gi j , v · nˇ i )  ,

(16)

ij

with

gi j =



1 in i j , 0 in



∂ i \ i j ,

for each segment  of each ij associated with subdomain ∂ i , while the basis function i = {χi , ϕi } is the solution of

(∇ · χi , w )i = (q, w )i ,

(17)

(K −1 χi , v )i = (ϕi , ∇ · v )i − (βi χi · nˇ i , v · nˇ i )∂ i \∂  − (g p , v · nˇ i ) ip ,

(18)

for each domain i , to take into account the contribution of the non-homogeneous source term as well as non-homogeneous pressure boundary conditions. With these precomputed sets of MMBFs, the coefficients of the linear combination (14) are taken as Ai j = ri,j ( ) (defined in Eq. (12)) and the approximate solution of the MuMM iterative procedure can be simply computed as







ui,h (α ) , pi,h (α ) =

 ij⊂

ri,j ( ) ij + i ,

(19)

∂ i

in such a way that the linear combination (19) can be used in the red-black iterative procedure to converge to the actual solution of the MuMM method {uih , pih }. Note that there is no linear algebra problem to be solved for each iteration α , only the evaluation of the linear combination (19), which is much cheaper, and the main advantage of this method. The final algorithm is summarized in Algorithm 1. Algorithm 1 Multiscale Mixed Method. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

Precompute the MMBF ij and i from (15)-(18) while convergence is not achieved do for all red domains do Evaluate {ui,h (α ) , pi,h (α ) } on from (19) using ri,j ( ) = ri,j (α −1 ) end for for all black domains do Evaluate {ui,h (α ) , pi,h (α ) } on from (19) using ri,j ( ) = ri,j (α ) end for α ←α+1 end while Evaluate the converged overall solution {uih , pih } in 

Note that the computational cost of Algorithm 1 resides in its Step 1, where a set of local boundary value problems are solved in each subdomain. The iterative procedure (Steps 2 to Step 10 of Algorithm 1) refer to the solution of a coarse interface problem and can be solved very fast (see [14,15] for a very efficient algorithm for the parallel solution of this interface problem). Thus the computational cost of Algorithm 1 can be accurately estimated by the number of local solutions in each subdomain that are solved in its Step 1. Remark 3. By using the Lagrange multipliers from the hybridized finite element formulation, one can compute the MuMM solution without having to reconstruct the flux and pressure fields in . Since Lagrange multipliers are defined on the edges of Th , {ui,h (α ) , pi,h (α ) } are reconstructed only on the skeleton , reducing even further the cost of achieving convergence. The final fields can be reconstructed after convergence of the iterative process. Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC 6

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

3. The multiscale perturbation method In the numerical solution of two-phase flow problems by an operator splitting method, if the MuMM is used to approximate the elliptic problem (along with some other method to solve the hyperbolic equation), a performance issue arises. Since saturation is updated by the hyperbolic solver, the solution of the elliptic equation in the next time (when the algorithm calls for an update of velocity and pressure) will have a new set of coefficients K = λ(sw )κ . This requires the recalculation of all MMBFs in the MuMM. Since the update of the saturation is usually linked to a movement of an oil-water interface throughout the domain, this change in saturation is usually concentrated around the interface. An approach to accelerate the overall method, commonly seen in the literature [2], is to recompute only the multiscale basis function of the subdomains that are crossed by the oilwater interface. This procedure, however, would obviously generate load imbalance among processes of a parallel solver. The new method introduced in this work seeks to reuse the precomputed MMBFs from previous time steps, instead of actually recomputing them, thus reducing computation cost and avoiding load imbalance of the overall algorithm in the case of a parallel implementation. To introduce the ideas of our method let us consider only the elliptic equation that arises in the modelling of twophase flow in a domain  ⊂ Rd , d = 2, 3, written in dimensionless form (see A.8). Let us also introduce two stationary problems, one referred as Pt0 (the time zero problem with known coefficient K0 ) and another referred as Pt f (the target problem with known coefficient Kf ). One can think of them as elliptic problems associated with fictitious times t0 and tf > t0 respectively, representing two different steps of the overall time-marching scheme. Our goal is to find a solution for Pt f in terms of the MMBFs computed for the solution of problem Pt0 without having to recompute the whole set of multiscale basis functions. We will also assume that an approximation for the solution of the target problem is given and these pressure and velocity fields are denoted by (pref , uref ). In the solution of multiphase flow problems these fields are typically computed by extrapolation in time of the corresponding quantities computed at earlier time steps. The time zero and target problems are given by:

⎧ u0 = −K0 ∇ p0 ⎪ ⎪ ⎪ ⎨∇ · u = q 0 Pt0 ⎪ p0 = g p ⎪ ⎪ ⎩ u ˇ =g u0 · n

and

⎧ u = −K f ∇ p ⎪ ⎪ ⎪ ⎨∇ · u = q Pt f : ⎪ p = gp ⎪ ⎪ ⎩ u ˇ =g u·n

in , in ,

(20)

∂ p, on ∂ u , on

in , in ,

(21)

∂ p, on ∂ u . on

In these equations, q = q(x ) is considered to be a function of space. The Dirichlet (gp ) and the Neumann boundary functions (gu ) are known. In the proposed method we write the target pressure and velocity fields of Pt f as a perturbation of the approximations of the reference solutions (pref , uref ) as follows:

p = pre f + δ p,

(22)

u = ure f + δ u.

(23)

3.1. The problem for the perturbations (δ p, δ u) Using (22) and (23) in the target problem we get:

⎧ ure f + δ u − (−K f ∇ pre f ) = −K f ∇δ p ⎪ ⎪ ⎪ ⎨ ∇ · (ure f + δ u ) = q ⎪ pre f + δ p = g p ⎪ ⎪ ⎩ (ure f + δ u ) · nˇ = gu

in , in ,

∂ p, on ∂ u . on

(24)

Considering ua = ure f + δ u − (−K f ∇ pre f ) we can rewrite (24) as follows:

⎧ ⎪ ⎪ua = −K f ∇δ p ⎪ ⎨∇ · u = q − ∇ · (−(K )∇ p ) a f re f ⎪δ p = gp − pre f ⎪ ⎪ ⎩ ˇ = gu − (−(K f )∇ pre f ) · n ˇ ua · n

in , in ,

∂ p, on ∂ u . on

(25)

Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

7

We remark that (25) can be simplified if the reference solution satisfies ure f = (−K f ∇ pre f ). This reduces the problem for (δ u, δ p) to:

⎧ δ u = −K f ∇δ p ⎪ ⎪ ⎪ ⎨∇ · δ u = q − ∇ · u re f p ⎪ δ p = g − p re f ⎪ ⎪ ⎩ δ u · nˇ = gu − ure f · nˇ

in , in ,

(26)

∂ p, on ∂ u . on

Note that the solution of both problems for (δ p, δ u) defined above would obviously require the same amount of work as the solution of the target problem given by (21). We now develop an approximate solution for (δ p, δ u) that can take advantage of MMBFs computed for the time zero problem. 3.2. The interface problem Here we develop an approximation for the flux interface values (the flux variables on ) by using perturbation theory (see [16]) for (δ p, δ u). The approximate MPM solution of (21) will use these interface flux values in setting Neumann boundary conditions for auxiliary local BVPs defined in Section 3.3. In this development we view the target permeability Kf as a perturbation of the time zero permeability K0 . Hence, given the fields Kf and K0 we can express Kf as,



where

 (K f − K0 )  = K0 +  Kp , (K f − K0 )

K f = K0 + (K f − K0 ) 

(27)

   = (K f − K0 ),

(28)

and

Kp =

(K − K0 ) (K − K0 )  f = f . (K f − K0 ) 

(29)

Here  is a dimensionless parameter which gives a measure of how much Kf differs from K0 (due to changes in λ(sw )). Typically   1. If this is not the case it is an indication that a smaller time step should be considered for the update of the elliptic equation. L2 norm is used to compute  . Kp is an auxiliary field which assists along with  to measure the perturbation that has occurred from K0 to Kf . For high-contrast, realistic permeability fields our experiments show that  ≈ 10−2 . We disregard the term that is proportional to  2 ≈ 10−4 : our experiments show that good accuracy is obtained by this approximation in comparison with a direct MuMM solution. We remark that the relative error of the velocity field computed by state of the art multiscale methods (with respect to the direct fine grid solution) for high contrast media is typically of the order of 10−1 (see [3]). Hence, consider the following expansion,

δ u = δ u0 + δ u1 ,

(30)

δ p = δ p0 + δ p1 .

(31)

0:

The BVPs associated with Substituting (30) and (31) into (25) we get,

⎧ ua0 = −K0 ∇δ p0 ⎪ ⎪ ⎨∇ · u = q − ∇ · (−K ∇ p ) a0 f re f 0 : δ p = ⎪ 0 gp − pre f ⎪ ⎩ ˇ = gu − (−K f ∇ pre f ) · n ˇ ua0 · n

in , in ,

∂ p, on ∂ u . on

(32)

Here, ua0 = ure f + δ u0 − (−K f ∇ pre f ). If we substitute (30) and (31) into (26) then, the corresponding BVP associated with  0 is,

⎧ δ u0 = −K0 ∇δ p0 ⎪ ⎪ ⎨∇ · δ u = q − ∇ · u 0 re f 0 : p δ p = g − p ⎪ re f ⎪ ⎩ 0 δ u0 · nˇ = gu − ure f · nˇ

The BVP associated with  1 :

1 :

⎧ ˆ 1 = −K0 ∇δ p1 u ⎪ ⎪ ⎨∇ · uˆ = −∇ · u˜ 1

δp = 0 ⎪ ⎪ ⎩ 1

ˇ =0 uˆ1 · n

0

in , in ,

∂ p, on ∂ u . on

(33)

in , in ,

∂ p, on ∂ u , on

(34)

Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC 8

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

˜ 0 = −K p ∇δ p0 and δ u1 = u ˜0 + u ˆ 1 . Note that (34) remains the same for the general case which is (25) as well as the where u simplified case which is given in (26). The solution of (33) and (34) are approximated by the MuMM, taking advantage of MMBFs defined for the zero time problem. Hence, the final solution is given by the linear combination defined in (19) with a Robin interface condition and H > h and a new i (for each subdomain) computed for the problem at hand. Here, the MMBFs (ij ) are computed and stored when the time zero problem Pt0 is solved. They are reused to find the solution of the discrete system corresponding to the zeroth and the first order of  derived in (33) and (34). This is the main advantage of MPM in terms of reduced computational cost. 3.3. Local update of the solution The fluxes defined in (3.2) in the skeleton of the domain decomposition are set as boundary conditions for local BVPs. We will denote the final flux and pressure obtained from the perturbation step as (ump , pmp ). The BVP for each subdomain can be stated as:

⎧ i u = −K f ∇ pi ⎪ ⎪ ⎪ i ⎪ ⎨∇ · u = q

ˇ i = ump · n ˇi PN : ui · n

⎪ ⎪ ⎪ pi = g p ⎪ ⎩

ˇ = gu ui · n

in i , in i , on i j ,

(35)

∂ ip , on ∂ iu . on

Here, i , i = 1, . . . , M denotes the number of subdomains in the domain . By solving (35) we obtain an approximate solution to the problem Pt f . The parallel algorithm for the MPM is shown next. 3.4. The parallel algorithm for MPM The parallel algorithm for the MPM, consisting of two stages (offline and online) and two steps for the online stage (interface flux calculation followed by local updates), is described in Algorithm 2. Algorithm 2 The MPM scheme to solve Pt f in parallel. Offline Stage: Compute (in parallel) a set of MMBFs ij and 0i for each subdomain for the problem Pt0 (20) using (15 - 18). Online Stage: 2: Interface problem (The Perturbation Step):

1:

a. Solve (32) using the Algorithm (1) discussed in section (2.1) (Skip Step 1 of Algorithm (1)). b. Solve (34) using the Algorithm (1) discussed in section (2.1) (Skip Step 1 of Algorithm (1)). 3:

Local update (MPM Solution): a. Compute (in parallel) the local updates defined in (35).

We remark upon the fact that Algorithm 2 has effectively a fixed cost. First, note that the offline stage of the calculations should not be considered for the purpose of assessing the cost of Algorithm 2. Thus, its cost is given by the solution of three local solutions in each subdomain (one local solution computed in Steps 2.a, 2.b and 3.a). In line with the discussion above on the computational cost of Algorithm 1 this estimate does not take into account the (negligible) cost of the solution of the coarse interface problems associated with Steps 2.a and 2.b of Algorithm 2. The reduced cost of the MPM (Algorithm 2) with respect to the MuMM (Algorithm 1) can be clearly seen: while the MuMM requires the solution of a large number of multiscale basis functions in each subdomain, the MPM needs only 3 local solutions. We further remark that the solution computed by MPM can be further refined, for instance using the primal-dual smoothing of [15]. Note that if we consider H > h a downscaling step is required to define velocity fields with continuous normal components. 4. Numerical results In this section, we consider two different examples to discuss the accuracy obtained by our method. The first example is the quarter of a 5-spot problem (the 5-spot problem is a classical benchmark example in reservoir engineering [12]; the quarter of a 5-spot problem is defined on a square domain with injection and production wells placed at opposite corners of the domain). In the second example we consider a slab geometry, and we discuss the results obtained by using the permeability data of layer 36 of the SPE10 problem. The quarter of a 5-spot problem is discussed with three different Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

9

Fig. 2. Log permeability field for the quarter of a 5-spot problem. Left: The time zero permeability (K0 ). Center: The target permeability field (Kf ) and the dark lines indicate the position where the α parameter of (36) has a discontinuity. Right: The permeability field K p = (K f − K0 )/ showing the perturbation between the time zero and the target permeability fields. Here, K f,max /K f,min = 105 . Thus, Kf has a jump of 0.5 units in the [0, 0.6] × [0.4, 1] region of the domain with 40 × 40 fine elements in the top left corner, when compared to K0 . Note that the geometry of the jump mimics an invading water front.

reference solutions to investigate the importance of a good reference solution in order to obtain improved accuracy. The examples mentioned above are tested with different H scales and the accuracy of each of the results obtained is measured by comparing it to the fine grid solution. 4.1. The quarter of a 5-spot problem Here we indicate the accuracy of the proposed method by considering a classical 5-spot geometry with a considerable perturbation introduced in the field K0 that mimics an invading water front in a reservoir containing initially oil. The computational domain for this problem is  = [0, 1] × [0, 1] with 64 × 64 fine grid elements. The global boundary conditions are of Dirichlet type in the top left and the bottom right corners. The remaining global boundary conditions are no flow (Neumann) boundary conditions (see Fig. A.6). The time zero and the target permeability fields K0 (x) and Kf (x) respectively are computed as

α (x )eδξ (x) ,

(36)

where K0 (x) considers

α (x ) = 0.3, x ∈ [0, 1] × [0, 1], and Kf (x) considers



α (x ) =

0.8, x ∈ [0, 0.6] × [0.4, 1], 0.3, otherwise.

Here, α (x) is used to control the magnitude of the perturbation between the two permeability fields. δ is chosen in such a way that the permeability contrast K f,max /K f,min = 105 and 108 , the field ξ (x) is a self similar Gaussian distribution having zero mean and the covariance function given by C (x, y ) = |x − y|−1/2 . Note that both the fields are in a dimensionless form (see Appendix A). The maximum value of the time zero permeability field, K0,max is used to obtain the dimensionless form for the time zero and the target permeability fields. The log permeability fields for K0 , Kf and Kp are shown in Fig. 2. The domain is divided into 4 × 4 subdomains each containing 16 × 16 fine grid elements. This makes the coarse scale to be H = 1/4 and fine scale h = H/16 for each subdomain. The proposed method is tested on two intermediate scales namely, H = 2h = H/8 and H = 4h = H/4 along with the fine scale h. In addition, we have tested the method using three different reference solutions. The first reference solution is the fine grid solution obtained from the permeability data that has the same units of jump as the target field but the region in the domain that it covers is 39 × 39 region of fine cells in the top left corner of the domain unlike the target permeability which has a jump in 40 × 40 region of the domain. The second reference solution is the fine grid solution obtained from the permeability data that has a jump spanning over 36 × 36 region in the top left corner of the domain with the same units of jump as the target permeability and the third reference solution is obtained by solving the time zero problem (Pt0 ) at fine scale h. The velocity field corresponding to the solution obtained from the perturbation step and the MPM by considering the 39 × 39 reference computed at the fine scale h is shown in Fig. 3. We can observe from the tables below (see Tables 1, 2, and 3) that the MPM solution (for all reference solutions considered in our studies) is computationally less expensive and more accurate than the corresponding MuMM solution. Concerning the comparison of the computational costs of MuMM and MPM note that, for the problem at Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC 10

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

Fig. 3. Velocity field for the quarter of a 5-spot problem. Left: The solution of the target problem Pt f (green) superimposed to the solution of the Perturbation step (magenta) and the final MPM solution (blue). Here the permeability contrast is K f,max /K f,min = 105 and (uref , pref ) is the solution of the 39 × 39 reference permeability data. Right: The zoomed version of the vector field highlighted with a black box in the picture on the left. We can note from this figure that there is a very good agreement between the fine grid solution and the solution obtained from MPM. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 1 Relative Flux Error for the quarter of a 5-spot problem obtained using the 39 × 39 reference solution. Here, K f,max /K f,min = 105 and 108 . H

Perturbation MPM MuMM

K f,max /K f,min = 105

K f,max /K f,min = 108

h

2h

4h

h

2h

4h

8.8×10−3 4.8×10−3 -

9.0×10−3 5.5×10−3 5.6×10−2

1.0×10−2 7.0×10−3 9.4×10−2

1.3×10−2 9.6×10−3 -

1.2×10−2 9.6×10−3 8.6×10−2

1.2×10−2 1.1×10−2 1.4×10−1

Table 2 Relative Flux Error for the quarter of a 5-spot problem obtained using the 36 × 36 reference solution. Here, K f,max /K f,min = 105 and 108 . H

Perturbation MPM MuMM

K f,max /K f,min = 105

K f,max /K f,min = 108

h

2h

4h

h

2h

4h

2.6×10−2 1.8×10−2 -

3.0×10−2 2.3×10−2 5.6×10−2

4.1×10−2 3.8×10−2 9.4×10−2

3.4×10−2 2.5×10−2 -

4.3×10−2 3.5×10−2 8.6×10−2

6.1×10−2 5.5×10−2 1.4×10−1

Table 3 Relative Flux Error for the quarter of a 5-spot problem obtained using the time zero solution as the reference solution. Here, K f,max /K f,min = 105 and 108 . H

K f,max /K f,min = 105 h

Perturbation MPM MuMM

K f,max /K f,min = 108 2h

−2

4.6×10 4.5 ×10−2 -

4h −2

5.0×10 4.9×10−2 5.6×10−2

h −2

5.6×10 5.8×10−2 9.4×10−2

2h −2

2.5×10 2.4×10−2 -

4h −2

3.7×10 3.6 ×10−2 8.6×10−2

5.5×10−2 5.5×10−2 1.4×10−1

hand, if H = h, 2h, and 4h the numbers of multiscale basis functions needed for each subdomain for the MuMM solution are, respectively, 256, 64 and 16. In line with our remarks about Algorithm 2 the MPM cost is fixed (three multiscale basis functions for each subdomain). Note also that MPM with a full set of multiscale basis functions (H = h) computed at time zero (columns 1 and 4 of Tables 1, 2 and 3) produces more accurate velocity fields than the ones produced by the MuMM with H = 2h, 4h. The velocity field produced by MPM with H = h has continuous normal components, and a downscaling step would not be needed to couple this field with a transport equation. Thus, MPM can take full advantage of an offline stage of computation, where the multiscale basis functions are computed. Note also that, as expected, the error increases as the quality of the reference solution decreases. We consider the L2 norm in the relative errors reported here. The approximate solution takes into account the contributions from the  0 and  1 terms. Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

JID: AMC

ARTICLE IN PRESS

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

11

Fig. 4. Log permeability field for the SPE10 layer 36. Top: The time zero permeability (K0 ). Center: The target permeability field (Kf ). Bottom: The permeability field K p = (K f − K0 )/ showing the perturbation between the time zero and the target permeability fields. Here, K f,max /K f,min = 105 and Kf has a jump of 0.5 units in the region lying to the left of the sine wave.

4.2. SPE10 layer 36 In this example we test our method on the permeability field provided by the Society of Petroleum Engineers (SPE) for the SPE10 project (https://www.spe.org/web/csp/) [17]. This is considered to be a benchmark problem in this field. The original 3D dataset consists of 220 × 60 × 85 cells spanned through several layers in the z direction. We consider a 2D dataset, the layer 36 of the SPE10 project that contains a high contrast channel to test our method. Our domain for this problem is  = [0, 11/3] × [0, 1] which is discretized into 220 × 60 cells and 11 × 3 subdomains with each subdomain containing 20 × 20 fine cells. Here H = 1/3 and we use three different scales which includes the fine scale h = H/20 and two intermediate scales namely H = 2h = H/10 and H = 4h = H/5. The top and bottom global boundary conditions are no flow and the left and the right boundaries have Dirichlet boundary conditions. There are no sources and sinks in this example. The time zero permeability K0 (x) and the target permeability field Kf (x) are defined in the same manner as in (36) where K0 (x) considers

α (x ) = 0.3, x ∈ [0, 11/3] × [0, 1], Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC 12

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx Table 4 Relative Flux Error for the SPE10 layer 36 permeability field. Here K f,max /K f,min = 105 . H Perturbation MPM MuMM

h

2h −2

2.9×10 6.5×10−3 -

4h −2

2.9×10 7.0×10−3 1.2×10−1

3.0×10−2 1.2×10−2 2.6×10−1

Fig. 5. Velocity field for the SPE10 problem. Top: The solution of the target problem Pt f (green) superimposed to the solution of the Perturbation step (magenta) and the final MPM solution (blue). Here the permeability contrast is K f,max /K f,min = 105 . Bottom: The zoomed version of the vector field highlighted with a black box in the picture on the Top. Note the very good agreement between the MPM and fine grid solutions in the case of a challenging, channelized permeability field. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

and Kf (x) considers,

α (x ) =

0.8, x ∈ region lying to the left of the sine wave (see Fig. 4 ), 0.3, otherwise.

The perturbation jump is of the same height as the quarter of a 5-spot problem but the geometry is different (see Fig. 4). The reference solution used in this example follows the geometry similar to that of the target permeability field. The only difference is the area of the jump region is a little bit smaller than the one used in the target permeability data. Table 4 shows the relative errors for the case where K f,max /K f,min = 105 . Again, L2 norm is used to compute all the relative errors. The velocity field in Fig. 5 shows that there is a good level of agreement between the fine grid solution and the solution obtained using MPM. We can conclude that the MPM can deal with high-contrast data like the one used on this example, and give us a better approximation to the fine grid solution even when using a H scale larger than the fine scale h. Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

13

5. Conclusions and future work We have described a new method (MPM, the multiscale perturbation method) with great potential to reduce the computational cost of the solution of multiphase porous media flows when approximated by multiscale methods. The procedure is non-intrusive (in that it can use existing codes for multiscale methods without modifications), and has been designed for parallel processing. Numerical results indicate that the proposed method produces accurate approximations for second order elliptic equations with high-contrast coefficients. The investigation of a parallel implementation of the method proposed here, along with its use in the solution of two-phase flows approximated by an operator splitting technique, is left for future work. Besides the use of MPM to approximate subsurface flows by operator splitting techniques we envision that MPM can also be important in speeding up the implicit solution of multiphase flows. Another research area that can benefit from our developments is uncertainty quantification (UQ). In Markov chain Monte Carlo methods with a random walk sampler [18] one has to solve several thousands of numerical simulations with absolute permeabilities that are similar to each other (see, for instance [19]). Again, the reuse of multiscale basis functions can be of great benefit to reduce computational cost in such studies. The authors are currently investigating these topics. Acknowlgedgments F. S. Sousa gratefully acknowledges the financial support received from the Brazilian oil company Petrobras and from the São Paulo Research Foundation FAPESP, CEPID-CeMEAI grant 2013/07375-0. F. Pereira was funded in part by NSF-DMS 1514808, a Science Without Borders/CNPq-Brazil grant and UT Dallas. F. Pereira and F. S. Sousa would like to thank Dr. Michael Vynnycky for instructive discussions about perturbation theory. Appendix A. Dimensionless form of the Poisson equation In order to understand the contributions made by various physical quantities towards the occurrence of a physical phenomenon, we need to have a common point of reference for comparing these different physical quantities. Nondimensionalization of all the physical quantities involved helps in this case. Consider the following second order elliptic P.D.E. with non trivial Dirichlet boundary conditions as shown in Fig. (A.1),

⎧ u(x ) = −λ(sw )κ (x )∇ p(x ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎨∇ · u(x )= q(x ) pin on in , ⎪ p ( x ) = ⎪ ⎪ pout on out , ⎪ ⎪ ⎩ ˇ =0 u (x ) · n

in , in , (A.1) otherwise.

Our goal here is to derive a dimensionless form of the problem described in (A.1). Here the pressure p has units m2

N , m2

source term q(x) is 1s , flux u(x, t) is ms , the total phase mobility λ(sw ) has units Ns and finally the permeability κ (x) is measured in Darcy Units whose value is 9.869233 × 10−13 m2 . These units are in the I.S. Now, we will make each quantity dimensionless using a combination of parameters so that the order of the resultant variable is one. We will assume that a saturation field is given (usually obtained from transport step of the operator spliting technique). Note: Each quantity with superscript ∗ denotes a dimensionless quantity.

x∗ =

x , L

(A.2)

Fig. A1. Dirichlet/Neumann Boundary Value Problem for the Poisson equation. Left: The dimensional form. Right: The dimensionless form.

Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023

ARTICLE IN PRESS

JID: AMC 14

[m3Gsc;January 24, 2020;10:27]

A. Ali, H. Mankad and F. Pereira et al. / Applied Mathematics and Computation xxx (xxxx) xxx

p∗ =

p − pout p − pout = , pin − pout p

(A.3)

κ∗ =

κ , κmax

(A.4)

λ∗ =

λ , λmax

(A.5)

λmax κmax  p

u=

q=

L

λmax κmax  p L2

u∗ .

(A.6)

q∗ .

(A.7)

We note that the units of

λmax κmax  p

L using (A.2) - (A.7) in (A.1), we have,

are

⎧ ∗ ∗ u (x ) = −λ∗ κ ∗ (x∗ )∇ ∗ p∗ (x∗ ) ⎪ ⎪ ⎪ ∗ ∗ ∗ ∗ ∗ ⎪ ⎪ ⎨∇ · u (x ) = q (x ) p∗in = 1 on in , ∗ ∗ ⎪ p ( x ) = ⎪ ∗ ⎪ p = 0 on out , ⎪ out ⎪ ⎩ ∗ ∗ ˇ =0 u (x ) · n

m s

giving us the units for u and hence making u∗ a dimensionless quantity. Now,

in , in , (A.8) otherwise,

which is the required nondimensional form of (A.1). Note that in order to simplify the notation we have not written any of the quantities in the Sections 3 and 4 with a superscript ∗ , but they are all indeed in dimensionless form. References [1] T. Hou, X.-H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys. 134 (1) (1997) 169–189. [2] P. Jenny, S. Lee, H. Tchelepi, Multi-scale finite-volume method for elliptic problems in subsurface flow simulation, J. Comput. Phys. 187 (1) (2003) 47–67. [3] R. Guiraldello, R. Ausas, F. Sousa, F. Pereira, G. Buscaglia, The multiscale robin coupled method for flows in porous media, J. Comput. Phys. 355 (2018) 1–21. [4] R. Guiraldello, R. Ausas, F. Sousa, F. Pereira, G. Buscaglia, Interface spaces for the multiscale robin coupled method in reservoir simulation, Math. Comput. Simul. 164 (2019) 103–119. [5] T. Arbogast, G. Pencheva, M. Wheeler, I. Yotov, A multiscale mortar mixed finite element method, Multisc. Model. Simul. 6 (1) (2007) 319. [6] R. Araya, C. Harder, D. Paredes, F. Valentin, Multiscale hybrid-Mixed method, SIAM J. Numer. Anal. 51 (6) (2013) 3505–3531. [7] A. Francisco, V. Ginting, F. Pereira, J. Rigelo, Design and implementation of a multiscale mixed method based on a nonoverlapping domain decomposition procedure, Math. Comput. Simul. 99 (2014) 125–138. [8] J. Douglas, R. Ewing, M. Wheeler, A time-discretization procedure for a mixed finite element approximation of miscible displacement in porous media, RAIRO. Anal. Numérique 17 (3) (1983) 249–265. [9] J. Douglas, F. Furtado, F. Pereira, On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs, Comput. Geosci. 1 (2) (1997) 155–190. [10] J. Douglas, F. Pereira, L. Yeh, A locally conservative Eulerian–Lagrangian numerical method and its application to nonlinear transport in porous media, Comput. Geosci. 4 (1) (20 0 0) 1–40. [11] F. Furtado, V. Ginting, F. Pereira, M. Presho, Operator splitting multiscale finite volume element method for two-phase flow with capillary pressure, Transp. Porous Media 90 (3) (2011) 927–947. [12] C. Chen, G. Huan, Y. Ma, Computational Methods for Multiphase Flows in Porous Media, SIAM, 2006. [13] R. LeVeque, Finite Volume Methods for Hyperbolic Problems, 31, Cambridge University Press, 2002. [14] C.P. Ferraz, A novel recursive formulation of multiscale mixed methods and relaxation modeling of flow in porous media., University of Campinas, 2019 Ph.D. Thesis. [15] H. Akbari, A.P. Engsig-Karup, V. Ginting, F. Pereira, A multiscale direct solver for the approximation of flows in high contrast porous media, J. Comput. Appl. Math. 359 (2019) 88–101. [16] C. Holmes, Introduction to Perturbation Techniques, 20, 2, Springer-Verlag New York, 2013. [17] M. Christie, M. Blunt, et al., Tenth SPE comparative solution project: a comparison of upscaling techniques, in: Proceedings of the SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 2001. [18] D. Gamerman, H. Lopes, Markov Chain Monte Carlo - Stochastic simulation for Bayesian inference, in: Texts in Statistical Science, 68, 2, Chapman & Hall/CRC, 2006. [19] V. Ginting, F. Pereira, M. Presho, S. Wo, Application of the two-stage Markov chain monte carlo method for characterization of fractured reservoirs using a surrogate flow model, Comput. Geosci. 15 (2011) 691–707.

Please cite this article as: A. Ali, H. Mankad and F. Pereira et al., The multiscale perturbation method for second order elliptic equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.125023