Probabilistic UQ using PDEs with Random Data: A Case Study Oliver - - PowerPoint PPT Presentation

▶
probabilistic uq using pdes with random data a case study
SMART_READER_LITE
LIVE PREVIEW

Probabilistic UQ using PDEs with Random Data: A Case Study Oliver - - PowerPoint PPT Presentation

Apr 10, 2023 117 likes •665 views

Institut fr Numerische Mathematik und Optimierung Probabilistic UQ using PDEs with Random Data: A Case Study Oliver Ernst December 12, 2011 RICAM Special Semester on Multiscale Simulation & Analysis in Energy and the Environment

slide-1
SLIDE 1

Institut für Numerische Mathematik und Optimierung

Probabilistic UQ using PDEs with Random Data: A Case Study

Oliver Ernst

December 12, 2011 RICAM Special Semester on Multiscale Simulation & Analysis in Energy and the Environment Workshop 4: Multiscale Problems & Stochastic Modelling

slide-2
SLIDE 2

Collaborators, Acknowledgments

Ingolf Busch, Björn Sprungk TU Freiberg, Institute of Numerical Analysis and Optimization Gerald van den Boogaart TU Freiberg, Institute of Stochastics

  • K. Andrew Cliffe

U of Nottingham, School of Mathematical Sciences Elisabeth Ullmann U of Bath, Department of Mathematical Sciences

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 1

slide-3
SLIDE 3

Uncertainty Quantification (UQ)

Goal: Account for the fact that numerical simulations in science and engineering are typically based on data which is uncertain. Enabling factors: mature, sophisticated numerical software, highly parallel computing hardware, growing abundance of data. Here: uncertainty modeled as probabilistic. (There are other ways.) Leads to elliptic BVP with random data. Case study: groundwater flow at WIPP Compare 3 different UQ methods for obtaining quantity of interest (contaminant particle travel time). Emphasis on probabilistic model.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 2

slide-4
SLIDE 4

Uncertainty Quantification (UQ)

Goal: Account for the fact that numerical simulations in science and engineering are typically based on data which is uncertain. Enabling factors: mature, sophisticated numerical software, highly parallel computing hardware, growing abundance of data. Here: uncertainty modeled as probabilistic. (There are other ways.) Leads to elliptic BVP with random data. Case study: groundwater flow at WIPP Compare 3 different UQ methods for obtaining quantity of interest (contaminant particle travel time). Emphasis on probabilistic model.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 2

slide-5
SLIDE 5

Uncertainty Quantification (UQ)

Goal: Account for the fact that numerical simulations in science and engineering are typically based on data which is uncertain. Enabling factors: mature, sophisticated numerical software, highly parallel computing hardware, growing abundance of data. Here: uncertainty modeled as probabilistic. (There are other ways.) Leads to elliptic BVP with random data. Case study: groundwater flow at WIPP Compare 3 different UQ methods for obtaining quantity of interest (contaminant particle travel time). Emphasis on probabilistic model.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 2

slide-6
SLIDE 6

Uncertainty Quantification (UQ)

Goal: Account for the fact that numerical simulations in science and engineering are typically based on data which is uncertain. Enabling factors: mature, sophisticated numerical software, highly parallel computing hardware, growing abundance of data. Here: uncertainty modeled as probabilistic. (There are other ways.) Leads to elliptic BVP with random data. Case study: groundwater flow at WIPP Compare 3 different UQ methods for obtaining quantity of interest (contaminant particle travel time). Emphasis on probabilistic model.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 2

slide-7
SLIDE 7

Outline

WIPP Groundwater Flow Problem Parametrization of Random Fields Solution methods for PDEs with Random Data Numerical Results Concluding Remarks

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 3

slide-8
SLIDE 8

Outline

WIPP Groundwater Flow Problem Parametrization of Random Fields Solution methods for PDEs with Random Data Numerical Results Concluding Remarks

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 4

slide-9
SLIDE 9

WIPP - Waste Isolation Pilot Plant

US DOE repository for radioactive waste situated in New Mexico. Fully operational since 1999. Extensive site characterization and performance assessment since 1976; certification/recertification by US EPA (every 5 years). Large amount of publicly available data.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 5

slide-10
SLIDE 10

WIPP Groundwater Flow Problem

Repository located at a depth of 655m within bedded evaporites, primarily halite (salt). The most transmissive rock in the region is the Culebra Dolomite. In the event of an accidental breach, Culebra would be the principal pathway for transport of radionuclides away from the repository.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 6

slide-11
SLIDE 11

WIPP Scenario

Uncertainty Quantification Problem

Release of radionuclides by means of a borehole drilled into the repository. Transport of radionuclides by groundwater through the Culebra Dolomite. Travel time from release point in the repository to the boundary of the region is an important quantity. Flow is two-dimensional to a good approximation. UQ Problem: permeablilty properties of subsurface available only at few observation points. Approach: Model lack of information stochastically, propagate random input data to travel time.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 7

slide-12
SLIDE 12

Groundwater Flow Model

Stationary Darcy flow q = −K∇p q : Darcy flux K : hydraulic conductivity p : hydraulic head mass conservation ∇· u = 0 u : pore velocity q = φu φ : porosity transmissivity T = Kb b : aquifer thickness particle transport ˙ x(t) = −T(x) bφ ∇p(x) x : particle position x(0) = x0 x0 : release location Quantity of interest: particle travel time to reach WIPP boundary (actually, its log).

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 8

slide-13
SLIDE 13

PDE with Random Coefficient

Primal form of Darcy equations: ∇·[T(x)∇p(x)] = 0, x ∈ D, p = p0 along ∂D. Model T as a random field (RF) T = T(x, ω), ω ∈ Ω, x ∈ D with respect to underlying probability space (Ω, A, P). Assumptions: T has finite mean and covariance T(x) = E [T(x, ·)] , x ∈ D, CovT (x, y) = E

  • T(x, ·) − T(x)

T(y, ·) − T(y)

  • ,

x, y ∈ D. T is lognormal, i.e., Z(x, ω) := log T(x, ω) is a Gaussian RF. CovZ is stationary, isotropic, and of Matérn type.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 9

slide-14
SLIDE 14

WIPP Data

6.05 6.1 6.15 6.2 x 10

5

3.57 3.575 3.58 3.585 3.59 3.595 x 10

6

UTM Easting UTM Northing above median below median WIPP site boundary

transmissivity measurements at 38 test wells head measurements, used to

  • btain boundary data via

statistical interpolation (kriging) constant layer thickness of b = 8m constant porosity of φ = 0.16 SANDIA Nat. Labs reports

[Caufman et al., 1990] [La Venue et al., 1990]

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 10

slide-15
SLIDE 15

Outline

WIPP Groundwater Flow Problem Parametrization of Random Fields Solution methods for PDEs with Random Data Numerical Results Concluding Remarks

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 11

slide-16
SLIDE 16

Parametrization of Random Fields

Merge transmissivity data with statistical model

Given: Measurements of transmissivity Assumption: Matérn covariance structure (contains parameters) Regression model Z(x) = k

i=1 βifi(x) for mean (incorporate

geological information) Procedure: (1) {fi}k

i=1 yield point estimates of parameters in Matérn covariance

function via restricted maximum likelihood estimation (RML). (2) Condition RF Z = log T on transmissivity measurements. (for Gauss RF coincides with kriging, BLUP, leads to low-rank modification of covariance operator.) (3) Mean of conditioned field is kriging interpolant. (4) Approximate log T by truncated Karhunen-Loève expansion.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 12

slide-17
SLIDE 17

Parametrization of Random Fields

Covariance operator

Covariance function of RF Z ∈ L∞(D) ⊗ L2

P(Ω)

c(x, y) = CovZ(x, y) := E

  • Z(x, ·) − Z(x)
  • Z(y, ·) − Z(y)
  • symmetric in x, y ∈ D, positive semidefinite, and continuous on D × D

if continuous along ‘diagonal’ {(x, x) : x ∈ D}. The covariance operator C = CZ : L2(D) → L2(D), (Cu)(x) =

  • D

u(y) c(x, y) dy is therefore selfadjoint, compact, nonnegative. Its eigenvalues {λm}m∈N form a nonincreasing sequence accumulating at most at 0.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 13

slide-18
SLIDE 18

Parametrization of Random Fields

Karhunen-Loève expansion

Exists sequence of RV {ξm}m∈N ⊂ L2

P(Ω),

E [ξm] = 0, E [ξkξm] = δk,m, such that expansion in eigenfunctions {Zm}m∈N Z(x, ω) = Z(x) +

∞

  • m=1
  • λm Zm(x) ξm(ω)

converges in L∞(D) ⊗ L2

P(Ω).

[Karhunen, 1947], [Loève, 1948]

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 14

slide-19
SLIDE 19

Parametrization of Random Fields

Karhunen-Loève expansion and variance

For normalized eigenfunctions Zm(x), VarZ(x) := c(x, x) =

∞

  • m=1

λmZm(x)2, Total variance:

  • D

VarZ(x) dx =

∞

  • m=1

λm (Zm, Zm)D

  • =1

= trace C. For constant variance (e.g., stationary RF), VarZ ≡ σ2 > 0,

  • m

λm = |D| σ2. Interpretation: M first covariance eigenmodes form best rank-M approximation to C in sense of retaining maximal amount of variance.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 15

slide-20
SLIDE 20

Parametrization of Random Fields

Truncated KL expansion

Truncate KL expansion after M leading terms: Z(M)(x, ω) = Z(x) +

M

  • m=1
  • λm Zm(x) ξm(ω).

Truncation error E

  • Z − Z(M)2

L2(D)

  • =

∞

  • m=M+1

λm. Choose M to retain sufficient fraction δ ∈ (0, 1) of total variance, i.e., E

  • Z − Z(M)2

L2(D)

  • E
  • Z2

L2(D)

  • =

∞

m=M+1 λm

∞

m=1 λm

= 1 −

M

m=1 λm

|D|σ2 < δ.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 16

slide-21
SLIDE 21

Matèrn Covariance Functions

c(x, y) = cθ(r) = σ2 2ν−1 Γ(ν)

  • 2√ν r

ρ

ν

Kν

  • 2√ν r

ρ

  • ,

r = x − y2 Kν : modified Bessel function of order ν Parameters θ = (σ2, ρ, ν) σ2 : variance ρ : correlation length ν : smoothness parameter

Special cases: ν = 1

2 :

c(r) = σ2 exp(− √ 2r/ρ) exponential covariance ν = 1 : c(r) = σ2

2r ρ

  • K1
  • 2r

ρ

  • Bessel covariance

ν → ∞ : c(r) = σ2 exp(−r2/ρ2) Gaussian covariance

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 17

slide-22
SLIDE 22

Matèrn Covariance Functions

c(x, y) = cθ(r) = σ2 2ν−1 Γ(ν)

  • 2√ν r

ρ

ν

Kν

  • 2√ν r

ρ

  • ,

r = x − y2 Kν : modified Bessel function of order ν Parameters θ = (σ2, ρ, ν) σ2 : variance ρ : correlation length ν : smoothness parameter

Special cases: ν = 1

2 :

c(r) = σ2 exp(− √ 2r/ρ) exponential covariance ν = 1 : c(r) = σ2

2r ρ

  • K1
  • 2r

ρ

  • Bessel covariance

ν → ∞ : c(r) = σ2 exp(−r2/ρ2) Gaussian covariance

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 17

slide-23
SLIDE 23

Matérn Covariance Functions

1 2 3 4 5 0.2 0.4 0.6 0.8 1 r c(r) ν=1/2 ν=1 ν=2 ν=3 Gauss

ρ = 1

1 2 3 4 5 0.2 0.4 0.6 0.8 1 r c(r) ν=1/2 ν=1 ν=2 ν=3 Gauss

ρ = 3

Smoothness of realizations: RF Z(·, ω) is s times mean-square differentiable if and only if ν > s.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 18

slide-24
SLIDE 24

Matérn Covariance Functions

Eigenvalue asymptotics

Eigenvalue Decay: the smoother the kernel, the faster {λm}m∈N → 0. More precisely: if D ⊂ Rd, then if the kernel function c is piecewise Hs : λm ≤ c1m−s/d piecewise smooth : λm ≤ c2m−s for any s > 0 piecewise analytic : λm ≤ c3 exp(−c4m1/d) for suitable constants c1, c2, c3, c4. Note: piecewise smoothness of kernel also leads to bounds on derivatives of eigenfunctions am in L∞(D). Proven e.g. in [Schwab & Todor, 2006], [Todor, 2006]

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 19

slide-25
SLIDE 25

Matérn Covariance Functions

Eigenvalue preasymptotics

Before asymptotic decay (determined by the smoothness ν of the kernel) sets in, there is a preasymptotic plateau determined by the correlation length parameter ρ.

10 10

1

10

2

10

−6

10

−4

10

−2

10

m λm

ν=1/2, ρ=1 ν=1/2, ρ=1/10 ν=1/2, ρ=1/50 ν=1, ρ=1 ν=1, ρ=1/10 ν=1, ρ=1/50

Eigenvalue decay, Matérn covariance kernel, D = [−1, 1].

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 20

slide-26
SLIDE 26

The Covariance Eigenvalue Problem

Numerical solution

Piecewise constant Galerkin discretization on triangles/tetrahedra. Results in generalized matrix eigenvalue problem Cx = λMx, where M can be chosen diagonal but C is generally dense. Construction/storage of and multiplication with Galerkin matrix via hierarchical matrix technique. Near-field Galerkin blocks using adapted quadrature. Thick-restart Lanczos method for dominant eigenpairs of given

  • rder, extended to block variant for multiple eigenpairs.

Conditioning of covariance on measured data via low-rank correction.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 21

slide-27
SLIDE 27

The Covariance Eigenvalue Problem

Performance

Bessel covariance (ν = 1) on D = [−1, 1]2, ρ = 1/4. Compute leading 300 eigenpairs on sequence of uniformly refined triangular grids, resulting in N = 104 to N = 3 014 656 DOF. Used Krylov space of dimension 500.

10

1

10

2

10

3

10

4

10

5

10

6

10

7

10

−2

10 10

2

10

4

10

6

10

8

N time matrix generation [s]

Bessel kernel Gaussian kernel N log N 10

1

10

2

10

3

10

4

10

5

10

6

10

7

10

−2

10 10

2

10

4

10

6

10

8

N time eigenvalue solve [s]

Bessel kernel Gaussian kernel N

hierarchical matrix generation eigenvalue calculation

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 22

slide-28
SLIDE 28

WIPP KL Modes

Conditioned on 38 transmissivity observations unconditioned, m = 1, 2, 9, 16 conditioned, m = 1, 2, 9, 16

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 23

slide-29
SLIDE 29

WIPP KL Truncation

Portion of captured variance

20 40 60 80 100 120 140 160 180 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

  • No. retained KL modes (M)

Portion of total variance ()

unconditioned conditioned

RML: ˆ σ2 = 3.57, ˆ ρ = 9865, ˆ ν = 0.58.

M δ 5 77.16 % 10 81.28 % 15 83.81 % 50 90.14 % 1000 99.11 %

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 24

slide-30
SLIDE 30

Restricted Maximum Likelihood Estimation

Goal: estimate covariance parameters θ = (σ2, ρ, ν) from observations ζ = {ζj}n

j=1 of Z = log T without introducing bias from estimation of

regression parameters β in model for mean Z(x) =

k

  • i=1

βifi(x) Random vector Z(ω) ∈ Rn at measurement points {xj}n

j=1.

Here Z = Fβ, [F]i,j = fj(xi). Once covariance parameters θ known obtain β by minimizing weighted least squares functional ζ − Fβ2

C −1

θ

:= (ζ − Fβ)⊤C −1

θ (ζ − Fβ)

[Cθ]i,j = cθ(xi, xj), i, j = 1, . . . , n.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 25

slide-31
SLIDE 31

Restricted Maximum Likelihood Estimation

Bias from estimation of β avoided by maximum likelihood estimation of transformed random vector Z ′ := PZ, P : Rn → range(F)⊥ orthogonal projection. Then, regardless of β, E

Z ′ = E

  • P(Fβ +

Z)

  • = 0.

Now maximize likelihood function: joint probability density of Z ′ as function of θ to obtain RML estimate for θ.

[Patterson & Thompson, 1971], [Kitanidis, 1983]

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 26

slide-32
SLIDE 32

Restricted Maximum Likelihood Estimation

Regression functions for WIPP

For x = (x1, x2) ∈ D used Z(x) = β1 + β2x1 + β3x2 + β4d(x) + β5i(x) d(x) : thickness of Culebra overburden at x, i(x) : indicator function for geological zone 1

/&

  • Distinct geological zones at WIPP (Source: WIPP CRA 2009)

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 27

slide-33
SLIDE 33

Restricted Maximum Likelihood Estimation

RML results for different regression functions

constant: σ2 = 3.56, ρ = 9 900, ν = 0.59 linear: σ2 = 1.09, ρ = 2 100, ν = 1.09 linear, overburden: σ2 = 1.15, ρ = 2 200, ν = 1.07 linear, overburden, zone: σ2 = 0.72, ρ = 1 500, ν = 1.22

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 28

slide-34
SLIDE 34

Kriging

Best unbiased linear prediction

Georges Matheron, 1970s: “le krigeage”, “statistical interpolation” Given mean zero random field Z = Z(x, ω) with known covariance structure, observed at n spatial locations Z(xj) = ζj, j = 1, . . . , n. Seek “linear predictor” ˆ Z(x, ω) := m(x)⊤Z(ω), m(x) =

  

m1(x) . . . mn(x)

   ,

Z(ω) =

  

Z(x1, ω) . . . Z(xn, ω)

  

which is unbiased, i.e. E

ˆ

Z(x, ·)

  • = E [Z(x, ·)] for all x ∈ D, and

minimizes the kriging error E

Z(x, ·) − m(x)⊤Z(·) 2

for each x ∈ D.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 29

slide-35
SLIDE 35

Kriging

Best unbiased linear prediction

This results in ˆ Z(x, ω) = c(x)⊤C −1ζ, where c(x) =

  

c(x, x1) . . . c(x, xn)

   ,

[C]i,j = c(xj, xn), ζ =

  

ζ1 . . . ζn

   .

Conditioned covariance function is then the covariance function of the residual field R = Z − ˆ Z: cR(x, y) = c(x)⊤C −1c(y).

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 30

slide-36
SLIDE 36

Kriging

Variants: simple kriging: Z ≡ 0.

  • rdinary kriging: constant

unknown mean universal kriging: non-constant mean obtained by linear regression (used here).

Right: Kriging contours of WIPP trans- missivity data.

6.05 6.1 6.15 6.2 x 10

5

3.57 3.575 3.58 3.585 3.59 3.595 x 10

6

UTM Easting UTM Northing 6.5 6 5.5 5 4.5 4 3.5 Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 31

slide-37
SLIDE 37

Outline

WIPP Groundwater Flow Problem Parametrization of Random Fields Solution methods for PDEs with Random Data Numerical Results Concluding Remarks

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 32

slide-38
SLIDE 38

Solution methods for PDEs with Random Data

Deterministic Parametric Representation

Parametrize input RF by vector of independent RV {ξm}M

m=1 =: ξ.

If ξm has density ρm and Ξm := ξm(Ω), then (Doob-Dynkin lemma) L2

P(Ω) ≃ L2 ρ(ΞM),

where ΞM := ×M

m=1Ξm, ρ = M

  • m=1

ρm. Z(x, ω) ← Z(x, ξ), p(x, ω) ← p(x, ξ), . . . BVP becomes purely deterministic with (possibly) high-dimensional parameter space: ∇·[T(x, ξ)∇p(x, ξ)] = 0, x ∈ D, ρ - a.e., p(x, ξ) = p0(x), x ∈ ∂D, ρ - a.e., where log T(x, ξ) = Z(x) +

M

  • m=1
  • λm Zm(x) ξm.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 33

slide-39
SLIDE 39

Travel-time computations

Generate sufficiently large ensemble of log-travel times s(ξ) = log t(ξ) Compute empirical CDF to quantify uncertainty in travel time. Three sampling methods: (1) Monte Carlo (MC) sampling of RV ξ → s(ξ). (NMC solutions of PDE) (2) Stochastic collocation (SC) → RF representation of velocity field uNSC(x, ξ), use this to sample s(ξ). (NSC solutions of PDE) (3) Gaussian process emulator: NDP MC samples of s(ξ) used to calibrate surrogate of mapping ξ → s(ξ), use this surrogate to sample s(ξ). (NDP solutions of PDE)

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 34

slide-40
SLIDE 40

Stochastic Collocation

For an elliptic BVP with M-dimensional random parameter ξ ∈ ΞM in variational form: find u(ξ) ∈ X such that a(u(ξ), v) = ℓ(v) ∀v ∈ X, (1) stochastic collocation constructs a polynomial approximation uN(ξ) from the solutions u(ξj) at N collocation points {ξj}N

j=1 ⊂ ΞN.

Variants: (1) Full tensor product on ΞM with multivariate Lagrange interpolation (curse of dimensionality). (2) Smolyak sparse grids based on zeros of Hermite polynomials with associated interpolant (used here). (3) Elliptic case: [Tatang & al., 1997], [Mathelin & Husseini, 2003], [Xiu &

Hesthaven, 2005], [Babuška & al., 2007] [Schwab & al., 2009–]

(4) Approximation and convergence result extended to lognormal coefficient: [E. & Sprung, 2011], [Charrier, 2010]

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 35

slide-41
SLIDE 41

Stochastic Collocation

Full tensor vs. Smolyak grids

−10 −8 −6 −4 −2 2 4 6 8 10 −10 −8 −6 −4 −2 2 4 6 8 10

#Nodes = 1089

−10 −8 −6 −4 −2 2 4 6 8 10 −10 −8 −6 −4 −2 2 4 6 8 10

#Nodes = 301

Full-Tensor and Smolyak sparse grids of level q = 5 for Gauss-Hermite abscissas (ni = 2i−1 + 1)

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 36

slide-42
SLIDE 42

Gaussian Process Emulators

An emulator is a statistical approximation to the output of a computer code y = f(x). Basic idea: (1) Represent the code, f(·) as a Gaussian stochastic process f(x, ω) =

n

  • i=1

βihi(x) + Z(x, ω). (2) Run model for sample of design inputs x and observe outputs y. (3) Condition GP on observed outputs y. (4) Emulator provides a distribution function for the output of the computer code. (5) Use emulator as a surrogate for computer model when performing MC analysis.

[Kennedy & O’Hagan, 2001], [Stone, 2011]

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 37

slide-43
SLIDE 43

WIPP FE Discretization

mixed FE discretization, lowest order RT elements for q, pcw. constants for p. fixed mesh, 7124 triangles, respects WIPP boundary flow divergence-free ⇒ discrete fluxes pcw. constant, makes particle trajectory calculation trivial.

0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 x 10

4

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 x 10

4

x y

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 38

slide-44
SLIDE 44

Outline

WIPP Groundwater Flow Problem Parametrization of Random Fields Solution methods for PDEs with Random Data Numerical Results Concluding Remarks

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 39

slide-45
SLIDE 45

Numerical Results

Computational Schemes

(1) Random transmissivity field T(x, ω) from data (2) KL expansion (3) Solve PDE with random diffusion coefficient to obtain realizations

  • f flux field q = q(x, ω)

(4) MC + empirical CDF for log travel time. Step (3) performed by straight MC stochastic collocation Gaussian process emulator

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 40

slide-46
SLIDE 46

Numerical Results

Collocation error

Error measures using MC reference calculation with NMC = 20, 000. expected value of collocation error for head E

  • ph − ph,refL2(D)

1/2

This is the norm in L2(D) ⊗ L2

P(Ω).

norm of expected value of collocation error for head E [p − ph,ref] L2(D) error in variance Var(ph) − Var(ph,ref)L2(D). Same for flux field and log travel time.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 41

slide-47
SLIDE 47

Stochastic Collocation Errors vs. MC solution

10 10

2

10

4

10

4

10

3

10

2

  • rel. error

head 10 10

2

10

4

10

5

10

4

10

3

  • rel. error

head mean 10 10

2

10

4

10

2

10

1

10

  • rel. error

head variance 10 10

2

10

4

10

0.5

10

0.3

10

0.1

  • rel. error

flux 10 10

2

10

4

10

2

10

1

10

  • rel. error

flux mean 10 10

2

10

4

10

  • rel. error

flux variance 10 10

2

10

4

10

3

10

2

10

1

  • rel. error

log tt # nodes 10 10

2

10

4

10

5

10

  • rel. error

log tt mean # nodes 10 10

2

10

4

10

3

10

2

10

1

10 #nodes

  • rel. error

log tt variance SC SC2 SC SC2 SC SC2 SC SC2 SC SC2 SC SC2 SC SC2 SC SC2 SC SC2

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 42

slide-48
SLIDE 48

Travel Time CDFs

M = 5

4 4.5 5 5.5 6 6.5 7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 log10 t ECDF of log travel time, 5 KL modes, 20 000 samples MC SC2 L=2 SC L=2 GE 100

4.35 4.4 4.45 4.5 4.55 4.6 4.65 0.02 0.04 0.06 0.08 0.1 log10 t zoom on early times Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 43

slide-49
SLIDE 49

Travel Time CDFs

M = 10

4 4.5 5 5.5 6 6.5 7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 log10 t ECDF of log travel time 10 KL modes, 20 000 samples MC SC2 L=2 SC L=2 GE 200

4.35 4.4 4.45 4.5 4.55 4.6 4.65 4.7 0.02 0.04 0.06 0.08 0.1 0.12 log10 t

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 44

slide-50
SLIDE 50

Kolmogorov-Smirnov Test

M = 5

Statistical test to determine whether the empirical CDFs were drawn from the same distribution as pure MC with 2K resp. 20K samples. Significance level: α = 0.05.

Surrogate PDE solves KS-test (2K) KS-test (20K) GE 50

  • ✗

100

  • ✗

SC L=0 1 ✗ ✗ SC L=1 11

  • ✗

SC L=2 71

  • SC L=3

351

  • SC2 L=0

1 ✗ ✗ SC2 L=1 11 ✗ ✗ SC2 L=2 61

  • ✗

SC2 L=3 241

  • Oliver Ernst (TU Freiberg)

UQ via PDEs with Random Data: A Case Study Linz, December 2011 45

slide-51
SLIDE 51

Kolmogorov-Smirnov Test

M = 15 Surrogate PDE solves KS-test (2K) KS-test (20K) GE 150 ✗ ✗ 300 ✗ ✗ SC L=0 1 ✗ ✗ SC L=1 31 ✗ ✗ SC L=2 511 ✗

  • SC L=3

5951

  • SC2 L=0

1 ✗ ✗ SC2 L=1 31 ✗ ✗ SC2 L=2 481

  • ✗

SC2 L=3 5021

  • Oliver Ernst (TU Freiberg)

UQ via PDEs with Random Data: A Case Study Linz, December 2011 46

slide-52
SLIDE 52

Neglected Variance

4 4.5 5 5.5 6 6.5 7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 M = 5 M = 10 M = 15 M = 50 M = 1000

Empirical CDF of log t based on KL approximations of different length.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 47

slide-53
SLIDE 53

Outline

WIPP Groundwater Flow Problem Parametrization of Random Fields Solution methods for PDEs with Random Data Numerical Results Concluding Remarks

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 48

slide-54
SLIDE 54

Concluding Remarks

Summary: 3 methods for PDEs with random data for WIPP UQ scenario. Emphasized geostatistical model from data. Collocation convergence analysis for lognormal case. Collocation and emulators show potential for outperforming MC. Covering most of the variance requires more powerful methods. In progress: Anisotropic and sparse tensor collocation methods MC error quantification (Berry-Esseen bounds) Outlook: Bayesian inversion of transmissivity field using SC approximation as surrogate.

Oliver Ernst (TU Freiberg) UQ via PDEs with Random Data: A Case Study Linz, December 2011 49