- Email: [email protected]
Design and implementation of a multiscale mixed method based on a nonoverlapping domain decomposition procedure
Accepted Manuscript Title: Design and Implementation of a Multiscale Mixed Method Based on a Nonoverlapping Domain Decomposition Procedure Author: A. Francisco V. Ginting F. Pereira J. Rigelo PII: DOI: Reference:
S0378-4754(13)00199-7 http://dx.doi.org/doi:10.1016/j.matcom.2013.04.022 MATCOM 3973
To appear in:
Mathematics and Computers in Simulation
Received date: Revised date: Accepted date:
29-11-2011 26-8-2012 3-4-2013
Please cite this article as: A. Francisco, V. Ginting, F. Pereira, J. Rigelo, Design and Implementation of a Multiscale Mixed Method Based on a Nonoverlapping Domain Decomposition Procedure, Mathematics and Computers in Simulation (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.022 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.
Design and Implementation of a Multiscale Mixed Method Based on a Nonoverlapping Domain Decomposition Procedure
Departamento de Engenharia Mecˆanica, Universidade Federal Fluminense, Av. dos Trabalhadores, 420, Volta Redonda, RJ 27255-125, Brazil 2
University of Wyoming, Department of Mathematics, 1000 E. University Ave., Laramie, WY 82071-3036, U.S.A.
University of Wyoming, Department of Mathematics and School of Energy Resources, 1000 E. University Ave., Laramie, WY 82071-3036, U.S.A.
cr
3
us
1
ip t
A. Francisco 1 , V. Ginting 2 , F. Pereira 3 , J. Rigelo 2
an
Abstract
te
d
M
We use a nonoverlapping iterative domain decomposition procedure based on the Robin interface condition to develop a new multiscale mixed method to compute the velocity field in heterogeneous porous media. Hybridized mixed finite elements are used for the spatial discretization of the equations. We define local, multiscale mixed basis functions to represent the discrete solutions in subdomains. Appropriate subspaces of the vector space spanned by these basis functions can be considered in the numerical approximations of heterogeneous porous media flow problems. The balance between numerical accuracy and numerical efficiency is determined by the choice of these subspaces. A detailed description of the numerical method is presented. Following that, numerical experiments are discussed to illustrate the important features of the new procedure and its comparison to the traditional fine grid simulations.
Ac ce p
Keywords: Multiscale Methods, Mixed Finite Elements, Domain Decomposition.
1. Introduction
We are concerned with the development of numerical procedures for the fast and accurate approximation of subsurface flows that can take advantage of heterogeneous (e.g., CPU-GPU clusters) processing units. Such units are relatively inexpensive and have larger computational power (by one or two orders of magnitude) than CPUs alone. In the case of second order elliptic equations, multiscale procedures have the potential to fit well within the above mentioned heterogeneous computational environments because GPUs can be efficiently used to solve many small local problems that are posed by such procedures. We present a new multiscale method based on a nonoverlapping iterative domain decomposition for mixed finite elements introduced in [10]. Domain decomposition techniques for mixed methods distinct from the one used here can be found in, e.g., [5, 8, 12, 13, 16, 15].
Preprint submitted to Mathematics and Computers in Simulation
August 31, 2013
Page 1 of 26
cr
ip t
The development of multiscale numerical procedures for second order elliptic equations arising in porous media flow problems has attracted the attention of several research groups. We refer the reader to the recent publications [7, 11] and [17] for references and also a very comprehensive comparison of several procedures of this type. More recent developments include [19] along with [14] among many others. In [14] the authors make use of the concept of a multiscale flux basis. A similar concept will be used in this work. In the work presented here we define an iterative domain decomposition algorithm taking advantage of multiple grids and a new family of local basis functions to approximate large problems efficiently. The new procedure has some distinctive properties:
an
us
• Three length scales are introduced in the definition of the procedure: the solution is sought in the finest scale (h); multiscale basis functions are defined in the subdomains associated with the coarsest scale (H). Flux conservation is directly imposed on an intermediate scale (H with h ≤ H ≤ H). If H = h the new method produces the mixed finite element solution in the finest scale.
M
• The approximation is locally conservative at the H scale through the use of mixed finite elements. Discontinuous coefficients and source terms can be naturally handled.
d
• A post-processing step is introduced to recover the solution in the fine grid, given the solution with flux conservation at the H scale. A downscaling procedure ensures local conservation at the h scale.
te
• The implementation of the iterative scheme, which is based on the Robin interface conditions imposed on the intermediate scale, is straightforward and converges fast.
Ac ce p
• The procedure fits well in heterogeneous processing units (CPU-GPU): all the local problems can be efficiently solved in GPUs. These local problems are relatively small and fit in the GPU memory. Moreover, the iterative procedure associated with the global problem in a coarse grid takes place only at the boundaries of the coarse grid blocks, and can also run in GPUs. This paper is organized as follows. In section §2 we describe the mixed formulation of the elliptic problem considered in this work, next we give a step-by-step formulation of the new multiscale mixed method (MuMM), first introducing a discretization by hybridized mixed finite elements and a modified Robin condition, followed by the description of the MuMM iteration. In section §3 the computational implementation is described: multiscale basis functions of Robin type are introduced and a post-processing step is discussed. Section §4 provides several numerical results for the new method. Finally, section §5 contains our concluding remarks.
2
Page 2 of 26
2. The Multiscale Mixed Method (MuMM) 2.1. The Pressure Equation and its Variational Formulation
∇ · u = q,
ip t
For simplicity in the presentation, we consider Ω ⊂ <2 , a bounded domain with a Lipschitz boundary ∂Ω. The authors are currently extending this work to problems in <3 . The pressure-velocity system in the single-phase porous media flows ([6]) is given by u = −K(x)∇p,
(2.1)
cr
and
p = pb
u · ν = ub
on ΓD ,
us
where K(x) is known as the absolute permeability (a positive definite tensor), u is the Darcy velocity, p is the pressure and q = q(x) represents sources and sinks. Typical boundary conditions that occur in subsurface flow problems are Dirichlet and Neumann, which are respectively expressed as on ΓN ,
an
where ∂Ω = ΓD ∪ ΓN , ΓD ∩ ΓN = ∅ and ν is the outward unit normal vector. Our numerical procedure is derived from the weak formulation of the above pressure-velocity system which is described globally. To do so, we first define the following spaces:
M
W (Ω) = L2 (Ω),
H(div; Ω) = {v ∈ (L2 (Ω))2 | div v ∈ L2 (Ω)}, Vg (Ω) = {v ∈ H(div; Ω) | v · ν = g on ΓN },
te
d
for some function g. The global weak form of the pressure-velocity system (2.1) is given by finding {p, u} ∈ W × Vub such that ( div u, w)Ω = (q, w)Ω , ∀w ∈ W,
(2.2)
(x)u, v)Ω − (p, div v)Ω + hpb , v · νiΓD = 0, ∀v ∈ V0 ,
(2.3)
Ac ce p
(K
−1
where (·, ·)Ω is the L2 (Ω) inner product and h·, ·iΓ is the L2 (Γ) inner product involving line integration over Γ. We begin with the description of two partitions of Ω: the first is a uniform fine Cartesian mesh Ωf ; the second is a coarse Cartesian mesh Ωc constructed as sets of elements in Ωf (see Fig. 1). The proposed method is based on partitions Ωγ into rectangles {Ωγj , j = 1, . . . , Mγ , } (see Fig. 1), such that Mγ [ γ γ Ω = Ωj ; Ωγj ∩ Ωγk = ∅, j 6= k, γ = c, f. j=1
Define Γ = ∂Ω and, for j = 1, . . . , Mγ : Dγ
Γj
= ∂Ωγj ∩ ΓD ,
Nγ
Γj
= ∂Ωγj ∩ ΓN ,
Γγjk = Γγkj = ∂Ωγj ∩ ∂Ωγk ,
γ = c, f.
3
Page 3 of 26
Figure 1: The fine Ωf (dashed lines) and coarse Ωc (solid lines) partitions of Ω along with the three spatial scales used in the definition of the new multiscale procedure.
us
cr
ip t
As indicated in Fig. 1, we refer to three length scales in the description of the new multiscale procedure: H, the mesh size for the coarse partition; h, the mesh size of an underlying fine grid and H, an intermediate length scale (h ≤ H ≤ H). The Multiscale Mixed Method gives an approximate solution to (2.2-2.3) in that the continuity of the pressure and normal component of the flux are imposed on the intermediate scale H (in terms of averages of fine scale Lagrange multipliers and normal components of fluxes, respectively) along the edges of the coarse partition Ωc . Flux conservation at the finest scale is obtained after a downscaling step. The new method is based on a iterative procedure similar to the one defined in [10] with a modified Robin interface condition. To formulate the hybridized mixed finite element method over Ωf , we start by defining the local finite dimensional spaces W (Ωfj ) ⊂ L2 (Ωfj ) and V (Ωfj ) ⊂ H(div; Ωfj ). The Neumann boundary N
for that we set rj ∈ P% (Γj f ) to satisfy Nf
Γj
N
= 0,
N
∀ϕ ∈ P% (Γj f ),
M
hub − rj , ϕi
an
value ub needs to be projected appropriately in each of the segments ∂Ωfj that intersect with ΓN , and
(2.4)
N
d
where P% (Γj f ) is the space of polynomial of fixed degree % in Γj f . Then we set
te
N Vrj (Ωfj ) = v ∈ V (Ωfj ) | v · ν = rj on Γj f . At this stage we introduce Lagrange multipliers `jk on the edges {Γfjk } (see for example[1]). We set
Ac ce p
Λ = `jk : `jk |Γf ∈ P% (Γfjk ) = Λjk , Γfjk 6= ∅ .
(2.5)
jk
We note that there are two values of P% assigned to Γfjk , i.e., Λjk and Λkj . The hybridized mixed finite element method for the pressure-velocity system (2.2-2.3) is given by seeking for each j = 1, 2..., Mf {pj , uj , `jk } ∈ W (Ωfj ) × Ve (Ωfj ) × Λ, such that ( div uj , w)Ωf = (q, w)Ωf , j
j
∀w ∈ W (Ωfj ),
(2.6)
(K −1 uj , v)Ωf − (pj , div v)Ωf +
X k6=j
j
j
h`jk , v · νiΓf + hpb , v · νiΓ jk
f D ∩∂Ωj
= 0, ∀v ∈ V0 (Ωfj ),
(2.7)
f where Ve (Ωfj ) = Vrj (Ωfj ) if ΓN 6= ∅ and Ve (Ωfj ) = V (Ωfj ) otherwise and νjf is the outward normal j
unit vector on Ωfj . We note that the above formulation is local in nature, i.e., one that is posed in every 4
Page 4 of 26
Ωfj . To transform it to a global algebraic system, the following continuity conditions are imposed: (2.8)
uj · νj + uk · νk = 0 on Γfjk ,
(2.9)
ip t
`jk = `kj on Γfjk ,
an
us
cr
where νj (and similarly for νk ) is the unit normal vector pointing outward of Ωfj . For a typical situation, imposing these conditions allow for the construction of a global linear equations governing the Lagrange multipliers (see for example, [1, 4]). In relative comparison to the traditional mixed finite element method, this hybridization technique has a clear advantage of reducing the system dimension. Another advantage is that the resulting global matrix from the hybridization is positive definite. Still, for problems such as multiphase flow in heterogeneous porous media, construction of a global matrix is prohibitive and may have to be avoided. In an attempt to provide an alternative along this direction, the authors of [10] developed an iterative procedure that is based on the local calculation (2.6) and (2.7) along with the imposition of a Robin transmission boundary conditions −χjk uj · νj + `jk = χjk uk · νk + `kj ,
on Γfjk ⊂ ∂Ωfj ,
(2.10)
−χjk uk · νk + `kj = χjk uj · νj + `jk ,
on Γfkj ⊂ ∂Ωfk ,
(2.11)
Ac ce p
2.2. Local problems on Ωcm
te
d
M
where χjk is a positive function on Γfjk . These Robin conditions are in fact equivalent to (2.8) and (2.9). In [10], the core of the iterative procedure was built by substituting `jk in (2.7) by the one in (2.10) with uk and `kj treated as known previous iteration values from the neighboring elements. The multiscale method we propose in this paper follows a closely related iterative procedure. The main distinction is that our method employs a red-black iterative procedure and an inclusion of information from intermediate scale through appropriate construction of local problems. This is the subject of the next subsection.
To describe the multiscale procedure for approximating the above formulation we construct a set of local problems on Ωcm , m = 1, · · · , Mc with Robin boundary conditions on ∂Ωcm . For each m = 1, · · · , Mc , we write the variational formulation of the pressure-velocity system (2.1) as to seek {pm , um } ∈ L2 (Ωcm ) × H(div; Ωcm ) such that ( div um , w)Ωcm = (q, w)Ωcm , ∀w ∈ L2 (Ωcm ),
(K
−1
m
(2.12)
m
(x)u , v)Ωcm − (p , div v)Ωcm
+ aum · ν m + g, v · ν m ∂Ωc = 0, ∀v ∈ H(div; Ωcm ), m
(2.13)
where we have imposed a Robin boundary condition on ∂Ωcm which is expressed in a generic format −aum · ν m + pm = g,
5
Page 5 of 26
for given functions a and g, and ν m the unit normal vector on ∂Ωcm pointing outward of Ωcm . We note that in the case Ωcm is adjacent to either of the global boundary, then we need to modify (2.13) appropriately to honor that particular global boundary.
ip t
Figure 2: Illustration of localization on Ωcm for the iterative procedure
us
cr
Fig. 2 outlines a schematic of the localization of (2.12) and (2.13) in Ωcm . Here, without loss of generality, we consider an internal Ωcm whose adjacent neighbors are Ωcn , where n = l, r, b, t. As described earlier, the key to the construction of the iterative procedure is the localization of (2.12) and (2.13) through a careful choice of Robin boundary condition on ∂Ωcm . In particular, the Robin boundary condition in (2.13) can be decomposed into X
m m au · ν + g, v · ν m ∂Ωc = aum · ν m + g, v · ν m Γc , (2.14) m
n=l,r,b,t
mn
an
where each Γcmn is as expressed in Fig. 2. Viewed this way, the natural choice for g in each term is a functional that depends on the neighboring element values, i.e., g c = g(un , pn ),
M
Γmn
d
where un and pn are velocity and trace pressure of Ωcn calculated at Γcmn . A detailed explanation of the modified Robin boundary condition is the subject of the next subsection. 2.3. The modified Robin condition
Ac ce p
te
To proceed, we set H and H such that r = H/H is an integer. As depicted in Fig. 2, a Γcmn is i , i = 1, · · · , r, each having length H. Thus decomposed into a collection of line segments Imn ∂Ωcm
=
[
Γcmn
=
n=l,r,b,t
[
r [
n=l,r,b,t
i=1
i Imn ,
i i = nψ = 4r. For each i, we define on Imn ⊂ Γcmn , which implies that dim Imn Z 1 i g i = g mn = (χimn un · ν n + pn )dl, Imn i H Imn
(2.15)
i . Hence, each term in the summation on where we assume that χmn is a positive constant on each Imn the right hand side of (2.14) can be expressed as
m m au · ν + g, v · ν m Γc
=
r X
mn
χimn um · ν m + g imn , v · ν m
i Imn
.
(2.16)
i=1
6
Page 6 of 26
We note that imposing this condition is equivalent to f m m i i −χimn um j · ν + pjk = g mn , ∀ Γjk ⊂ Imn ,
(2.17)
f m m where Γfjk = ∂Ωfjm ∩ ∂Ωfkn (see Fig. 2), um j = u Ωf , and pjk is the trace pressure at Γjk .
ip t
jm
As for the numerical approximation of (2.12-2.13), we use the hybridized mixed finite element approximation that is akin to (2.6) and (2.7).
cr
2.4. The MuMM iteration
M
an
us
We now define an iterative procedure that localizes the calculations and gives an approximate solution to the global problem (2.6-2.7). Consider the domain Ω and label its coarse partitions as “red” and “black” in an ordinary checkerboard order. The iterative process is listed in Algorithm 1 where we have assumed to use the hybridized mixed finite element approximation for solving the local problem (2.12) and (2.13). In essence, either red or black calculation in Algorithm 1 is represented by the schematic in Fig. 2. For example, in the red calculation, Ωcm is the red subdomain, while Ωcn , n = l, r, b, t are the black subdomains from which the boundary information are gathered. Similar setting is also true for black calculation. We remark that in the case of H = H = h and Ωf = Ωc the convergence of this iterative scheme has been established in [10]. Algorithm 1 MuMM Iteration
te
d
Set initial guess for u(0) , `(0) . for α = 1 until convergence do Red Subdomain: • For each red element, compute p(α) , u(α) , `(α) by solving (2.12-2.13)
Ac ce p
• Use the (α − 1) iterate from the black subdomains nearest neighbors in the right hand side of (2.14) Black Subdomain: • For each black element, compute p(α) , u(α) , `(α) by solving (2.12-2.13) • Use α iterate from the red subdomains nearest neighbors in the right hand side of (2.14) end for
3. Computational Implementation
In our numerical experiments we shall adopt the lowest order Raviart-Thomas elements and we consider this space in the discussion of the new scheme. The natural degrees of freedom on the element Ωfj , where the elements are squares of side length h, are the (constant) value pj of the pressure over Ωfj , which can be interpreted as the value at the center of the element, and the four (constant) values ujβ , β = l, r, b, t, of the outward normal component of the total flux across the edges of the element.
7
Page 7 of 26
us
cr
ip t
The Lagrange multipliers `jk and `kj on Γfjk are also constants. We shall also assume that the absolute permeability K(x) and the external flow rate q to be constant on each element. However, it is obvious that Algorithm 1 can easily become an expensive calculation especially when the dimension of the local calculation in Ωcm increases significantly. This is even more pronounced when it is applied to multiphase flow simulations that require repeated calculations of pressure and velocity and hence repeated algorithm calls. For that matter, it is imperative to devise an innovative way to be incorporated in Algorithm 1 that allows for efficient calculation and yet maintains an acceptable accuracy. Our strategy in this direction is to construct a set of multiscale basis functions associated with Ωcm . Calculation of these basis functions is a one time preprocessing step and thus is reasonable to expend. Furthermore, with these basis functions in place, the iteration procedure of Algorithm 1 is only conducted on ∂Ωcm . In the next subsection, we describe the aforementioned basis functions.
M
an
3.1. Mixed Multiscale Basis Functions of Robin Type As mentioned, the multiscale basis functions are associated with Ωcm , m = 1, ·i· · , Mc . Moreover, m the number of multiscale basis functions for each Ωc depends on how many Imn are employed. Thus the total number of basis functions for each Ωm c is nψ = 4r, where as described earlier, r = H/H. This means there are Mc nψ of multiscale basis functions constructed on the whole domain Ω. As in Algorithm 1, the local problem on Ωm c for the multiscale basis functions areclosed with a i Robin boundary condition. One Imn is tied to one basis function, which we denote by ψ imn , φimn . These functions satisfy variational formulation that is similar to (2.12) and (2.13) (with q = 0). Furthermore, the Robin boundary condition is set to satisfy (3.1)
is as depicted in Fig. 3, namely
te
where
i gmn
d
i −χimn ψ imn · ν m + φimn = gmn in ∂Ωcm ,
( i 1, in Imn = i . 0, in ∂Ωcm \Imn
Ac ce p
i gmn
Taking all these into account we arrive at the following variational formulation: find the divergence i,H = ψ imn , φimn satisfying free Mixed Multiscale Basis Functions (MMBFs) Ψmn ( div ψ imn , w)Ωcm = 0, ∀w ∈ L2 (Ωcm ), (K
(3.2)
−1
(x)ψ imn , v)Ωcm − (φimn , div v)Ωcm
i + χimn ψ imn · ν m + gmn , v · ν m ∂Ωc m
= 0, ∀v ∈ H(div; Ωcm ).
(3.3)
We note that as in Algorithm 1, this variational formulation can be approximated using the hybridized mixed finite element. If x ∈ ΓD or x ∈ ΓN this construction is modified to take trivial Dirichlet or Neumann conditions. Moreover, to the set of basis functions associated with a Ωcm that is adjacent to the external boundary, an extra MMBF has to be computed, with the prescribed external boundary condition and zero value for the Robin condition in the remaining boundaries. We remark that in [14] a similar set of basis functions was defined with distinct boundary conditions. 8
Page 8 of 26
i H Figure 3: One line segment Imn where the Robin boundary condition for ψmi is imposed.
S Ωcm =
X
r X
Aimn Ψi,H mn ,
(3.4)
cr
n=l,r,b,t i=1
ip t
3.2. The MuMM iteration revisited In performing the MuMM iteration, the multiscale basis functions described above have to be pre-computed. Within the iteration updating the local solution for each Ωcm consists of computing
an
us
where S Ωcm indicates the approximation of um , pm . The coefficients Aimn in (3.4) are given by the (previous iteration) nearest neighbor modified Robin condition (see (2.15)); there is no linear algebra problem to be solved in each iteration. Note that in a computational implementation of (3.4) the sum should be evaluated only at the boundaries of each Ωcm . After convergence of the MuMM iteration a downscaling step to be discussed next will produce the solution inside Ωcm . We remark that:
M
• The approximate solution for the global problem (2.2-2.3) is then given by Mc X
S=
S Ωcm .
(3.5)
m=1
d
• Flux conservation is maintained in the H scale.
Ac ce p
te
• The balance between numerical accuracy and numerical efficiency is determined by the choice i,h of span{Ψi,H mn } ⊂ span{Ψmn }. The closer H is to h the larger is the number of MMBFs to be computed and the closer the MuMM solution is to the fine grid solution. 3.3. Downscaling As mentioned above, after convergence the flux conservation is maintained in the H scale. In order to have the flux conservation in the h scale, we employ a two-step procedure. We first solve local Neumann problems to recover the solution inside each element of Ωc . The boundary condition for these problems specified in the h scale is set to be the average flux produced by the MuMM iteration at the boundaries between two neighbor elements of Ωc . Next we solve a family of local Neumann problems in a staggered grid. Fig. 4 illustrates this procedure. After this step we recover an approximate solution in the fine scale. The computational cost of the downscaling is equivalent to the contruction of one MMBF, in both cases we solve one elliptic problem for each subdomain. R EMARK : Although the procedure has been defined for uniform structured grids, we believe that it could be implemented on more complicated (non-uniform or unstructured) grids. Such implementations are, however, outside the scope of this work.
9
Page 9 of 26
Figure 4: The solid boundary corresponds to a cell of Ωc and indicates the region where a local Neumann problem is solved in step-1 of the downscaling procedure. In the second step of the procedure the local solution is computed in the region indicated by a dashed boundary.
ip t
4. Numerical Results
φ
an
us
cr
In order to show that MuMM is a viable alternative to solve reservoir simulation problems, we present numerical results for the simulation of two-dimensional fluid flows associated with two cases, namely, quarter of a 5-spot problem and a slab geometry on a channelized field. In all simulations, we compare the MuMM with the usual traditional simulation on the mesh in which the fine scale permeability field is resolved. Aspects of comparison include the accuracy of the approximate velocity, efficiency of MuMM in terms of number of iterations to reach a convergence tolerance, and the robustness of MuMM when incorporated in a single phase flow simulation. With respect to this single phase flow simulation, the model governs the physical transport of fluids that uses the approximate velocity as an input to the equation ∂c + u · ∇c = 0. ∂t
(4.1)
te
d
M
This equation is solved using the genuinely multidimensional second-order central scheme for conservation laws described in described in [18]. Here φ is the porosity and c is the concentration which is defined on the fine mesh. We consider a constant c(x, t = 0) for initial condition. In practice, it is of interest to gather the concentration at a dimensionless time PVI, whose value is Z t Z PVI = Vp −1 u · ν dl dτ, (4.2) 0
∂Ωout
Ac ce p
for a given dimensional time t. Here Vp is the total pore-volume of the reservoir. 4.1. Example 1: The quarter of a 5-spot problem This is a standard problem in reservoir simulation and thus it is an important example to evaluate the performance of the MuMM. In this simulation the global fine grid is of 128 × 128 elements in the region (0, 1) × (0, 1). The global boundary condition is of Neumann type (with no flow on the boundary), and there is an injection through a well at the lower left corner of the computational domain, and a production, at the same rate, through a well in the upper right corner. The experiment uses a 2-D heterogeneous permeability whose log value field ξ(x) is shown in Fig. 5. We set k(x) = exp(δ ξ(x)), where δ is chosen in order to set a high permeability ratio kmax /kmin = 106 . The field ξ(x) is a self similar (Gaussian), (see [2, 3]), assumed to be stationary and isotropic, and is characterized by its (zero) mean and covariance function C(x, y) given by: C(x, y) = |x − y|−0.5 .
(4.3)
10
Page 10 of 26
Figure 5: Log permeability field used in the quarter of a 5-spot problem.
ip t
Figure 6: Comparison of concentration surface profiles in Example 1, where the velocity is from direct fine grid simulation (left) and from MuMM (right). The time in PVI is, from top to bottom, 0.32 and 0.65.
k uMuMM − ufine k j j m
M
an
us
cr
In Fig. 6 we show concentration profiles at two different PVIs, in which the velocity used in the left plots is from the direct fine grid calculation using the iterative procedure from [10], while its counterpart in the right plots is from the MuMM with coarse grid of 4 × 4 elements (i.e., H = 1/4, and thus making each coarse element Ωcm to possess 32 × 32 fine elements). We use H = 1/16, which gives dim{Ψi,H mn } = 4H/H = 4 · 4 = 16, i.e., there are 16 MMBFs per coarse element. An additional basis function has to be constructed for lower left and upper right coarse element in which the injection and production well is respectively located. Each of these extra basis functions is defined with trivial Robin conditions that takes the source/sink terms into account. The comparison in Fig. 6 indicates that the MuMM produces a good quality numerical approximation of the solution in the fine grid. As described, the aforementioned simulations take the approximate velocity on the fine grid as an input, so maintaining a satisfactory level of its accuracy is imperative. Further investigation of the relative error spatial profile for MuMM approximation depicted in Fig. 7 confirms the previous finding. The relative error for each fine element Ωfj , j = 1, ..., Mf is evaluated as , where m = max k ufine k, j
(4.4)
d
j
Ac ce p
te
where k . k is the usual Euclidean norm. Simulation to create this figure uses a fixed H = 1/16 and varying H = 1/16, 1/8, 1/4, respectively having 8 × 8, 16 × 16, and 32 × 32 fine elements in each coarse element. Consequently, we maintain a fixed total number of 64 MMBFs on the whole domain. We note that the maximum relative error in each of the plots from left to right is respectively 0.15, 0.08, and 0.019. This result suggests that when the total number of MMBFs remains fixed on the whole domain, the accuracy is better when the H scale is coarser but at the expense of having a larger degrees of freedom in each multiscale basis function calculation. To complete the study of the MuMM for the quarter of 5-spot problem, we present in Table 1 the number of iterations of the method for various ratios of permeability along with two choices of coarse grids. Note that the number of iterations increases slowly as a function of the strength of the permeability field. The criterion to stop iteration is based on how close the flux conservation is satisfied in the H scale.
Figure 7: The relative errors of velocity obtained from MuMM in Example 1 for various H: 1/16 (left), 1/8 (middle), 1/4 (right). The value of H is fixed at 1/16.
11
Page 11 of 26
103 104 105 106
Coarse grid of 8×8 56 63 74 91
Coarse grid of 4×4 24 26 29 34
cr
kmax /kmin
ip t
Table 1: Number of iterations in MuMM for various ratios of permeability and coarse grid settings in Example 1.
us
4.2. Example 2: Slab geometry on a channelized field
Ac ce p
te
d
M
an
The complex features of a channelized permeability field makes this example very challenging for testing multiscale methods, so that it is important to consider such a case in the evaluation of the MuMM performance. We use Ω = (0, 3 23 ) × (0, 1) which is discretized into 220 × 60 fine elements on which the permeability data is posed. The MuMM employs a coarse grid of 11 × 3 implying H = 1/3. This way, each coarse element possesses 20 × 20 fine elements. The left and right global boundary conditions are of Dirichlet type, and the top and bottom global boundary conditions are of homogeneous Neumann type (no flow). Fig. 8 shows the spatial profile ξ(x), which is the log permeability of layer 36 of SPE10 benchmark1 . The permeability is given by k(x) = exp(δ ξ(x)), where we adjust the parameter δ to give us permeability ratios of interest to our study. In the first case we consider the permeability ratio of 176 and three choices for the H scale, namely, 1/3, 1/6, and 1/12, respectively giving 4, 8 and 16 MMBFs in each coarse element. In Fig. 9 we compare tracer injection profiles for these three choices of H with the fine grid solution (computed using the iterative procedure from [10]). We note that with only 16 MMBFs per coarse element, the MuMM iteration can recover a very similar profile to the one obtained from a direct fine grid simulation. Furthermore, Fig. 10 compares the tracer cut curves, which is the fraction of the tracer in the produced fluid represented as Z cu · ν dl ∂Ωout F (t) = 1 − Z u · ν dl ∂Ωout
where ∂Ωout denotes outflow boundary. The tracer cut is plotted against the PVI given in (4.2). Next, to illustrate the importance of downscaling step within MuMM, we use 16 MMBFs per coarse element (H = H/4) and kmax /kmin = 104 . We calculate the velocity using MuMM with and without the downscaling step and make comparison against the direct fine grid simulation whose profile is plotted in Fig. 11. We note on the dominance of the flow in the regions that exhibit the channel 1
More information at http://www.spe.org/csp
12
Page 12 of 26
Figure 8: The log permeability field from the SPE10 model (layer 36)
ip t
Figure 9: Tracer concentrations of (from top to bottom): MuMM for H=H , H=H/2, H=H/4 and the fine grid, in Example 2. Here H = 1/3 and time (in PVI) is 0.77.
an
us
cr
feature. In particular, we zoom on the region confined by the white lines (see Fig. 11) and compare the respective velocities in Fig. 12. As can be seen from that figure, the downscaling procedure improves the numerical solution. Finally we conclude the analysis of MuMM for Example 2 with Table 2, which shows the iterations of MuMM for various permeability ratios. The cost per iteration is the cost of computing (3.4) for each subdomain. These results indicate the efficiency of the MuMM to capture the features of the channelized permeability field. As in Example 1, the criterion to stop iteration is based on how close the flux conservation is satisfied. R EMARK : In line with [9], in expressions (2.10 -2.11) we set
2Kj Kk ch H , where Kjk = , H Kj + Kk Kjk
M
χjk =
and c is a dimensionless constant. Within the modified Robin condition (2.15) we set
te
d
χimn =
ch i
K jk
i
Ac ce p
H set in the line for the log normal permeability field, where K jk is the average value among Kjk i . For the channelized field we set χi segment Imn mn to be the maximum value among χjk previously i . determined in the line segment Imn
5. Concluding Remarks
In this paper, we have laid out a framework for efficient numerical approximation of the velocity in the presence of heterogeneous media using a mixed method that is incorporated into a nonoverlapping domain decomposition procedure. The main attraction of the proposed framework is its immediate parallelizability owing to the construction of the multiscale basis functions within each coarse elements. On the global level, an iterative procedure is devised which only involves calculation on the edges of the coarse elements. This iteration avoids the construction and solving of linear algebra.
Figure 10: Tracer cut curve profiles for MuMM and the fine grid from Example 2.
13
Page 13 of 26
kmax /kmin 102 103 104
cr
MuMM Iterations 41 42 88
ip t
Table 2: Number of iterations using different ratios of permeability in Example 2:
us
Figure 11: The velocity field with the log permeability on the background. See Fig. 12 for comparison of MuMM with fine grid simulations on the zoomed region (a box confined with the white lines).
an
After the downscaling step we recover a good approximation to the solution in the fine scale, and mass conservation holds also in that scale. Sources and sinks are naturally incorporated in the procedure, by considering an extra MMBF (constructed with trivial Robin conditions) that takes into account the source or sink in the coarse partition element where it is located. The authors are currently implementing the 3D version of the procedure described here in a CPUGPU cluster. Moreover, the numerical analysis of the MuMM is underway.
M
Acknowledgments
Ac ce p
te
d
The authors wish to thank Jim Douglas, Jr. for some discussions about the method presented in the this paper and Dr. Marcos Mendes for his help with some numerical simulations. V. Ginting and F. Pereira are supported in part by DOE DE-FE0004832 and by grants from the Center for Fundamentals of Subsurface Flow of the School of Energy Resources of the University of Wyoming, under contract numbers WYDEQ49811GNTG and WYDEQ49811PER (Ginting) and WYDEQ49811GNTG, WYDEQ49811PER and WYDEQ49811FRTD (Pereira). Ginting is also supported in part by grants from US Department of Energy (DE-SC0004982) and the National Science Foundation (DMS-1016283). Pereira is also supported in part by the Wyoming State Legislature, account number 1100 20352 2012. A. Francisco was supported by grant from CAPES/FULBRIGHT Program (BEX 2984/09-1). J. Rigelo is supported by DOE DE-FE0004832. References
[1] D. N. Arnold, F. Brezzi, Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, RAIRO Mod´el. Math. Anal. Num´er. 19 (1) (1985) 7–32.
Figure 12: A zoomed window selected from the domain in Fig. 11 above (from left to right): fine grid, MuMM iteration before and after downscaling. Note the importance of the downscaling procedure to correct the flux within the channel.
14
Page 14 of 26
[2] M. R. Borges, F. Furtado, F. Pereira, H. P. Amaral Souto, Scaling analysis for the tracer flow problem in self-similar permeability fields, Multiscale Model. Simul. 7 (3) (2008) 1130–1147.
ip t
[3] M. R. Borges, F. Pereira, H. P. Amaral Souto, Efficient generation of multi-scale random fields: a hierarchical approach, Int. J. Numer. Methods Biomed. Eng. 26 (2) (2010) 176–189.
cr
[4] G. Chavent, J. Roberts, A unified physical presentation of mixed, mixed-hybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems, Adv. Water Resour. 14 (6) (1991) 329–348.
us
[5] W. Chen, M. Gunzburger, F. Hua, X. Wang, A parallel Robin-Robin domain decomposition method for the Stokes-Darcy system, SIAM J. Numer. Anal. 49 (3) (2011) 1064–1084.
an
[6] Z. Chen, G. Huan, Y. Ma, Computational methods for multiphase flows in porous media, Computational Science & Engineering, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2006. [7] J. Chu, Y. Efendiev, V. Ginting, T. Hou, Flow based oversampling technique for multiscale finite element methods, Adv. Water Resour. 31 (4) (2008) 599–608.
d
M
[8] L. C. Cowsar, M. F. Wheeler, Parallel domain decomposition method for mixed finite elements for elliptic partial differential equations, in: Fourth International Symposium on Domain Decomposition Methods for Partial Differential Equations (Moscow, 1990), SIAM, Philadelphia, PA, 1991, pp. 358–372.
te
[9] J. Douglas, Jr., F. Furtado, F. Pereira, On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs, Comput. Geosci. 1 (2) (1997) 155–190.
Ac ce p
[10] J. Douglas, Jr., P. J. Paes-Leme, J. E. Roberts, J. P. Wang, A parallel iterative procedure applicable to the approximate solution of second order partial differential equations by mixed finite element methods, Numer. Math. 65 (1) (1993) 95–108. [11] Y. Efendiev, T. Y. Hou, Multiscale finite element methods, vol. 4 of Surveys and Tutorials in the Applied Mathematical Sciences, Springer, New York, 2009, theory and applications. [12] R. E. Ewing, J. Wang, Analysis of the Schwarz algorithm for mixed finite elements methods, RAIRO Mod´el. Math. Anal. Num´er. 26 (6) (1992) 739–756. [13] R. E. Ewing, J. Wang, Analysis of multilevel decomposition iterative methods for mixed finite element methods, RAIRO Mod´el. Math. Anal. Num´er. 28 (4) (1994) 377–398. [14] B. Ganis, I. Yotov, Implementation of a mortar mixed finite element method using a multiscale flux basis, Comput. Methods Appl. Mech. Engrg. 198 (49-52) (2009) 3989–3998.
15
Page 15 of 26
[15] R. Glowinski, W. Kinton, M. F. Wheeler, Acceleration of domain decomposition algorithms for mixed finite elements by multi-level methods, in: Third International Symposium on Domain Decomposition Methods for Partial Differential Equations (Houston, TX, 1989), SIAM, Philadelphia, PA, 1990, pp. 263–289.
ip t
[16] R. Glowinski, M. F. Wheeler, Domain decomposition and mixed finite element methods for elliptic problems, in: First International Symposium on Domain Decomposition Methods for Partial Differential Equations (Paris, 1987), SIAM, Philadelphia, PA, 1988, pp. 144–172.
cr
[17] V. Kippe, J. E. Aarnes, K.-A. Lie, A comparison of multiscale methods for elliptic problems in porous media flow, Comput. Geosci. 12 (3) (2008) 377–398.
us
[18] F. Pereira, A. Rahunanthan, A semi-discrete central scheme for the approximation of two-phase flows in three space dimensions, Math. Comput. Simulation 81 (10) (2011) 2296–2306.
Ac ce p
te
d
M
an
[19] M. F. Wheeler, T. Wildey, G. Xue, Efficient algorithms for multiscale modeling in porous media, Numer. Linear Algebra Appl. 17 (5) (2010) 771–785.
16
Page 16 of 26
d
te
Ac ce p us
an
M
cr
H
h
H
Page 17 of 26
ip t
1
ip t
Imn 2
Ac ce p
te
d
M
an
us
cr
Imn
Page 18 of 26
ip t cr us an M d te Ac ce p
0
0
1
Imn 1
0
0
0 0
0
1 Figure 3: One line segment Imn where the Robin boundary condition for ψ
Page 19 of 26
Page 20 of 26
d
te
Ac ce p us
an
M
cr
ip t
ip t cr us an
120
M
100
80
te
40
d
60
20
40
60
80
100
120
Ac ce p
20
Figure 5: Permeability field used in example 1 with permeability ratio of kmax /kmin = 106 . The red and blue colors refer to the highest and lowest permeability values, respectively .
13
Page 21 of 26
ip t cr us an M d te
Ac ce p
Figure 6: Concentration surface profiles of MuMM (left side) and the fine grid (right side), in Example 1. Note the good agreement between the solutions for a large value of kmax /kmin = 106
14
Page 22 of 26
−3
x 10 18 120
0.02
120 16
ip t
0.018 100
100
14
0.016 0.014
80
12
80
0.012 0.01
10
60
8
cr
60
0.008 40
40 0.006 0.004
20
20
40
60
80
100
120
20
40
60
80
100
4 2
120
Ac ce p
te
d
M
an
20
us
0.002
6
Page 23 of 26
ip t cr us an M d te Ac ce p Figure 8: The permeability field from the SPE10 model (layer 36) with permeability ratio of 176 used in the second example.
16
Page 24 of 26
ip t cr us an M d te Ac ce p
Figure 9: Tracer concentrations of: MuMM for H=H (upper left), H=H/2 (upper right), H=H/4 (bottom left) and the fine grid (bottom right), in Example 2.
Figure 10: Tracer cut curve profiles for MuMM and the fine grid from Example 2. A good agreement is observed for H ≤ H/2.
17
Page 25 of 26
ip t cr us an M d te Ac ce p
Figure 11: The permeability field together with the velocity field for (from top to bottom): fine grid, MuMM iteration before and after downscaling.
Figure 12: A window selected from the domain in Fig. 11 above (from left to right): fine grid, MuMM iteration before and after downscaling. Note the importance of the downscaling procedure to correct the flux within the channel.
18
Page 26 of 26
S0378-4754(13)00199-7 http://dx.doi.org/doi:10.1016/j.matcom.2013.04.022 MATCOM 3973
To appear in:
Mathematics and Computers in Simulation
Received date: Revised date: Accepted date:
29-11-2011 26-8-2012 3-4-2013
Please cite this article as: A. Francisco, V. Ginting, F. Pereira, J. Rigelo, Design and Implementation of a Multiscale Mixed Method Based on a Nonoverlapping Domain Decomposition Procedure, Mathematics and Computers in Simulation (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.022 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.
Design and Implementation of a Multiscale Mixed Method Based on a Nonoverlapping Domain Decomposition Procedure
Departamento de Engenharia Mecˆanica, Universidade Federal Fluminense, Av. dos Trabalhadores, 420, Volta Redonda, RJ 27255-125, Brazil 2
University of Wyoming, Department of Mathematics, 1000 E. University Ave., Laramie, WY 82071-3036, U.S.A.
University of Wyoming, Department of Mathematics and School of Energy Resources, 1000 E. University Ave., Laramie, WY 82071-3036, U.S.A.
cr
3
us
1
ip t
A. Francisco 1 , V. Ginting 2 , F. Pereira 3 , J. Rigelo 2
an
Abstract
te
d
M
We use a nonoverlapping iterative domain decomposition procedure based on the Robin interface condition to develop a new multiscale mixed method to compute the velocity field in heterogeneous porous media. Hybridized mixed finite elements are used for the spatial discretization of the equations. We define local, multiscale mixed basis functions to represent the discrete solutions in subdomains. Appropriate subspaces of the vector space spanned by these basis functions can be considered in the numerical approximations of heterogeneous porous media flow problems. The balance between numerical accuracy and numerical efficiency is determined by the choice of these subspaces. A detailed description of the numerical method is presented. Following that, numerical experiments are discussed to illustrate the important features of the new procedure and its comparison to the traditional fine grid simulations.
Ac ce p
Keywords: Multiscale Methods, Mixed Finite Elements, Domain Decomposition.
1. Introduction
We are concerned with the development of numerical procedures for the fast and accurate approximation of subsurface flows that can take advantage of heterogeneous (e.g., CPU-GPU clusters) processing units. Such units are relatively inexpensive and have larger computational power (by one or two orders of magnitude) than CPUs alone. In the case of second order elliptic equations, multiscale procedures have the potential to fit well within the above mentioned heterogeneous computational environments because GPUs can be efficiently used to solve many small local problems that are posed by such procedures. We present a new multiscale method based on a nonoverlapping iterative domain decomposition for mixed finite elements introduced in [10]. Domain decomposition techniques for mixed methods distinct from the one used here can be found in, e.g., [5, 8, 12, 13, 16, 15].
Preprint submitted to Mathematics and Computers in Simulation
August 31, 2013
Page 1 of 26
cr
ip t
The development of multiscale numerical procedures for second order elliptic equations arising in porous media flow problems has attracted the attention of several research groups. We refer the reader to the recent publications [7, 11] and [17] for references and also a very comprehensive comparison of several procedures of this type. More recent developments include [19] along with [14] among many others. In [14] the authors make use of the concept of a multiscale flux basis. A similar concept will be used in this work. In the work presented here we define an iterative domain decomposition algorithm taking advantage of multiple grids and a new family of local basis functions to approximate large problems efficiently. The new procedure has some distinctive properties:
an
us
• Three length scales are introduced in the definition of the procedure: the solution is sought in the finest scale (h); multiscale basis functions are defined in the subdomains associated with the coarsest scale (H). Flux conservation is directly imposed on an intermediate scale (H with h ≤ H ≤ H). If H = h the new method produces the mixed finite element solution in the finest scale.
M
• The approximation is locally conservative at the H scale through the use of mixed finite elements. Discontinuous coefficients and source terms can be naturally handled.
d
• A post-processing step is introduced to recover the solution in the fine grid, given the solution with flux conservation at the H scale. A downscaling procedure ensures local conservation at the h scale.
te
• The implementation of the iterative scheme, which is based on the Robin interface conditions imposed on the intermediate scale, is straightforward and converges fast.
Ac ce p
• The procedure fits well in heterogeneous processing units (CPU-GPU): all the local problems can be efficiently solved in GPUs. These local problems are relatively small and fit in the GPU memory. Moreover, the iterative procedure associated with the global problem in a coarse grid takes place only at the boundaries of the coarse grid blocks, and can also run in GPUs. This paper is organized as follows. In section §2 we describe the mixed formulation of the elliptic problem considered in this work, next we give a step-by-step formulation of the new multiscale mixed method (MuMM), first introducing a discretization by hybridized mixed finite elements and a modified Robin condition, followed by the description of the MuMM iteration. In section §3 the computational implementation is described: multiscale basis functions of Robin type are introduced and a post-processing step is discussed. Section §4 provides several numerical results for the new method. Finally, section §5 contains our concluding remarks.
2
Page 2 of 26
2. The Multiscale Mixed Method (MuMM) 2.1. The Pressure Equation and its Variational Formulation
∇ · u = q,
ip t
For simplicity in the presentation, we consider Ω ⊂ <2 , a bounded domain with a Lipschitz boundary ∂Ω. The authors are currently extending this work to problems in <3 . The pressure-velocity system in the single-phase porous media flows ([6]) is given by u = −K(x)∇p,
(2.1)
cr
and
p = pb
u · ν = ub
on ΓD ,
us
where K(x) is known as the absolute permeability (a positive definite tensor), u is the Darcy velocity, p is the pressure and q = q(x) represents sources and sinks. Typical boundary conditions that occur in subsurface flow problems are Dirichlet and Neumann, which are respectively expressed as on ΓN ,
an
where ∂Ω = ΓD ∪ ΓN , ΓD ∩ ΓN = ∅ and ν is the outward unit normal vector. Our numerical procedure is derived from the weak formulation of the above pressure-velocity system which is described globally. To do so, we first define the following spaces:
M
W (Ω) = L2 (Ω),
H(div; Ω) = {v ∈ (L2 (Ω))2 | div v ∈ L2 (Ω)}, Vg (Ω) = {v ∈ H(div; Ω) | v · ν = g on ΓN },
te
d
for some function g. The global weak form of the pressure-velocity system (2.1) is given by finding {p, u} ∈ W × Vub such that ( div u, w)Ω = (q, w)Ω , ∀w ∈ W,
(2.2)
(x)u, v)Ω − (p, div v)Ω + hpb , v · νiΓD = 0, ∀v ∈ V0 ,
(2.3)
Ac ce p
(K
−1
where (·, ·)Ω is the L2 (Ω) inner product and h·, ·iΓ is the L2 (Γ) inner product involving line integration over Γ. We begin with the description of two partitions of Ω: the first is a uniform fine Cartesian mesh Ωf ; the second is a coarse Cartesian mesh Ωc constructed as sets of elements in Ωf (see Fig. 1). The proposed method is based on partitions Ωγ into rectangles {Ωγj , j = 1, . . . , Mγ , } (see Fig. 1), such that Mγ [ γ γ Ω = Ωj ; Ωγj ∩ Ωγk = ∅, j 6= k, γ = c, f. j=1
Define Γ = ∂Ω and, for j = 1, . . . , Mγ : Dγ
Γj
= ∂Ωγj ∩ ΓD ,
Nγ
Γj
= ∂Ωγj ∩ ΓN ,
Γγjk = Γγkj = ∂Ωγj ∩ ∂Ωγk ,
γ = c, f.
3
Page 3 of 26
Figure 1: The fine Ωf (dashed lines) and coarse Ωc (solid lines) partitions of Ω along with the three spatial scales used in the definition of the new multiscale procedure.
us
cr
ip t
As indicated in Fig. 1, we refer to three length scales in the description of the new multiscale procedure: H, the mesh size for the coarse partition; h, the mesh size of an underlying fine grid and H, an intermediate length scale (h ≤ H ≤ H). The Multiscale Mixed Method gives an approximate solution to (2.2-2.3) in that the continuity of the pressure and normal component of the flux are imposed on the intermediate scale H (in terms of averages of fine scale Lagrange multipliers and normal components of fluxes, respectively) along the edges of the coarse partition Ωc . Flux conservation at the finest scale is obtained after a downscaling step. The new method is based on a iterative procedure similar to the one defined in [10] with a modified Robin interface condition. To formulate the hybridized mixed finite element method over Ωf , we start by defining the local finite dimensional spaces W (Ωfj ) ⊂ L2 (Ωfj ) and V (Ωfj ) ⊂ H(div; Ωfj ). The Neumann boundary N
for that we set rj ∈ P% (Γj f ) to satisfy Nf
Γj
N
= 0,
N
∀ϕ ∈ P% (Γj f ),
M
hub − rj , ϕi
an
value ub needs to be projected appropriately in each of the segments ∂Ωfj that intersect with ΓN , and
(2.4)
N
d
where P% (Γj f ) is the space of polynomial of fixed degree % in Γj f . Then we set
te
N Vrj (Ωfj ) = v ∈ V (Ωfj ) | v · ν = rj on Γj f . At this stage we introduce Lagrange multipliers `jk on the edges {Γfjk } (see for example[1]). We set
Ac ce p
Λ = `jk : `jk |Γf ∈ P% (Γfjk ) = Λjk , Γfjk 6= ∅ .
(2.5)
jk
We note that there are two values of P% assigned to Γfjk , i.e., Λjk and Λkj . The hybridized mixed finite element method for the pressure-velocity system (2.2-2.3) is given by seeking for each j = 1, 2..., Mf {pj , uj , `jk } ∈ W (Ωfj ) × Ve (Ωfj ) × Λ, such that ( div uj , w)Ωf = (q, w)Ωf , j
j
∀w ∈ W (Ωfj ),
(2.6)
(K −1 uj , v)Ωf − (pj , div v)Ωf +
X k6=j
j
j
h`jk , v · νiΓf + hpb , v · νiΓ jk
f D ∩∂Ωj
= 0, ∀v ∈ V0 (Ωfj ),
(2.7)
f where Ve (Ωfj ) = Vrj (Ωfj ) if ΓN 6= ∅ and Ve (Ωfj ) = V (Ωfj ) otherwise and νjf is the outward normal j
unit vector on Ωfj . We note that the above formulation is local in nature, i.e., one that is posed in every 4
Page 4 of 26
Ωfj . To transform it to a global algebraic system, the following continuity conditions are imposed: (2.8)
uj · νj + uk · νk = 0 on Γfjk ,
(2.9)
ip t
`jk = `kj on Γfjk ,
an
us
cr
where νj (and similarly for νk ) is the unit normal vector pointing outward of Ωfj . For a typical situation, imposing these conditions allow for the construction of a global linear equations governing the Lagrange multipliers (see for example, [1, 4]). In relative comparison to the traditional mixed finite element method, this hybridization technique has a clear advantage of reducing the system dimension. Another advantage is that the resulting global matrix from the hybridization is positive definite. Still, for problems such as multiphase flow in heterogeneous porous media, construction of a global matrix is prohibitive and may have to be avoided. In an attempt to provide an alternative along this direction, the authors of [10] developed an iterative procedure that is based on the local calculation (2.6) and (2.7) along with the imposition of a Robin transmission boundary conditions −χjk uj · νj + `jk = χjk uk · νk + `kj ,
on Γfjk ⊂ ∂Ωfj ,
(2.10)
−χjk uk · νk + `kj = χjk uj · νj + `jk ,
on Γfkj ⊂ ∂Ωfk ,
(2.11)
Ac ce p
2.2. Local problems on Ωcm
te
d
M
where χjk is a positive function on Γfjk . These Robin conditions are in fact equivalent to (2.8) and (2.9). In [10], the core of the iterative procedure was built by substituting `jk in (2.7) by the one in (2.10) with uk and `kj treated as known previous iteration values from the neighboring elements. The multiscale method we propose in this paper follows a closely related iterative procedure. The main distinction is that our method employs a red-black iterative procedure and an inclusion of information from intermediate scale through appropriate construction of local problems. This is the subject of the next subsection.
To describe the multiscale procedure for approximating the above formulation we construct a set of local problems on Ωcm , m = 1, · · · , Mc with Robin boundary conditions on ∂Ωcm . For each m = 1, · · · , Mc , we write the variational formulation of the pressure-velocity system (2.1) as to seek {pm , um } ∈ L2 (Ωcm ) × H(div; Ωcm ) such that ( div um , w)Ωcm = (q, w)Ωcm , ∀w ∈ L2 (Ωcm ),
(K
−1
m
(2.12)
m
(x)u , v)Ωcm − (p , div v)Ωcm
+ aum · ν m + g, v · ν m ∂Ωc = 0, ∀v ∈ H(div; Ωcm ), m
(2.13)
where we have imposed a Robin boundary condition on ∂Ωcm which is expressed in a generic format −aum · ν m + pm = g,
5
Page 5 of 26
for given functions a and g, and ν m the unit normal vector on ∂Ωcm pointing outward of Ωcm . We note that in the case Ωcm is adjacent to either of the global boundary, then we need to modify (2.13) appropriately to honor that particular global boundary.
ip t
Figure 2: Illustration of localization on Ωcm for the iterative procedure
us
cr
Fig. 2 outlines a schematic of the localization of (2.12) and (2.13) in Ωcm . Here, without loss of generality, we consider an internal Ωcm whose adjacent neighbors are Ωcn , where n = l, r, b, t. As described earlier, the key to the construction of the iterative procedure is the localization of (2.12) and (2.13) through a careful choice of Robin boundary condition on ∂Ωcm . In particular, the Robin boundary condition in (2.13) can be decomposed into X
m m au · ν + g, v · ν m ∂Ωc = aum · ν m + g, v · ν m Γc , (2.14) m
n=l,r,b,t
mn
an
where each Γcmn is as expressed in Fig. 2. Viewed this way, the natural choice for g in each term is a functional that depends on the neighboring element values, i.e., g c = g(un , pn ),
M
Γmn
d
where un and pn are velocity and trace pressure of Ωcn calculated at Γcmn . A detailed explanation of the modified Robin boundary condition is the subject of the next subsection. 2.3. The modified Robin condition
Ac ce p
te
To proceed, we set H and H such that r = H/H is an integer. As depicted in Fig. 2, a Γcmn is i , i = 1, · · · , r, each having length H. Thus decomposed into a collection of line segments Imn ∂Ωcm
=
[
Γcmn
=
n=l,r,b,t
[
r [
n=l,r,b,t
i=1
i Imn ,
i i = nψ = 4r. For each i, we define on Imn ⊂ Γcmn , which implies that dim Imn Z 1 i g i = g mn = (χimn un · ν n + pn )dl, Imn i H Imn
(2.15)
i . Hence, each term in the summation on where we assume that χmn is a positive constant on each Imn the right hand side of (2.14) can be expressed as
m m au · ν + g, v · ν m Γc
=
r X
mn
χimn um · ν m + g imn , v · ν m
i Imn
.
(2.16)
i=1
6
Page 6 of 26
We note that imposing this condition is equivalent to f m m i i −χimn um j · ν + pjk = g mn , ∀ Γjk ⊂ Imn ,
(2.17)
f m m where Γfjk = ∂Ωfjm ∩ ∂Ωfkn (see Fig. 2), um j = u Ωf , and pjk is the trace pressure at Γjk .
ip t
jm
As for the numerical approximation of (2.12-2.13), we use the hybridized mixed finite element approximation that is akin to (2.6) and (2.7).
cr
2.4. The MuMM iteration
M
an
us
We now define an iterative procedure that localizes the calculations and gives an approximate solution to the global problem (2.6-2.7). Consider the domain Ω and label its coarse partitions as “red” and “black” in an ordinary checkerboard order. The iterative process is listed in Algorithm 1 where we have assumed to use the hybridized mixed finite element approximation for solving the local problem (2.12) and (2.13). In essence, either red or black calculation in Algorithm 1 is represented by the schematic in Fig. 2. For example, in the red calculation, Ωcm is the red subdomain, while Ωcn , n = l, r, b, t are the black subdomains from which the boundary information are gathered. Similar setting is also true for black calculation. We remark that in the case of H = H = h and Ωf = Ωc the convergence of this iterative scheme has been established in [10]. Algorithm 1 MuMM Iteration
te
d
Set initial guess for u(0) , `(0) . for α = 1 until convergence do Red Subdomain: • For each red element, compute p(α) , u(α) , `(α) by solving (2.12-2.13)
Ac ce p
• Use the (α − 1) iterate from the black subdomains nearest neighbors in the right hand side of (2.14) Black Subdomain: • For each black element, compute p(α) , u(α) , `(α) by solving (2.12-2.13) • Use α iterate from the red subdomains nearest neighbors in the right hand side of (2.14) end for
3. Computational Implementation
In our numerical experiments we shall adopt the lowest order Raviart-Thomas elements and we consider this space in the discussion of the new scheme. The natural degrees of freedom on the element Ωfj , where the elements are squares of side length h, are the (constant) value pj of the pressure over Ωfj , which can be interpreted as the value at the center of the element, and the four (constant) values ujβ , β = l, r, b, t, of the outward normal component of the total flux across the edges of the element.
7
Page 7 of 26
us
cr
ip t
The Lagrange multipliers `jk and `kj on Γfjk are also constants. We shall also assume that the absolute permeability K(x) and the external flow rate q to be constant on each element. However, it is obvious that Algorithm 1 can easily become an expensive calculation especially when the dimension of the local calculation in Ωcm increases significantly. This is even more pronounced when it is applied to multiphase flow simulations that require repeated calculations of pressure and velocity and hence repeated algorithm calls. For that matter, it is imperative to devise an innovative way to be incorporated in Algorithm 1 that allows for efficient calculation and yet maintains an acceptable accuracy. Our strategy in this direction is to construct a set of multiscale basis functions associated with Ωcm . Calculation of these basis functions is a one time preprocessing step and thus is reasonable to expend. Furthermore, with these basis functions in place, the iteration procedure of Algorithm 1 is only conducted on ∂Ωcm . In the next subsection, we describe the aforementioned basis functions.
M
an
3.1. Mixed Multiscale Basis Functions of Robin Type As mentioned, the multiscale basis functions are associated with Ωcm , m = 1, ·i· · , Mc . Moreover, m the number of multiscale basis functions for each Ωc depends on how many Imn are employed. Thus the total number of basis functions for each Ωm c is nψ = 4r, where as described earlier, r = H/H. This means there are Mc nψ of multiscale basis functions constructed on the whole domain Ω. As in Algorithm 1, the local problem on Ωm c for the multiscale basis functions areclosed with a i Robin boundary condition. One Imn is tied to one basis function, which we denote by ψ imn , φimn . These functions satisfy variational formulation that is similar to (2.12) and (2.13) (with q = 0). Furthermore, the Robin boundary condition is set to satisfy (3.1)
is as depicted in Fig. 3, namely
te
where
i gmn
d
i −χimn ψ imn · ν m + φimn = gmn in ∂Ωcm ,
( i 1, in Imn = i . 0, in ∂Ωcm \Imn
Ac ce p
i gmn
Taking all these into account we arrive at the following variational formulation: find the divergence i,H = ψ imn , φimn satisfying free Mixed Multiscale Basis Functions (MMBFs) Ψmn ( div ψ imn , w)Ωcm = 0, ∀w ∈ L2 (Ωcm ), (K
(3.2)
−1
(x)ψ imn , v)Ωcm − (φimn , div v)Ωcm
i + χimn ψ imn · ν m + gmn , v · ν m ∂Ωc m
= 0, ∀v ∈ H(div; Ωcm ).
(3.3)
We note that as in Algorithm 1, this variational formulation can be approximated using the hybridized mixed finite element. If x ∈ ΓD or x ∈ ΓN this construction is modified to take trivial Dirichlet or Neumann conditions. Moreover, to the set of basis functions associated with a Ωcm that is adjacent to the external boundary, an extra MMBF has to be computed, with the prescribed external boundary condition and zero value for the Robin condition in the remaining boundaries. We remark that in [14] a similar set of basis functions was defined with distinct boundary conditions. 8
Page 8 of 26
i H Figure 3: One line segment Imn where the Robin boundary condition for ψmi is imposed.
S Ωcm =
X
r X
Aimn Ψi,H mn ,
(3.4)
cr
n=l,r,b,t i=1
ip t
3.2. The MuMM iteration revisited In performing the MuMM iteration, the multiscale basis functions described above have to be pre-computed. Within the iteration updating the local solution for each Ωcm consists of computing
an
us
where S Ωcm indicates the approximation of um , pm . The coefficients Aimn in (3.4) are given by the (previous iteration) nearest neighbor modified Robin condition (see (2.15)); there is no linear algebra problem to be solved in each iteration. Note that in a computational implementation of (3.4) the sum should be evaluated only at the boundaries of each Ωcm . After convergence of the MuMM iteration a downscaling step to be discussed next will produce the solution inside Ωcm . We remark that:
M
• The approximate solution for the global problem (2.2-2.3) is then given by Mc X
S=
S Ωcm .
(3.5)
m=1
d
• Flux conservation is maintained in the H scale.
Ac ce p
te
• The balance between numerical accuracy and numerical efficiency is determined by the choice i,h of span{Ψi,H mn } ⊂ span{Ψmn }. The closer H is to h the larger is the number of MMBFs to be computed and the closer the MuMM solution is to the fine grid solution. 3.3. Downscaling As mentioned above, after convergence the flux conservation is maintained in the H scale. In order to have the flux conservation in the h scale, we employ a two-step procedure. We first solve local Neumann problems to recover the solution inside each element of Ωc . The boundary condition for these problems specified in the h scale is set to be the average flux produced by the MuMM iteration at the boundaries between two neighbor elements of Ωc . Next we solve a family of local Neumann problems in a staggered grid. Fig. 4 illustrates this procedure. After this step we recover an approximate solution in the fine scale. The computational cost of the downscaling is equivalent to the contruction of one MMBF, in both cases we solve one elliptic problem for each subdomain. R EMARK : Although the procedure has been defined for uniform structured grids, we believe that it could be implemented on more complicated (non-uniform or unstructured) grids. Such implementations are, however, outside the scope of this work.
9
Page 9 of 26
Figure 4: The solid boundary corresponds to a cell of Ωc and indicates the region where a local Neumann problem is solved in step-1 of the downscaling procedure. In the second step of the procedure the local solution is computed in the region indicated by a dashed boundary.
ip t
4. Numerical Results
φ
an
us
cr
In order to show that MuMM is a viable alternative to solve reservoir simulation problems, we present numerical results for the simulation of two-dimensional fluid flows associated with two cases, namely, quarter of a 5-spot problem and a slab geometry on a channelized field. In all simulations, we compare the MuMM with the usual traditional simulation on the mesh in which the fine scale permeability field is resolved. Aspects of comparison include the accuracy of the approximate velocity, efficiency of MuMM in terms of number of iterations to reach a convergence tolerance, and the robustness of MuMM when incorporated in a single phase flow simulation. With respect to this single phase flow simulation, the model governs the physical transport of fluids that uses the approximate velocity as an input to the equation ∂c + u · ∇c = 0. ∂t
(4.1)
te
d
M
This equation is solved using the genuinely multidimensional second-order central scheme for conservation laws described in described in [18]. Here φ is the porosity and c is the concentration which is defined on the fine mesh. We consider a constant c(x, t = 0) for initial condition. In practice, it is of interest to gather the concentration at a dimensionless time PVI, whose value is Z t Z PVI = Vp −1 u · ν dl dτ, (4.2) 0
∂Ωout
Ac ce p
for a given dimensional time t. Here Vp is the total pore-volume of the reservoir. 4.1. Example 1: The quarter of a 5-spot problem This is a standard problem in reservoir simulation and thus it is an important example to evaluate the performance of the MuMM. In this simulation the global fine grid is of 128 × 128 elements in the region (0, 1) × (0, 1). The global boundary condition is of Neumann type (with no flow on the boundary), and there is an injection through a well at the lower left corner of the computational domain, and a production, at the same rate, through a well in the upper right corner. The experiment uses a 2-D heterogeneous permeability whose log value field ξ(x) is shown in Fig. 5. We set k(x) = exp(δ ξ(x)), where δ is chosen in order to set a high permeability ratio kmax /kmin = 106 . The field ξ(x) is a self similar (Gaussian), (see [2, 3]), assumed to be stationary and isotropic, and is characterized by its (zero) mean and covariance function C(x, y) given by: C(x, y) = |x − y|−0.5 .
(4.3)
10
Page 10 of 26
Figure 5: Log permeability field used in the quarter of a 5-spot problem.
ip t
Figure 6: Comparison of concentration surface profiles in Example 1, where the velocity is from direct fine grid simulation (left) and from MuMM (right). The time in PVI is, from top to bottom, 0.32 and 0.65.
k uMuMM − ufine k j j m
M
an
us
cr
In Fig. 6 we show concentration profiles at two different PVIs, in which the velocity used in the left plots is from the direct fine grid calculation using the iterative procedure from [10], while its counterpart in the right plots is from the MuMM with coarse grid of 4 × 4 elements (i.e., H = 1/4, and thus making each coarse element Ωcm to possess 32 × 32 fine elements). We use H = 1/16, which gives dim{Ψi,H mn } = 4H/H = 4 · 4 = 16, i.e., there are 16 MMBFs per coarse element. An additional basis function has to be constructed for lower left and upper right coarse element in which the injection and production well is respectively located. Each of these extra basis functions is defined with trivial Robin conditions that takes the source/sink terms into account. The comparison in Fig. 6 indicates that the MuMM produces a good quality numerical approximation of the solution in the fine grid. As described, the aforementioned simulations take the approximate velocity on the fine grid as an input, so maintaining a satisfactory level of its accuracy is imperative. Further investigation of the relative error spatial profile for MuMM approximation depicted in Fig. 7 confirms the previous finding. The relative error for each fine element Ωfj , j = 1, ..., Mf is evaluated as , where m = max k ufine k, j
(4.4)
d
j
Ac ce p
te
where k . k is the usual Euclidean norm. Simulation to create this figure uses a fixed H = 1/16 and varying H = 1/16, 1/8, 1/4, respectively having 8 × 8, 16 × 16, and 32 × 32 fine elements in each coarse element. Consequently, we maintain a fixed total number of 64 MMBFs on the whole domain. We note that the maximum relative error in each of the plots from left to right is respectively 0.15, 0.08, and 0.019. This result suggests that when the total number of MMBFs remains fixed on the whole domain, the accuracy is better when the H scale is coarser but at the expense of having a larger degrees of freedom in each multiscale basis function calculation. To complete the study of the MuMM for the quarter of 5-spot problem, we present in Table 1 the number of iterations of the method for various ratios of permeability along with two choices of coarse grids. Note that the number of iterations increases slowly as a function of the strength of the permeability field. The criterion to stop iteration is based on how close the flux conservation is satisfied in the H scale.
Figure 7: The relative errors of velocity obtained from MuMM in Example 1 for various H: 1/16 (left), 1/8 (middle), 1/4 (right). The value of H is fixed at 1/16.
11
Page 11 of 26
103 104 105 106
Coarse grid of 8×8 56 63 74 91
Coarse grid of 4×4 24 26 29 34
cr
kmax /kmin
ip t
Table 1: Number of iterations in MuMM for various ratios of permeability and coarse grid settings in Example 1.
us
4.2. Example 2: Slab geometry on a channelized field
Ac ce p
te
d
M
an
The complex features of a channelized permeability field makes this example very challenging for testing multiscale methods, so that it is important to consider such a case in the evaluation of the MuMM performance. We use Ω = (0, 3 23 ) × (0, 1) which is discretized into 220 × 60 fine elements on which the permeability data is posed. The MuMM employs a coarse grid of 11 × 3 implying H = 1/3. This way, each coarse element possesses 20 × 20 fine elements. The left and right global boundary conditions are of Dirichlet type, and the top and bottom global boundary conditions are of homogeneous Neumann type (no flow). Fig. 8 shows the spatial profile ξ(x), which is the log permeability of layer 36 of SPE10 benchmark1 . The permeability is given by k(x) = exp(δ ξ(x)), where we adjust the parameter δ to give us permeability ratios of interest to our study. In the first case we consider the permeability ratio of 176 and three choices for the H scale, namely, 1/3, 1/6, and 1/12, respectively giving 4, 8 and 16 MMBFs in each coarse element. In Fig. 9 we compare tracer injection profiles for these three choices of H with the fine grid solution (computed using the iterative procedure from [10]). We note that with only 16 MMBFs per coarse element, the MuMM iteration can recover a very similar profile to the one obtained from a direct fine grid simulation. Furthermore, Fig. 10 compares the tracer cut curves, which is the fraction of the tracer in the produced fluid represented as Z cu · ν dl ∂Ωout F (t) = 1 − Z u · ν dl ∂Ωout
where ∂Ωout denotes outflow boundary. The tracer cut is plotted against the PVI given in (4.2). Next, to illustrate the importance of downscaling step within MuMM, we use 16 MMBFs per coarse element (H = H/4) and kmax /kmin = 104 . We calculate the velocity using MuMM with and without the downscaling step and make comparison against the direct fine grid simulation whose profile is plotted in Fig. 11. We note on the dominance of the flow in the regions that exhibit the channel 1
More information at http://www.spe.org/csp
12
Page 12 of 26
Figure 8: The log permeability field from the SPE10 model (layer 36)
ip t
Figure 9: Tracer concentrations of (from top to bottom): MuMM for H=H , H=H/2, H=H/4 and the fine grid, in Example 2. Here H = 1/3 and time (in PVI) is 0.77.
an
us
cr
feature. In particular, we zoom on the region confined by the white lines (see Fig. 11) and compare the respective velocities in Fig. 12. As can be seen from that figure, the downscaling procedure improves the numerical solution. Finally we conclude the analysis of MuMM for Example 2 with Table 2, which shows the iterations of MuMM for various permeability ratios. The cost per iteration is the cost of computing (3.4) for each subdomain. These results indicate the efficiency of the MuMM to capture the features of the channelized permeability field. As in Example 1, the criterion to stop iteration is based on how close the flux conservation is satisfied. R EMARK : In line with [9], in expressions (2.10 -2.11) we set
2Kj Kk ch H , where Kjk = , H Kj + Kk Kjk
M
χjk =
and c is a dimensionless constant. Within the modified Robin condition (2.15) we set
te
d
χimn =
ch i
K jk
i
Ac ce p
H set in the line for the log normal permeability field, where K jk is the average value among Kjk i . For the channelized field we set χi segment Imn mn to be the maximum value among χjk previously i . determined in the line segment Imn
5. Concluding Remarks
In this paper, we have laid out a framework for efficient numerical approximation of the velocity in the presence of heterogeneous media using a mixed method that is incorporated into a nonoverlapping domain decomposition procedure. The main attraction of the proposed framework is its immediate parallelizability owing to the construction of the multiscale basis functions within each coarse elements. On the global level, an iterative procedure is devised which only involves calculation on the edges of the coarse elements. This iteration avoids the construction and solving of linear algebra.
Figure 10: Tracer cut curve profiles for MuMM and the fine grid from Example 2.
13
Page 13 of 26
kmax /kmin 102 103 104
cr
MuMM Iterations 41 42 88
ip t
Table 2: Number of iterations using different ratios of permeability in Example 2:
us
Figure 11: The velocity field with the log permeability on the background. See Fig. 12 for comparison of MuMM with fine grid simulations on the zoomed region (a box confined with the white lines).
an
After the downscaling step we recover a good approximation to the solution in the fine scale, and mass conservation holds also in that scale. Sources and sinks are naturally incorporated in the procedure, by considering an extra MMBF (constructed with trivial Robin conditions) that takes into account the source or sink in the coarse partition element where it is located. The authors are currently implementing the 3D version of the procedure described here in a CPUGPU cluster. Moreover, the numerical analysis of the MuMM is underway.
M
Acknowledgments
Ac ce p
te
d
The authors wish to thank Jim Douglas, Jr. for some discussions about the method presented in the this paper and Dr. Marcos Mendes for his help with some numerical simulations. V. Ginting and F. Pereira are supported in part by DOE DE-FE0004832 and by grants from the Center for Fundamentals of Subsurface Flow of the School of Energy Resources of the University of Wyoming, under contract numbers WYDEQ49811GNTG and WYDEQ49811PER (Ginting) and WYDEQ49811GNTG, WYDEQ49811PER and WYDEQ49811FRTD (Pereira). Ginting is also supported in part by grants from US Department of Energy (DE-SC0004982) and the National Science Foundation (DMS-1016283). Pereira is also supported in part by the Wyoming State Legislature, account number 1100 20352 2012. A. Francisco was supported by grant from CAPES/FULBRIGHT Program (BEX 2984/09-1). J. Rigelo is supported by DOE DE-FE0004832. References
[1] D. N. Arnold, F. Brezzi, Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, RAIRO Mod´el. Math. Anal. Num´er. 19 (1) (1985) 7–32.
Figure 12: A zoomed window selected from the domain in Fig. 11 above (from left to right): fine grid, MuMM iteration before and after downscaling. Note the importance of the downscaling procedure to correct the flux within the channel.
14
Page 14 of 26
[2] M. R. Borges, F. Furtado, F. Pereira, H. P. Amaral Souto, Scaling analysis for the tracer flow problem in self-similar permeability fields, Multiscale Model. Simul. 7 (3) (2008) 1130–1147.
ip t
[3] M. R. Borges, F. Pereira, H. P. Amaral Souto, Efficient generation of multi-scale random fields: a hierarchical approach, Int. J. Numer. Methods Biomed. Eng. 26 (2) (2010) 176–189.
cr
[4] G. Chavent, J. Roberts, A unified physical presentation of mixed, mixed-hybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems, Adv. Water Resour. 14 (6) (1991) 329–348.
us
[5] W. Chen, M. Gunzburger, F. Hua, X. Wang, A parallel Robin-Robin domain decomposition method for the Stokes-Darcy system, SIAM J. Numer. Anal. 49 (3) (2011) 1064–1084.
an
[6] Z. Chen, G. Huan, Y. Ma, Computational methods for multiphase flows in porous media, Computational Science & Engineering, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2006. [7] J. Chu, Y. Efendiev, V. Ginting, T. Hou, Flow based oversampling technique for multiscale finite element methods, Adv. Water Resour. 31 (4) (2008) 599–608.
d
M
[8] L. C. Cowsar, M. F. Wheeler, Parallel domain decomposition method for mixed finite elements for elliptic partial differential equations, in: Fourth International Symposium on Domain Decomposition Methods for Partial Differential Equations (Moscow, 1990), SIAM, Philadelphia, PA, 1991, pp. 358–372.
te
[9] J. Douglas, Jr., F. Furtado, F. Pereira, On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs, Comput. Geosci. 1 (2) (1997) 155–190.
Ac ce p
[10] J. Douglas, Jr., P. J. Paes-Leme, J. E. Roberts, J. P. Wang, A parallel iterative procedure applicable to the approximate solution of second order partial differential equations by mixed finite element methods, Numer. Math. 65 (1) (1993) 95–108. [11] Y. Efendiev, T. Y. Hou, Multiscale finite element methods, vol. 4 of Surveys and Tutorials in the Applied Mathematical Sciences, Springer, New York, 2009, theory and applications. [12] R. E. Ewing, J. Wang, Analysis of the Schwarz algorithm for mixed finite elements methods, RAIRO Mod´el. Math. Anal. Num´er. 26 (6) (1992) 739–756. [13] R. E. Ewing, J. Wang, Analysis of multilevel decomposition iterative methods for mixed finite element methods, RAIRO Mod´el. Math. Anal. Num´er. 28 (4) (1994) 377–398. [14] B. Ganis, I. Yotov, Implementation of a mortar mixed finite element method using a multiscale flux basis, Comput. Methods Appl. Mech. Engrg. 198 (49-52) (2009) 3989–3998.
15
Page 15 of 26
[15] R. Glowinski, W. Kinton, M. F. Wheeler, Acceleration of domain decomposition algorithms for mixed finite elements by multi-level methods, in: Third International Symposium on Domain Decomposition Methods for Partial Differential Equations (Houston, TX, 1989), SIAM, Philadelphia, PA, 1990, pp. 263–289.
ip t
[16] R. Glowinski, M. F. Wheeler, Domain decomposition and mixed finite element methods for elliptic problems, in: First International Symposium on Domain Decomposition Methods for Partial Differential Equations (Paris, 1987), SIAM, Philadelphia, PA, 1988, pp. 144–172.
cr
[17] V. Kippe, J. E. Aarnes, K.-A. Lie, A comparison of multiscale methods for elliptic problems in porous media flow, Comput. Geosci. 12 (3) (2008) 377–398.
us
[18] F. Pereira, A. Rahunanthan, A semi-discrete central scheme for the approximation of two-phase flows in three space dimensions, Math. Comput. Simulation 81 (10) (2011) 2296–2306.
Ac ce p
te
d
M
an
[19] M. F. Wheeler, T. Wildey, G. Xue, Efficient algorithms for multiscale modeling in porous media, Numer. Linear Algebra Appl. 17 (5) (2010) 771–785.
16
Page 16 of 26
d
te
Ac ce p us
an
M
cr
H
h
H
Page 17 of 26
ip t
1
ip t
Imn 2
Ac ce p
te
d
M
an
us
cr
Imn
Page 18 of 26
ip t cr us an M d te Ac ce p
0
0
1
Imn 1
0
0
0 0
0
1 Figure 3: One line segment Imn where the Robin boundary condition for ψ
Page 19 of 26
Page 20 of 26
d
te
Ac ce p us
an
M
cr
ip t
ip t cr us an
120
M
100
80
te
40
d
60
20
40
60
80
100
120
Ac ce p
20
Figure 5: Permeability field used in example 1 with permeability ratio of kmax /kmin = 106 . The red and blue colors refer to the highest and lowest permeability values, respectively .
13
Page 21 of 26
ip t cr us an M d te
Ac ce p
Figure 6: Concentration surface profiles of MuMM (left side) and the fine grid (right side), in Example 1. Note the good agreement between the solutions for a large value of kmax /kmin = 106
14
Page 22 of 26
−3
x 10 18 120
0.02
120 16
ip t
0.018 100
100
14
0.016 0.014
80
12
80
0.012 0.01
10
60
8
cr
60
0.008 40
40 0.006 0.004
20
20
40
60
80
100
120
20
40
60
80
100
4 2
120
Ac ce p
te
d
M
an
20
us
0.002
6
Page 23 of 26
ip t cr us an M d te Ac ce p Figure 8: The permeability field from the SPE10 model (layer 36) with permeability ratio of 176 used in the second example.
16
Page 24 of 26
ip t cr us an M d te Ac ce p
Figure 9: Tracer concentrations of: MuMM for H=H (upper left), H=H/2 (upper right), H=H/4 (bottom left) and the fine grid (bottom right), in Example 2.
Figure 10: Tracer cut curve profiles for MuMM and the fine grid from Example 2. A good agreement is observed for H ≤ H/2.
17
Page 25 of 26
ip t cr us an M d te Ac ce p
Figure 11: The permeability field together with the velocity field for (from top to bottom): fine grid, MuMM iteration before and after downscaling.
Figure 12: A window selected from the domain in Fig. 11 above (from left to right): fine grid, MuMM iteration before and after downscaling. Note the importance of the downscaling procedure to correct the flux within the channel.
18
Page 26 of 26