Stochastic multiscale modeling of subsurface and surface flows. Part - - PowerPoint PPT Presentation

Stochastic multiscale modeling of subsurface and surface flows. Part I: Multiscale mortar mixed finite elements for Darcy flow

Ivan Yotov

Department of Mathematics, University of Pittsburgh

KAUST WEP Workshop January 30–February 1, 2010

Joint work with Todd Arbogast, Gergina Pencheva, Sunil Thomas, and Mary F.

Wheeler, The University of Texas at Austin; Benjamin Ganis, University of

Pittsburgh

Department of Mathematics, University of Pittsburgh 1

Energy and environment

• Ground water and surface water contamination
• Hydrocarbon energy production

Department of Mathematics, University of Pittsburgh 2

Reservoir rock

Department of Mathematics, University of Pittsburgh 3

Mathematical and numerical challenges for modeling porous media

• Multiscale - space and time
• Multiphysics - aquifer, surface water, waterflood, CO2, polymer, geomechanics
• Multiphase - gas, aqueous, multiple flowing phases
• Highly nonlinear coupled PDE systems - advection, reaction, diffusion/dispersion,

capillary effects

• Complex geology and geometry - faults, fractures, layers
• Very large scale computations

– millions of unknowns – parallel computing

Department of Mathematics, University of Pittsburgh 4

Multiblock approach for multiphysics problems

Fault Water flood Well Fully implicit model Time splitting scheme CO

2 flood

• complex geometry and gridding
• multinumerics
• multiscale resolution
• parallel algorithms

Department of Mathematics, University of Pittsburgh 5

X Y Z

Groundwater flow in a faulted aquifer Numerical grids and low conductivity barriers

Department of Mathematics, University of Pittsburgh 6

X Y Z

Contaminated groundwater flow DNAPL concentration at 44 days

Department of Mathematics, University of Pittsburgh 7

Outline

• Background and motivation
• A multiscale mortar mixed finite element method
• A priori error estimates
• A domain decomposition algorithm
• A posteriori error estimates
• Mortar and subdomain adaptivity
• A relationship between multiscale mortar MFE methods and subgrid upscaling

methods

• A multiscale flux basis formulation
• Extension to two-phase flow

Department of Mathematics, University of Pittsburgh 8

Single phase flow model

u = −K∇p in Ω ⊂ Rd (d = 2, 3) (Darcy’s law) ∇ · u = f in Ω (conservation of mass) u · n = 0

• n ∂Ω

(no flow BC) Variational mixed formulation H(div; Ω) = {v : v ∈ (L2(Ω))d, ∇ · v ∈ L2(Ω)} V = {v ∈ H(div; Ω) : v · n = 0 on ∂Ω} W = L2

0(Ω) = {w ∈ L2(Ω) :

• Ω

w dx = 0}. Find u ∈ V, p ∈ W such that (K−1u, v) = (p, ∇ · v), v ∈ V, (∇ · u, w) = (f, w), w ∈ W.

Department of Mathematics, University of Pittsburgh 9

The mixed finite element method

Th - finite element partition Vh × Wh ⊂ V × W - mixed finite element spaces

• ✁

pressure velocity

Find uh ∈ Vh, ph ∈ Wh such that (K−1uh, v) = (ph, ∇ · v), v ∈ Vh, (∇ · uh, w) = (f, w), w ∈ Wh.

• Simultaneous (accurate) approximation of pressure and velocity
• Local mass conservation: for each element E,

w =

• 1 on E,

0 otherwise = ⇒

• E

∇ · uh =

• E

q.

• Continuity of normal flux across element faces: for each e = ∂E1 ∩ ∂E2,

uh|E1 · ne = uh|E2 · ne.

Department of Mathematics, University of Pittsburgh 10

Error estimates and convergence testing

Theorem[Ingram-Wheeler-Xue-Y.]: u − uh + ∇ · (u − uh) + p − ph ≤ Ch

pres 2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2

Computed pressure and velocity

errp 0.007 0.0065 0.006 0.0055 0.005 0.0045 0.004 0.0035 0.003 0.0025 0.002 0.0015 0.001 0.0005

Pressure and velocity error

1/h p − ph Rh

p

u − uh Rh

u

∇ · (u − uh) Rh

∇·u

8 0.120E0 0.164E0 0.188E0 16 0.605E-1 1.0 0.834E-1 1.0 0.941E-1 1.0 32 0.304E-1 1.0 0.417E-1 1.0 0.470E-1 1.0 64 0.152E-1 1.0 0.208E-1 1.0 0.235E-1 1.0

Department of Mathematics, University of Pittsburgh 11

Motivation for multiscale modeling: flow in heterogeneous porous media

Heterogeneous permeability varies on a fine scale. Full fine scale grid resolution ⇒ large, highly coupled system of equations ⇒ solution is computationally intractable

• Variational Multiscale Method

– Hughes et al; Brezzi – Mixed FEM: Arbogast et al

• Multiscale Finite Elements

– Hou, Wu, Cai, Efendiev et al – Mixed FEM: Chen and Hou; Aarnes et al New approach: based on domain decomposition and mortar finite elements More flexible - easy to improve global accuracy by refining the local mortar grid where needed

Department of Mathematics, University of Pittsburgh 12

Multiscale finite element/subgrid upscaling methods

Lǫu = f ⇒ u ∈ V : a(u, v) = (f, v) ∀ v ∈ V. Multiscale approximation: H - coarse grid, h ≈ ǫ - fine grid (subgrid) VH,h = VH + V ′

h

Basis for V ′

h(E): φE h,i, i = 1, . . . , NE,

aE(φE

H,i + φE h,i, vh) = 0

∀ vh ∈ Vh(E)

H

• h

Multiscale solution: uH,h ∈ VH,h, a(uH,h, vH,h) = (f, vH,h) ∀ vH,h ∈ VH,h

Department of Mathematics, University of Pittsburgh 13

Multiblock formulation for single phase flow

¯ Ω = ∪n

i=1¯

Ωi; Γij = ∂Ωi ∩ ∂Ωj On each block Ωi: u = −K∇p in Ωi ∇ · u = q in Ωi u · n = 0

• n ∂Ωi ∩ ∂Ω

On each interface Γij: pi = pj

• n Γij

[u · n]ij = 0

• n Γij

where pi = p|∂Ωi [u · n]ij ≡ u|Ωi · n − u|Ωj · n

Department of Mathematics, University of Pittsburgh 14

Multiblock discretization spaces

Vh =

n

• i=1

Vh,i, Wh =

n

• i=1

Wh,i, Mh =

• 0≤i<j≤n

Mh,ij λh|Γij ∈ Mh,ij,

• Γij

[uh · n]ijµ = 0, µ ∈ Mh,ij. Subdomain grids do not need to match.

Department of Mathematics, University of Pittsburgh 15

The mortar mixed finite element method

Find uh ∈ Vh, ph ∈ Wh, λh ∈ Mh s.t. for 1 ≤ i ≤ n (K−1uh, v)Ωi − (ph, ∇ · v)Ωi + λh, v · niΓi = 0, v ∈ Vh,i, (∇ · uh, w)Ωi = (q, w)Ωi, w ∈ Wh,i,

n

• i=1

uh · ni, µΓi = 0, µ ∈ Mh. Stability, optimal convergence, superconvergence: Y.(1996,97), Arbogast, Cowsar, Wheeler, Y. (2000)

Department of Mathematics, University of Pittsburgh 16

Two-scale formulation: mortar upscaling

Two-scale problem:

H

• mortar dof

h

• Each block is an element of the coarse grid.
• Each block is discretized on the fine scale.
• A coarse mortar space on each interface.
• Result: Effective solution, fine scale on subdomains, coarse scale flux matching

Department of Mathematics, University of Pittsburgh 17

Multiscale mortar mixed finite element method

Allow for different scales and polynomial approximations on interfaces and subdo- mains. Assume Pk ⊂ Vh,i, Pl ⊂ Wh,i, Pm ⊂ MH, m ≥ k + 1 Find uh ∈ Vh, ph ∈ Wh, λH ∈ MH s.t. for 1 ≤ i ≤ n (K−1uh, v)Ωi − (ph, ∇ · v)Ωi + λH, v · niΓi = 0, v ∈ Vh,i, (∇ · uh, w)Ωi = (q, w)Ωi, w ∈ Wh,i,

n

• i=1

uh · ni, µΓi = 0, µ ∈ MH. Stability assumption: µ0,Γi,j ≤ C(Qh,iµ0,Γi,j + Qh,jµ0,Γi,j), µ ∈ MH, 1 ≤ i < j ≤ n.

Department of Mathematics, University of Pittsburgh 18

An approximation result

Weakly continuous velocities: Vh,0 =

• v ∈ Vh :

n

• i=1

v|Ωi · ni, µΓi = 0 ∀ µ ∈ MH

• .

Equivalent formulation: find uh ∈ Vh,0 and ph ∈ Wh such that (K−1uh, v) − (ph, ∇ · v) = 0, v ∈ Vh,0, (∇ · uh, w) = (q, w), w ∈ Wh Interpolation operator Π0 : V → Vh,0 such that (∇ · (Π0q − q), w)Ω = 0, w ∈ Wh. Π0q − q0 ≤ C

n

• i=1

(qr,Ωihr + qr+1/2,ΩihrH1/2), 1 ≤ r ≤ k + 1

Department of Mathematics, University of Pittsburgh 19

A priori error estimates

Theorem: u − uh ≤ C(Hm+1/2 + hk+1), ∇ · (u − uh) ≤ Chl+1 |||u − uh||| ≤ C(Hm+1/2 + hk+1H1/2) p − ph ≤ C(Hm+3/2 + hk+1H + hl+1) |||p − ph||| ≤ CHu − uhH(div) Balance u − uh or |||p − ph||| error terms (with l = k): H = h

k+1 m+1/2 ⇒ u − uh ≤ Chk+1,

|||p − ph||| ≤ Ch

k+1+

k+1 m+1/2

For example, RT0, k = 0, and quadratic mortars, m = 2, H = h2/5 : u − uh ≤ Ch, p − ph ≤ Ch, |||p − ph||| ≤ Ch1+2/5

Department of Mathematics, University of Pittsburgh 20

Parallel domain decomposition

Two types of subdomain problems:

BC BC q 1 q2 BC BC BC BC

gH(µ) =

n

• i=1

¯ uh,i · ni, µΓi

λ

aH,i(λ, µ) = −u∗

h,i(λ) · ni, µΓi

aH(λ, µ) =

n

• i=1

aH,i(λ, µ) The solution (uh, ph, λH) to the original problem satisfies AHλH = gH

• r

aH(λH, µ) = g(µ), ∀µ ∈ MH, with uh = u∗

h(λH) + ¯

uh, ph = p∗

h(λH) + ¯

ph.

Department of Mathematics, University of Pittsburgh 21

Interface iteration

Lemma The interface operator AH : MH → MH is symmetric and positive semi- definite. aH,i(λ, µ) = (K−1u∗

h,i(λ), u∗ h,i(µ)).

AH : λH → [u∗

h(λH) · n] is a Steklov-Poincare operator.

Apply the Conjugate Gradient method for AHλH = gH. Computing the action of the operator (needed at each CG iteration):

• Given mortar data λH ∈ MH, project onto subdomain grids:

λH → QhiλH

• Solve local problems in parallel with boundary data QhiλH
• Project fluxes onto the mortar space and compute the jump:

uh,i · ni → QT

hiuh,i · ni,

AHλH = QT

h1uh,1 · n1 + QT h2uh,2 · n2

Department of Mathematics, University of Pittsburgh 22

Numerical experiments

m H ||p − ph|| ||u − uh|| |||p − ph||| |||u − uh||| |||p − λH||| full K diag K 2 h1/2 1 1 1.5 1.25 1.25 1.5 1 2h 1 1 2 1.5 1.5 2 Table 1: Theoretical convergence rates for quadratic and linear mortars. Example 1: p(x, y) = x3y4 + x2 + sin(xy)cos(y), K =

• (x + 1)2 + y2

sin(xy) sin(xy) (x + 1)2

• .

Example 2: p(x, y) = x2y3 + cos(xy) 2x+9

20

2 y3 + cos 2x+9

20 y

, K =

• I,

0 ≤ x ≤ 1/2, 10 ∗ I, 1/2 ≤ x ≤ 1

Department of Mathematics, University of Pittsburgh 23

Computed solution for Example 1

• A. Discontinuous quadratic mortars
• B. Discontinuous linear mortars

Computed pressure (shade) and velocity (arrows).

Department of Mathematics, University of Pittsburgh 24

Convergence rates for Example 1

1/h iter. ||p − ph|| ||u − uh|| |||p − ph||| |||u − uh||| |||p − λH||| 4 8 2.64E-1 2.03E-1 4.62E-2 2.13E-2 4.45E-2 16 13 6.37E-2 4.86E-2 2.83E-3 1.81E-3 2.72E-3 64 15 1.59E-2 1.21E-2 1.75E-4 1.60E-4 1.69E-4 256 16 3.98E-3 3.03E-3 1.09E-5 1.77E-5 1.08E-5 rate O(h1.01) O(h1.01) O(h2.01) O(h1.71) O(h2.00) Continuous quadratic mortars on non-matching grids 1/h iter. ||p − ph|| ||u − uh|| |||p − ph||| |||u − uh||| |||p − λH||| 4 4 2.63E-1 2.04E-1 4.54E-2 2.35E-2 4.55E-2 8 7 1.28E-1 9.82E-2 1.14E-2 7.32E-3 1.13E-2 16 13 6.37E-2 4.86E-2 2.82E-3 2.23E-3 2.83E-3 32 18 3.18E-2 2.43E-2 7.01E-4 6.95E-4 7.05E-4 64 23 1.59E-2 1.21E-2 1.75E-4 2.24E-4 1.76E-4 128 23 7.95E-3 6.06E-3 4.36E-5 7.47E-5 4.40E-5 256 24 3.98E-3 3.03E-3 1.09E-5 2.54E-5 1.09E-5 rate O(h1.01) O(h1.01) O(h2.00) O(h1.65) O(h2.00) Continuous linear mortars on non-matching grids

Department of Mathematics, University of Pittsburgh 25

Error in the computed solution for Example 1

• A. Discontinuous quadratic mortars
• B. Discontinuous linear mortars

Error in computed pressure (shade) and velocity (arrows).

Department of Mathematics, University of Pittsburgh 26

Iterative convergence for Example 2

1/h BalCG CG cond. iter. cond. iter. 4 1.83E+0 5 1.45E+1 8 16 2.49E+1 13 4.29E+1 20 64 2.33E+1 14 1.27E+2 29 256 2.96E+1 15 3.63E+2 45 Continuous quadratic mortars on matching grids 1/h BalCG CG cond. iter. cond. iter. 4 1.79E+1 5 3.91E+1 8 8 1.78E+1 8 3.74E+1 11 16 2.50E+1 13 3.82E+1 19 32 3.68E+1 19 6.60E+1 26 64 4.71E+1 23 1.30E+2 34 128 5.96E+1 24 2.58E+2 51 256 7.29E+1 24 5.16E+2 72 Continuous linear mortars on matching grids

Department of Mathematics, University of Pittsburgh 27

A posteriori error estimates

• Estimate the error by computable quantities
• Use the error estimator to dynamically adapt the grids

E ∈ Th : ω2

E = K−1uh + ∇ph2 Eh2 E + f − ∇ · uh2 Eh2 E + λH − ph2 ∂E∩ΓhE,

τ ∈ T Γ,H : ω2

τ = [uh · n]2 τH3 τ,

˜ ω2

E = h−2 E ω2 E,

˜ ω2

τ = H−2 τ ω2 τ.

Theorem (Upper bounds): p − ph2 ≤ C   

• E∈Th

ω2

E +

• τ∈T Γ,h

ω2

τ

   , u − uh2

H(div) ≤ C

  

• E∈Th

˜ ω2

E +

• τ∈T Γ,h

˜ ω2

τ

   .

Department of Mathematics, University of Pittsburgh 28

Residual-based estimates: lower bounds

Theorem:

• E∈Th

ω2

E +

• τ∈T Γ,H

ω2

τ ≤ C

 p − ph2 +

• E∈Th

h2

Eu − uh2 H(div;E)

+

• τ∈T Γ,H

Hτλ − λH2

τ +

• τ∈T Γ,H

h−1

E,τH3 τu − uh2 H(div;Eτ)

 

• E∈Th

˜ ω2

E +

• τ∈T Γ,H

˜ ω2

τ ≤ C

 

E∈Th

h−2

E p − ph2 E + u − uh2 H(div)

  . Efficient (lower) and reliable (upper) estimate for p − ph: C1  

E∈Th

ω2

E +

• τ∈T Γ,h

ω2

τ

  ≤ p − ph2 ≤ C2  

E∈Th

ω2

E +

• τ∈T Γ,H

ω2

τ

  .

Department of Mathematics, University of Pittsburgh 29

Adaptive mesh refinement algorithm

• 1. Solve the problem on a coarse (both subdomain and mortar) grid.
• 2. For each subdomain Ωi

(a) Compute ωi =  

E∈Th,i

ω2

E +

• τ∈T Γi,h

ω2

τ

 

1/2

(b) If ωi > .5 max1≤j≤n ωj, refine Th,i.

• 3. For each interface Γi,j, if either Ωi or Ωj has been refined m times, refine Th,i,j.
• 4. Solve the problem on the refined grid. If either the desired error tolerance or the

maximum refinement level has been reached, exit; otherwise go to step 2.

Department of Mathematics, University of Pittsburgh 30

Numerical experiments

Example 3: 2D problem with a boundary layer p(x, y) = 1000 x y e−10(x2+y2), K = I Dirichlet BCs; Continuous quadratic mortars Example 4: 2D problem with highly oscilating tensor K =    (105 − 100 sin(20πx) sin(20πy)) ∗ I, 0 ≤ x, y ≤ 1/2 or 1/2 ≤ x, y ≤ 1 (105 − 100 sin(2πx) sin(2πy)) ∗ I,

• therwise

BCs: p|x=0 = 1, p|x=1 = 0, no flow on the rest of the boundary Discontinuous quadratic mortars Multiblock decomposition: 6 × 6 subdomains Coarse grid: 2 × 2 on each subdomain and adaptive mesh refinement

Department of Mathematics, University of Pittsburgh 31

Computed solution for Example 3

Pressure on the fourth grid level

pres 18.00 17.00 16.00 15.00 14.00 13.00 12.00 11.00 10.00 9.00 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00

Discontinuous quadratic mortars

pres 18.00 17.00 16.00 15.00 14.00 13.00 12.00 11.00 10.00 9.00 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00

Discontinuous linear mortars

Department of Mathematics, University of Pittsburgh 32

Computed solution for Example 4

Magnitude of the velocity on the fifth grid level

speed 140.00 130.00 120.00 110.00 100.00 90.00 80.00 70.00 60.00 50.00 40.00 30.00 20.00

Discontinuous quadratic mortars

speed 140.00 130.00 120.00 110.00 100.00 90.00 80.00 70.00 60.00 50.00 40.00 30.00 20.00

Discontinuous linear mortars

Department of Mathematics, University of Pittsburgh 33

Relation to subgrid upscaling methods

Subgrid upscaling (Arbogast et. al.): V0

h,2

V0

h,1

VH Vh = V0

h,1 + V0 h,2 + VH

Multiscale mortar MFE: Vh,2 Vh,1 MH Vh,0 = Vh,1+Vh,2 : [v·n], µ = 0, µ ∈ MH Multiscale mortar MFE: fine scale velocity approximation on each coarse edge; coarse scale non-conforming error. The two methods are related to the two (dual) non-overlapping domain decomposition formulations of Glowinski and Wheeler.

Department of Mathematics, University of Pittsburgh 34

Multiscale mortar method: domain decomposition

BC BC q 1 q2 BC BC BC BC

gH(µ) =

n

• j=1

¯ uh,j · nj, µΓj

λ

aH,j(λ, µ) = −u∗

h,j(λ) · nj, µΓj

aH(λ, µ) =

n

• j=1

aH,j(λ, µ) The solution (uh, ph, λH) to the original problem satisfies the interface problem AHλH = gH

• r

aH(λH, µ) = g(µ), ∀µ ∈ MH, with uh = u∗

h(λH) + ¯

uh, ph = p∗

h(λH) + ¯

ph. CG on the interface: on each iteration compute AHλH = −[u∗

h(λH)·n] by solving

subdomain problems in parallel.

Department of Mathematics, University of Pittsburgh 35

Note on the implementation

It is possible to solve AHλH = gH by solving for the discrete Green’s functions u∗

h(µj H) for each mortar basis function µj H ∈ MH

and forming AH explicitly. In this case cost is comparable to subgrid upscaling and multiscale FEM. The iterative approach is more efficient as long as the number of inerface iterations is less than the number of mortar degrees of freedom per subdomain.

Department of Mathematics, University of Pittsburgh 36

Multiscale flux basis

• φ(k)

H,j

NH,j

k=1 : basis for the mortar space MH,j on Γj.

λH,j =

NH,j

• k=1

λ(k)

H,jφ(k) H,j

Multiscale flux basis: ψ(k)

H,j = AH,jφ(k) H,j,

k = 1, . . . , NH,j Computing ψ(k)

H,j requires solving a subdo-

main problem in Dj with Dirichlet bound- ary data φ(k)

H,j.

! " " "#####"#####"

!"#"$%&'( )* +',-.// 0'123/1'

Using the pre-computed multiscale flux basis in the CG iteration: AH,jλH,j = AH,j  

NH,j

• k=1

λ(k)

H,jφ(k) H,j

  =

NH,j

• k=1

λ(k)

H,jAH,jφ(k) H,j = NH,j

• k=1

λ(k)

H,jψ(k) H,j.

Department of Mathematics, University of Pittsburgh 37

Numerical experiments with multiscale flux basis

• Two types of grid refinement are performed

– Type A: Each subdomain has 10×10 grid; Number of Subdomains is increased

✲

– Type B: Global fine grid is fixed at 120×120 elements; Number of Subdomains is increased

✲

• Three numerical methods are compared

– Method D1 : Original Multiscale Mortar Method Implementation, No pre- conditioner – Method D2 : Original Multiscale Mortar Method Implementation, Balancing Preconditioner – Method D3 : New Multiscale Flux Basis Implementation, No preconditioner

Department of Mathematics, University of Pittsburgh 38

Table 2: Type A - Continuous Linear Mortars (3x1 grids) Method D1 Method D2 Method D3 Subdomains CGIter Solves CGIter Solves CGIter Solves 2x2x1 = 4 11 14 11 41 11 19 3x3x1 = 9 19 22 18 64 19 35 4x4x1 = 16 26 29 22 76 26 35 5x5x1 = 25 31 34 24 82 31 35 6x6x1 = 36 38 41 24 82 38 35 7x7x1 = 49 43 46 25 85 43 35 8x8x1 = 64 50 53 24 82 50 35 Table 3: Type B - Continuous Linear Mortars (3x1 grids) Method D1 Method D2 Method D3 Subdomains CGIter Solves CGIter Solves CGIter Solves 2x2x1 = 4 11 14 12 44 11 19 3x3x1 = 9 21 24 19 67 21 35 4x4x1 = 16 29 32 24 82 29 35 5x5x1 = 25 34 37 26 88 34 35 6x6x1 = 36 40 43 26 88 40 35 7x7x1 = 49 47 50 26 88 47 35 8x8x1 = 64 53 56 26 88 53 35

Department of Mathematics, University of Pittsburgh 39

Extension to two-phase flow

On each subdomain Ωi: Uα = −kα(Sα)K µα ρα(∇Pα − ραg∇D) (Darcy’s law) ∂(φραSα) ∂t + ∇ · Uα = qα (conservation of mass) On each interface Γij: Pα|Ωi = Pα|Ωj, [Uα · n]ij = 0. On each Ωi and Γij: Sw + Sn = 1, pc(Sw) = Pn − Pw.

Department of Mathematics, University of Pittsburgh 40

Domain decomposition

Interface operator BH : MH → MH For λ = (P M

n , P M w ) ∈ MH

BH(λ) = ([UM

n (λ)], [UM w (λ)]),

where UM

α (λ) is the mortar projection of the solution Uα(λ) · n

to subdomain problems with Dirichlet boundary data λ. The original problem is equivalent to solving for λ ∈ MH such that BH(λ) = 0.

Department of Mathematics, University of Pittsburgh 41

Comutational experiment

SPE Comparative Solution Upscaling Project; Oil-water displacement in a hori- zontal cross-section 1200 × 2200 [ft]

5000.00 4736.89 4473.79 4210.68 3947.58 3684.47 3421.37 3158.26 2895.16 2632.05 2368.95 2105.84 1842.74 1579.63 1316.53 1053.42 790.32 527.21 264.11 1.00

Permeability: 2 ∗ 10−3 − 2 ∗ 105 Computational grid 60 × 220 25 blocks: 5 × 5

Department of Mathematics, University of Pittsburgh 42

Computed oil pressure profiles at 2951 days

Three runs: fine grid (1 block), multiblock run with a single linear mortar element

• n each interface (upscaled), multiblock run with refined mortars (6 elements) near

the wells (adapted mortar).

6050 6043 6035 6028 6021 6013 6006 5998 5991 5984 5976 5969 5962 5954 5947 5939 5932 5925 5917 5910

Fine grid solution

6050 6043 6035 6028 6021 6013 6006 5998 5991 5984 5976 5969 5962 5954 5947 5939 5932 5925 5917 5910

Upscaled solution

6050 6043 6035 6028 6021 6013 6006 5998 5991 5984 5976 5969 5962 5954 5947 5939 5932 5925 5917 5910

Adapted mortar solution

Department of Mathematics, University of Pittsburgh 43

Computed oil concentration profiles at 2951 days

COIL 43.00 39.67 36.33 33.00 29.67 26.33 23.00 19.67 16.33 13.00

Fine grid solution

COIL 43.00 39.67 36.33 33.00 29.67 26.33 23.00 19.67 16.33 13.00

Upscaled solution

COIL 43.00 39.67 36.33 33.00 29.67 26.33 23.00 19.67 16.33 13.00

Adapted mortar solution Adapted mortar increases production well rates accuracy by a factor of 2 at the cost of increasing CPU time by 50%.

Department of Mathematics, University of Pittsburgh 44

Comparison of linear and quadratic mortars for two-phase flow

TCOFX 1594.19 1192.77 892.43 667.72 499.59 373.79 279.67 209.25 156.56 117.14 87.64 65.58 49.06 36.71 27.47 20.55 15.38 11.50 8.61 6.44

• A. Permeability field
• B. Grids on the coarsest level

POIL 4631.48 4448.48 4265.48 4082.48 3899.48 3716.48 3533.48 3350.48 3167.48 2984.48 2801.48 2618.48 2435.48 2252.49 2069.49 1886.49 1703.49 1520.49 1337.49 1154.49

• A. Oil pressure

SWAT 0.770 0.743 0.716 0.689 0.662 0.635 0.608 0.581 0.554 0.528 0.501 0.474 0.447 0.420 0.393 0.366 0.339 0.312 0.285 0.258

• B. Water saturation

Time [days] Oil Production Rate [stb/day]

100 200 300 400 500 600 700 50 100 150

• Cont. Linears Mortars
• Cont. Quadratic Mortars

Time [days] Water/Oil Ratio

100 200 300 400 500 600 700 0.02 0.04 0.06 0.08 0.1

• Cont. Linears Mortars
• Cont. Quadratic Mortars

linear mortars quadratic mortars mortar GMRES mortar GMRES 4 51.2 2 42.3 8 72.7 3 54.4 16 92.6 4 63.9 32 155.2 6 75.9 Comparison of recovery curves.

Department of Mathematics, University of Pittsburgh 45

Summary

• Mortar methods are related to multiscale methods
• The multiscale mortar MFE method can be viewed as a generalization of

subgrid upscaling methods.

• Different scales and polynomial degree approximations on interfaces and

subdomains.

• Effective solution: fine subdomain resolution with coarse-grid flux matching
• Optimal fine scale convergence
• Multiscale flux basis implementation can increase efficiency

Department of Mathematics, University of Pittsburgh 46