Multiscale mortar mixed finite element methods for flow in porous - - PowerPoint PPT Presentation

Multiscale mortar mixed finite element methods for flow in porous media

Ivan Yotov

Department of Mathematics, University of Pittsburgh

High Demensional Approximation Seminar University of South Carolina, March 17, 2008

Joint work with Todd Arbogast, Gergina Pencheva, and Mary F. Wheeler, The University of Texas at Austin

Department of Mathematics, University of Pittsburgh 1

Outline

• 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

• Extension to two-phase flow
• Extension to stochastic PDEs

Department of Mathematics, University of Pittsburgh 2

Reservoir rock

Department of Mathematics, University of Pittsburgh 3

Motivation: 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 4

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 5

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 · ν = 0

• n ∂Ωi ∩ ∂Ω

On each interface Γij: pi = pj

• n Γij

[u · ν]ij = 0

• n Γij

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

Department of Mathematics, University of Pittsburgh 6

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 · ν]ijµ = 0, µ ∈ Mh,ij. Subdomain grids do not need to match.

Department of Mathematics, University of Pittsburgh 7

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 · νiΓi = 0, v ∈ Vh,i, (∇ · uh, w)Ωi = (q, w)Ωi, w ∈ Wh,i,

n

• i=1

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

Department of Mathematics, University of Pittsburgh 8

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 9

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 · νiΓi = 0, v ∈ Vh,i, (∇ · uh, w)Ωi = (q, w)Ωi, w ∈ Wh,i,

n

• i=1

uh · νi, µΓ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 10

An approximation result

Weakly continuous velocities: Vh,0 =

• v ∈ Vh :

n

• i=1

v|Ωi · νi, µΓ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 11

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 12

Parallel domain decomposition

Two types of subdomain problems:

BC BC q 1 q2 BC BC BC BC

gH(µ) =

n

• i=1

¯ uh,i · νi, µΓi

λ

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

h,i(λ) · νi, µΓ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 13

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) · ν] 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 · νi → QT

hiuh,i · νi,

AHλH = QT

h1uh,1 · ν1 + QT h2uh,2 · ν2

Department of Mathematics, University of Pittsburgh 14

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 15

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 16

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 17

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 18

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 19

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 · ν]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 20

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 21

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 22

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 23

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 24

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 25

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·ν], µ = 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 26

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 27

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α · ν]ij = 0. On each Ωi and Γij: Sw + Sn = 1, pc(Sw) = Pn − Pw.

Department of Mathematics, University of Pittsburgh 28

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α(λ) · ν

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 29

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 30

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 31

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 32

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 33

Multiscale methods for stochastic PDEs

Joint with B. Ganis, H. Klie, M. F. Wheeler, T. Wildey, and D. Zhang

• Uncertainty in rock properties and sources: model as stochastic functions
• Goal: compute mean (expected) value of the solution
• Solution variance provides measure of uncertainty

D ⊂ Rd, d = 2, 3; Ω - stochastic event space with probability measure P u = −K(x, ω)∇p in D ⊗ Ω ∇ · u = q in D ⊗ Ω p = 0

• n ∂D ⊗ Ω

K(x, ω) - stochastic (random) tensor function ⇒ p(x, ω), u(x, ω) Expected value of a random variable ξ(ω): E[ξ] =

• Ω

ξ(ω)dP(ω) =

• R

yρ(y) dy,

Department of Mathematics, University of Pittsburgh 34

Stochastic modeling of elliptic problems

• Deb, Babuska, Oden (2001): stochastic finite element method
• Babuska, Tempone, Zouraris (2004): hp stochastic FEM
• Xiu, Karniadakis (2003): polynomial chaos
• Zhang, Lu (2004): Karhunen-Loeve for moment equations
• Babuska, Nobile, Tempone (2007); Xiu, Hesthaven (2005); Li, Zhang (2007):

stochastic collocation

• Nobile, Tempone, Webster (2007): sparse grid stochastic collocation
• Lu, Zhang (2007): non-stationary random coefficient

Department of Mathematics, University of Pittsburgh 35

Karhunen-Loeve decomposition

Y = ln(K); Covariance: CY (x, ¯ x) = E[Y (x, ω)Y (¯ x, ω)] - s.p.d. CY (x, ¯ x) =

∞

• n=1

λnfn(x)fn(¯ x); (fm, fn)L2(D) = δmn KL expansion: Y (x, ω) = E[Y ](x) +

∞

• n=1

ξn(ω) √ λnfn(x) ξn - independent Gaussian random variables; ρn(y) =

1 √ 2π exp(−y2/2)

Department of Mathematics, University of Pittsburgh 36

Finite noise assumption

Truncate the KL expansion: Y (x, ω) = E[Y ](x) +

N

• n=1

ξn(ω) √ λnfn(x) Γn = ξn(Ω); Γ = N

n=1 Γn ⊂ RN

Y (x, ω) = Y (x, ξ1(ω), . . . , ξN(ω)) = Y (x, y1, . . . , yN) ρ : Γ → R+ is a joint p.d.f.: ρ(y) = N

n=1 ρn(yn)

P([ξ1, ξ2, . . . , ξN] ∈ γ ⊂ Γ) =

• γ

ρ(y) dy Doob-Dynkn Lemma: u(x, ω) = u(x, ξ1(ω), . . . , ξN(ω)) = u(x, y1, . . . , yN) p(x, ω) = p(x, ξ1(ω), . . . , ξN(ω)) = p(x, y1, . . . , yN)

Department of Mathematics, University of Pittsburgh 37

Weak formulation

V = H(div; D) ⊗ L2(Γ), W = L2(D) ⊗ L2(Γ) v2

V =

• Γ

ρ(y)

• D
• v · v + (∇ · v)2

dx dy = E

• v2

H(div;D)

• w2

W =

• Γ

ρ(y)

• D

w2 dx dy = E

• w2

L2(D)

• Find u(x, y) ∈ V, p(x, y) ∈ W such that
• Γ

(K−1u, v)L2(D)ρ(y) dy =

• Γ

(p, ∇ · v)L2(D)ρ(y) dy, ∀v ∈ V

• Γ

(∇ · u, w)L2(D)ρ(y) dy =

• Γ

(q, w)L2(D)ρ(y) dy, ∀w ∈ W.

Department of Mathematics, University of Pittsburgh 38

Dicretization approaches

• Stochastic FEM: discretize D ⊗ Γ ⊂ Rd+N - high dimensional problem:

Intrusive method

• Monte Carlo Similation:

discretize D and sample Γ randomly; solve a deterministic discrete problem in D for each stochastic realization: Non-intrusive method Slow convergence: Error(p − pM

MC) = O(1/

√ M)

• Stochastic collocation: discretize D and sample Γ in specially chosen (collo-

cation) points: Non-intrusive method with faster convergence.

Department of Mathematics, University of Pittsburgh 39

Stochastic collocation

Semidiscrete formulation Mixed finite element spaces on D: Vh(D) × Wh(D) ⊂ H(div; D) × L2(D)

velocity Lagrange pressure multplier

Find uh : Γ → Vh(D) and ph : Γ → Wh(D) such that a.e. y ∈ Γ, (K−1uh, vh)L2(D) = (ph, ∇ · vh)L2(D), ∀vh ∈ Vh(D) (1) (∇ · uh, wh)L2(D) = (q, wh)L2(D), ∀wh ∈ Wh(D). (2) Collocation points: yk ∈ Γ, k = 1, . . . , M Interpolant: Assume there exists a unique Ipu ∈ Pp(Γ) (polynomial of degree p) such that Ipu(yk) = u(yk), k = 1, . . . , M. Stochastic collocation solution: uh,p ∈ Vh,p := Vh(D) ⊗ Pp(Γ), ph,p ∈ Wh,p := Wh(D) ⊗ Pp(Γ) such that uh,p = Ipuh, ph,p = Ipph.

Department of Mathematics, University of Pittsburgh 40

Stochastic collocation, cont.

Lagrange representation: uh,p(x, y) =

M

• k=1

uh(x, yk)lk(y), ph,p(x, y) =

M

• k=1

ph(x, yk)lk(y), where {lk} is the Lagrange basis: lj(yk) = δjk. To compute uh,p and ph,p, need to solve the deterministic system (1)–(2) for M stochastic realizations yk ∈ Γ.

Department of Mathematics, University of Pittsburgh 41

Tensor-product collocation

Pp(Γ) = Pp1(Γ1) ⊗ · · · ⊗ PpN(ΓN), p = (p1, . . . , pN) Total number of collocation points: M = dim Pp(Γ) = (p1 + 1) × · · · × (pN + 1) Collocation points on Γn: yn,1, . . . yn,pn+1 - the pn + 1 zeros of the orthogonal polynomial qpn+1 with respect to the inner product

• Γn

u(y)v(y)ρn(y) dy Let In

pnu ∈ Ppn(Γn) be the one-dimensional interpolant:

In

pnu(yn,k) = u(yn,k),

k = 1, . . . , pn + 1. Interpolant on Γ: Ipu ∈ Pp(Γ), Ip = I1

p1 ⊗ · · · ⊗ IN pN

Department of Mathematics, University of Pittsburgh 42

Error analysis

Error decomposition: u − uh,pV ≤ u − uhV + uh − uh,pV = u − uhV + uh − IpuhV From deterministic MFEM error analysis, u − uhV =

• Γ

ρ(y)u − uh2

H(div;D) dy

1/2 ≤ Chr+1

• Γ

ρ(y)u2

Hr+1(D) dy

1/2 = Chr+1uHr+1(D)⊗L2(Γ) Interpolation bound on Γ (Babuska, Nobile, Tempone 2007): uh − IpuhV ≤ C

N

• n=1

e−cn√pn, where cn depends on the smoothness of u in Γ.

Department of Mathematics, University of Pittsburgh 43

Error estimate for the pressure

Theorem ˆ p − ph,pW ≤ C(hu − uhV + ph − IpphW) The proof is based on a duality argument: −∇ · K∇ϕ = ˆ p − ph,p, in D, ϕ = 0,

• n ∂D.

Department of Mathematics, University of Pittsburgh 44

Computing expected values

E[ph,p](x) =

• Γ

ph,p(x, y)ρ(y) dy =

• Γ

Ipph(x, y)ρ(y) dy =

• Γ

M

• k=1

ph(x, yk)lk(y)ρ(y) dy =

M

• k=1

wkph(x, yk), where wk =

• Γ

lk(y)ρ(y) dy.

Department of Mathematics, University of Pittsburgh 45

Numerical experiments: flow from left to right

• The method has been implemented in the parallel flow simulators Parcel (Darcy

flow) and SDF (coupled Stokes-Darcy flow).

• Covariance eigenfunctions are computed using a code provided by Don Zhang
• 6 terms in the KL expansion, 4 collocation points in each direction: M = 46 =

4096

P 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

One realization

P 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Mean solution

Department of Mathematics, University of Pittsburgh 46

Flow from left to right: variance

varP 0.0115 0.011 0.0105 0.01 0.0095 0.009 0.0085 0.008 0.0075 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 variance

varU 0.3 0.29 0.28 0.27 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0.19 0.18 0.17 0.16 0.15 0.14 0.13 0.12 0.11 0.1 0.09 0.08

Horizontal vel variance

varV 0.06 0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005

Vertical vel variance

Department of Mathematics, University of Pittsburgh 47

Flow from left to right: convergence

• Convergence in physical space: fixed number of collocation points 54; refine

spatial grid grid Flux Error L2 Error

• Vel. L2 Error

10x10 9.22149E-04 1.11495E-04 9.26733E-04 20x20 2.33581E-04 (3.9479) 2.73432E-05 (4.0776) 2.46538E-04 (3.7590) 40x40 5.59873E-05 (4.1720) 6.50603E-06 (4.2028) 5.99880E-05 (4.1098) 80x80 1.17766E-05 (4.7541) 1.30309E-06 (4.9927) 1.21597E-05 (4.9333)

• Convergence in stochastic space: fixed grid, collocation points 24, 34, 44, 54

colloc points Flux Error L2 Error

• Vel. L2 Error

24 1.1966265E-02 9.8357964E-02 5.85160957E-02 34 2.9685762E-04 4.8594160E-03 1.42482552E-03 44 1.8623395E-05 7.1925780E-04 7.98492911E-05

Department of Mathematics, University of Pittsburgh 48

Quarter five spot

permX 2.262 2.167 2.072 1.977 1.883 1.788 1.693 1.598 1.503 1.409 1.314 1.219 1.124 1.030 0.935 0.840 0.745 0.650

A perm realization

P 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

• 0.05
• 0.1
• 0.15
• 0.2
• 0.25
• 0.3
• 0.35
• 0.4
• 0.45

A solution realization

P 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

• 0.05
• 0.1
• 0.15
• 0.2
• 0.25
• 0.3
• 0.35
• 0.4
• 0.45

Mean solution

varP 0.0044 0.0042 0.004 0.0038 0.0036 0.0034 0.0032 0.003 0.0028 0.0026 0.0024 0.0022 0.002 0.0018 0.0016 0.0014 0.0012 0.001 0.0008 0.0006 0.0004 0.0002

Pressure variance

varU 0.005 0.0046 0.0042 0.0038 0.0034 0.003 0.0026 0.0022 0.0018 0.0014 0.001 0.0006 0.0002

Horizontal vel variance

varV 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

Vertical vel variance

Department of Mathematics, University of Pittsburgh 49

Discontinuous permeability

permX 220 190 160 130 100 70 40 10 1.86207 1.65517 1.44828 1.24138 1.03448 0.827586 0.62069 0.413793 0.206897

A perm realization

P 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

A solution realization

P 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

Mean solution

varP 0.0044 0.0042 0.004 0.0038 0.0036 0.0034 0.0032 0.003 0.0028 0.0026 0.0024 0.0022 0.002 0.0018 0.0016 0.0014 0.0012 0.001 0.0008 0.0006 0.0004 0.0002

Pressure variance

varU 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100

Horizontal vel variance

varV 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100

Vertical vel variance

Department of Mathematics, University of Pittsburgh 50

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
• Extension to stochastic PDEs

Current and future work

• Sparse grid stochastic collocation
• Piecewise stochastic collocation
• Different KL expansions in different subdomains
• A posteriori error analysis for stochastic multiscale methods

Department of Mathematics, University of Pittsburgh 51