W HAT IS GOING TO HAPPEN TO THE POLAR ICE CAPS ? They affect - - PowerPoint PPT Presentation

▶
w hat is going to happen to the polar ice caps
SMART_READER_LITE
LIVE PREVIEW

W HAT IS GOING TO HAPPEN TO THE POLAR ICE CAPS ? They affect - - PowerPoint PPT Presentation

Jan 16, 2023 317 likes •1.08k views

S CALABLE S OLVERS FOR F ORWARD AND I NVERSE P ROBLEMS IN G EOPHYSICS Toby Isaac tisaac@uchicago.edu University of Chicago January 11, 2016 CAAM Department Rice University T. Isaac (U. Chicago) Scalable Forward & Inverse Solvers January

slide-1
SLIDE 1

SCALABLE SOLVERS FOR FORWARD AND INVERSE PROBLEMS IN GEOPHYSICS

Toby Isaac tisaac@uchicago.edu

University of Chicago

January 11, 2016 CAAM Department Rice University

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 1 / 58

slide-2
SLIDE 2

1

MOTIVATION

2

LARGE-SCALE BAYESIAN INFERENCE FOR PARAMETER FIELDS

3

SCALABLE SOLVERS FOR NONLINEAR STOKES EQUATIONS

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 2 / 58

slide-3
SLIDE 3

Motivation Bayesian Inference Stokes Solvers References

MY INTERESTS

1 Adaptive mesh refinement data structures and algorithms for

large-scale scientific computing (slides on my website)

2 Scientific software development:

p4est (including hybrid extension p6est), PETSc (including stand-alone plugins like DofColumns), Deal.II, mangll, . . . 2016 Scientific Software Days, February 25-26, Austin! scisoftdays.org

3 Uncertainty quantification, Bayesian inversion, PDE-constrained

  • ptimization

4 Scalable, efficient solvers for implicit PDEs

(3) and (4) today.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 3 / 58

slide-4
SLIDE 4

Motivation Bayesian Inference Stokes Solvers References

WHAT IS GOING TO HAPPEN TO THE POLAR ICE CAPS?

They affect temperature, sea level (70 m equivalent), ocean currents, etc. Quantities of Interest (QoIs): q visibleearth.nasa.gov

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 4 / 58

slide-5
SLIDE 5

Motivation Bayesian Inference Stokes Solvers References

WHAT IS GOING TO HAPPEN TO THE POLAR ICE CAPS?

They affect temperature, sea level (70 m equivalent), ocean currents, etc. Quantities of Interest (QoIs): q earthobservatory.nasa.gov

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 4 / 58

slide-6
SLIDE 6

Motivation Bayesian Inference Stokes Solvers References

WHAT IS GOING TO HAPPEN TO THE POLAR ICE CAPS?

They affect temperature, sea level (70 m equivalent), ocean currents, etc. Quantities of Interest (QoIs): q earthobservatory.nasa.gov

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 4 / 58

slide-7
SLIDE 7

Motivation Bayesian Inference Stokes Solvers References

WHAT IS GOING TO HAPPEN TO THE POLAR ICE CAPS?

They affect temperature, sea level (70 m equivalent), ocean currents, etc. Quantities of Interest (QoIs): q earthobservatory.nasa.gov

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 4 / 58

slide-8
SLIDE 8

Motivation Bayesian Inference Stokes Solvers References

SCIENTISTS GATHER LOTS OF DATA

Observations/Data: dobs

(NASA Earth Observatory)

Ice cores / boreholes, Ice-penetrating radar (+stratigraphy), GRACE gravity field measurements, Interferometric synthetic aperture radar (InSAR) / lidar, etc.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 5 / 58

slide-9
SLIDE 9

Motivation Bayesian Inference Stokes Solvers References

SCIENTISTS GATHER LOTS OF DATA

Observations/Data: dobs Ice cores / boreholes, Ice-penetrating radar (+stratigraphy), GRACE gravity field measurements, Interferometric synthetic aperture radar (InSAR) / lidar, etc.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 5 / 58

slide-10
SLIDE 10

Motivation Bayesian Inference Stokes Solvers References

SCIENTISTS GATHER LOTS OF DATA

Observations/Data: dobs

(UT Center for Space Research)

Ice cores / boreholes, Ice-penetrating radar (+stratigraphy), GRACE gravity field measurements, Interferometric synthetic aperture radar (InSAR) / lidar, etc.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 5 / 58

slide-11
SLIDE 11

Motivation Bayesian Inference Stokes Solvers References

SCIENTISTS GATHER LOTS OF DATA

Observations/Data: dobs

(NASA/JPL-Caltech)

Ice cores / boreholes, Ice-penetrating radar (+stratigraphy), GRACE gravity field measurements, Interferometric synthetic aperture radar (InSAR) / lidar, etc.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 5 / 58

slide-12
SLIDE 12

Motivation Bayesian Inference Stokes Solvers References

HOW CAN WE USE DATA TO GET QOIS?

Pure data analysis Use statistical tools and known values of q to understand/project the relationship between dobs and predicted values of q: E[q|dobs]; Var(q|dobs); etc. Pure deterministic inversion Use unified models of processes explaining dobs and q (united by parameter set m) to invert for m, then project: Fobs : m → d; FQoI : m → q; q := FQoI(“F −1

  • bs(dobs)”).

Model-based Uncertainty Quantification Use models to better inform statistical relationships between dobs and q: E[q = FQoI(m)|dobs]; Var[. . . ]; etc.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 6 / 58

slide-13
SLIDE 13

Motivation Bayesian Inference Stokes Solvers References

PDE-CONSTRAINED FOBS

The observations dobs belong to a vector space RNd. The observations are made on the state w ∈ X via b : X → RNd w is implicitly defined by the forward (PDE) model A parameterized by m: A(w; m) = 0. Fobs : m → b(w) s.t. A(w; m) = 0.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 7 / 58

slide-14
SLIDE 14

Motivation Bayesian Inference Stokes Solvers References

PDE-CONSTRAINED FOBS: ICE SHEET PROBLEM

[Rignot et al., 2011a] [Creyts and Schoof, 2009] Our dobs is the MEaSUREs surface velocity dataset for Antarctica [Rignot et al., 2011b]. Our model parameter m is log β, where β is the sliding coefficient between ice and substrate (units stress velocity−1):

phenomenological; highly spatially variable.

Our state variables w are the velocity u and pressure p in the ice sheet. Our model A is a nonlinear Stokes boundary value problem, where m is a coefficient field for a Robin-type boundary condition.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 8 / 58

slide-15
SLIDE 15

Motivation Bayesian Inference Stokes Solvers References

ICE SHEET DYNAMICS (INSTANTANEOUS)

Balance of linear momentum and mass −∇ · [2µ(T, u) D(u) − Ip] = ρg, [D(u) = 1

2(∇u + ∇uT )]

∇ · u = 0, Constitutive relations µ(T, u) = 1

2A(T)(D(u)II + ǫ)

1−n 2n

[D(u)II = 1

2tr(D(u)2)]

Boundary conditions σn|Γtop = 0, u · n|Γbase = 0, T

(σn + β(m)u)|Γbase = 0

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 9 / 58

slide-16
SLIDE 16

Motivation Bayesian Inference Stokes Solvers References

ICE SHEET DYNAMICS (INSTANTANEOUS)

Balance of linear momentum and mass −∇ · [2µ(T, u) D(u) − Ip] = ρg, [D(u) = 1

2(∇u + ∇uT )]

∇ · u = 0, Constitutive relations shear thinning with second invariant µ(T, u) = 1

2A(T)(D(u)II + ǫ)

1−n 2n

[D(u)II = 1

2tr(D(u)2)]

Boundary conditions σn|Γtop = 0, u · n|Γbase = 0, T

(σn + β(m)u)|Γbase = 0

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 9 / 58

slide-17
SLIDE 17

Motivation Bayesian Inference Stokes Solvers References

ICE SHEET DYNAMICS (INSTANTANEOUS)

Balance of linear momentum and mass −∇ · [2µ(T, u) D(u) − Ip] = ρg, [D(u) = 1

2(∇u + ∇uT )]

∇ · u = 0, Constitutive relations µ(T, u) = 1

2A(T)(D(u)II + ǫ)

1−n 2n

[D(u)II = 1

2tr(D(u)2)]

Boundary conditions stress-free surface σn|Γtop = 0, u · n|Γbase = 0, T

(σn + β(m)u)|Γbase = 0

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 9 / 58

slide-18
SLIDE 18

Motivation Bayesian Inference Stokes Solvers References

ICE SHEET DYNAMICS (INSTANTANEOUS)

Balance of linear momentum and mass −∇ · [2µ(T, u) D(u) − Ip] = ρg, [D(u) = 1

2(∇u + ∇uT )]

∇ · u = 0, Constitutive relations µ(T, u) = 1

2A(T)(D(u)II + ǫ)

1−n 2n

[D(u)II = 1

2tr(D(u)2)]

Boundary conditions basal sliding σn|Γtop = 0, u · n|Γbase = 0, T

(σn + β(m)u)|Γbase = 0

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 9 / 58

slide-19
SLIDE 19

Motivation Bayesian Inference Stokes Solvers References

GENERAL BAYESIAN INVERSION FRAMEWORK

Assume for the moment a finite-dimensional parameter space Mh. The likelihood of dobs given m is πlike(dobs| Fobs(m)). Bayes’ law: the posterior probability of the parameter given the

  • bservations is

πpost(m|dobs) = πprior(m)πlike(dobs| Fobs(m))

  • Mh πprior(m)πlike(dobs| Fobs(m)) dm.
  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 10 / 58

slide-20
SLIDE 20

Motivation Bayesian Inference Stokes Solvers References

GENERAL BAYESIAN INVERSION FRAMEWORK

Thus, e.g., E[FQoI(m)|dobs] =

  • Mh FQoI(m)πprior(m)πlike(dobs| Fobs(m)) dm
  • Mh πprior(m)πlike(dobs| Fobs(m)) dm

. Expectations require integration over Mh (high-dimensional) . . . . . . w.r.t. πpost(m|dobs) (implicitly-defined, cannot be sampled directly). General solution: Metropolis-Hasting. Constructs a Markov chain of samples {mi} such that lim

N→∞

N

i=1 FQoI(mi)

N = E[FQoI(m)|dobs] a.s.

Good performance requires drawing samples from a proposal distribution similar to πpost(m) to generate a sample chain.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 11 / 58

slide-21
SLIDE 21

Motivation Bayesian Inference Stokes Solvers References

GAUSSIAN ASSUMPTIONS

Assume Gaussian prior knowledge of the parameter, πprior(m) = N(mprior, Cprior). Assume the likelihood is a conditional probability given by πobs(dobs| Fobs(m)) = N(Fobs(m), Cobs). Cobs ∈ RNd×Nd incorporates noise in the observations and model error. The posterior probability has the form πpost(m|dobs) ∝ exp{− 1

2(dobs − Fobs(m)2 C−1

  • bs + m − mprior2

C−1

prior)

  • J (m):=

}. J (m) is a typical objective function in PDE-constrained optimization.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 12 / 58

slide-22
SLIDE 22

Motivation Bayesian Inference Stokes Solvers References

LAPLACE’S APPROXIMATION

ˆ πpost(m) := N(mMAP, H(mMAP)−1)

exp(−J (m)) ⇓

⇑ − log(πpost(m))

πpost ˆ πpost

A quadratic fit at the maximum a posteriori likelihood (MAP) point mMAP. In addition to a proposal density, also provides one-shot QoI estimates: E[q] ≈ FQoI(mMAP); Var(q) ≈ DFQoIH(mMAP)−1DF ∗

QoI.

Requires the Hessian H(m) := D2J (m).

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 13 / 58

slide-23
SLIDE 23

Motivation Bayesian Inference Stokes Solvers References

COMPUTING WITH ˆ πPOST(m)

What does it mean to compute with ˆ πpost(m)?

1 Evaluate ˆ

πpost(m) at a given m?

2 Compute mean mMAP? 3 Compute covariance Cpost := H(mMAP)−1? 4 Compute precision C−1

post := H(mMAP)?

5 Compute covariance-/precision-vector products? 6 Draw samples from ˆ

πpost(m)?

Which are necessary?

For Metropolis-Hastings, at least (1, up to a constant) and (6).

For (1), evaluate J (m) → evalute Fobs(m) → solve A(w; m) = 0. For (6) classical approach: symmetric factorization Cpost = LL∗, y ← mMAP + Lz, z ∼ N(0, I). At least (2) is necessary.

Which are intractable for large problems?

(3) and (4). We can’t compute L if we can’t compute Cpost: how to sample?

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 14 / 58

slide-24
SLIDE 24

Motivation Bayesian Inference Stokes Solvers References

COMPUTING THE MAP POINT

Armijo-Newton-Eisenstat-Walker-Steihaug(-Hestenes-Stiefel-Krylov-. . . ): Armijo: Update is mi+1 = mi + α ˆ m, Where ˆ m is descent direction for J (m) and α ∈ (0, 1) ensures sufficient decrease. Newton: Compute ˆ m by H(mi) ˆ m = −DJ (m). Eisenstat-Walker: Adaptive tolerance for inexact linear solver based

  • n nonlinear convergence history.

Steihaug: Conjugate gradient method with early termination for negative curvature (ensure descent direction). Additional requirements: gradient DJ (m), Hessian-vector product, (s.p.d. preconditioner for H(m)−1?)

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 15 / 58

slide-25
SLIDE 25

Motivation Bayesian Inference Stokes Solvers References

MAP POINT: ANTARCTICA

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 16 / 58

slide-26
SLIDE 26

Motivation Bayesian Inference Stokes Solvers References

GRADIENT DJ (m) VIA ADJOINTS

Let L(m, w, z) be the Lagrangian of J (m), L(m, w, z) := 1

2dobs − b(w)2 C−1

  • bs
  • :=Jmisfit(w)

+ 1

2m − mprior2 C−1

prior + z∗A(w; m),

and let w and z solve Lz(m, w, z) = 0 ⇔ A(w; m) = 0, Lw(m, w, z)∗ = 0 ⇔ Aw(w; m)∗z

  • linear PDE

= −DJmisfit(w). Then the gradient is DJ (m) = Lm(m, w, z).

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 17 / 58

slide-27
SLIDE 27

Motivation Bayesian Inference Stokes Solvers References

ADJOINT Aw(w; m)∗z = −DJMISFIT(m) IN STRONG FORM

Balance of linear momentum and mass −∇ ·

  • 2µ′(T, u) D(v) − Iq
  • = 0,

∇ · v = 0, Constitutive relations µ′(T, u) = µ(T, u)(I ⊗ I + 1−n

2n (D(u)II + ǫ)−1D(u) ⊗ D(u)),

Boundary conditions σn|Γtop = −fC−1

  • bs(u − dobs),

u · n|Γbase = 0, T

(σn + β(m)u)|Γbase = 0.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 18 / 58

slide-28
SLIDE 28

Motivation Bayesian Inference Stokes Solvers References

HESSIAN VIA ADJOINTS

Second partial derivatives of L: B := D2Jmisfit = b∗C−1

  • bsb,

C := Am, D := z∗Aww, E := z∗Amw, F := z∗Amm. Then, given ˆ m,

1 Solve Aw(u)ˆ

w = −C ˆ m.

2 Solve Aw(u)∗ˆ

z = −(D + B)ˆ w − E∗ ˆ m.

3 H(m) ˆ

m = (C−1

prior + F) ˆ

m + C∗ˆ z + E ˆ w. Or, in other words, H(m) = C−1

prior + F + C∗A−∗ w (D + B)A−1 w C − C∗A−∗ w E∗ − EA−1 w C

=: C−1

prior + Hmisfit(m).

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 19 / 58

slide-29
SLIDE 29

Motivation Bayesian Inference Stokes Solvers References

TO SUMMARIZE

Everything reduces to

1 Compute J (m)

[Solve A(w; m) = 0.]

2 Compute DJ (m)

[Solve A∗

wz systems.] 3 Compute H(m) ˆ

m [Solve Aw ˆ w, A∗

wˆ

z systems.]

4 Preconditioning H(m)−1. 5 Sampling a Gaussian with covariance Cpost = H(mMAP)−1.

Recall H(m) = C−1

prior + Hmisfit(m).

In [Isaac et al., 2015b], we demonstrate for the Antarctic ice sheet Bayesian inference problem that Cprior is (asymptotically) good for (4), (5) can be efficiently done for Cpost using low-rank updates provided it can be efficiently done for Cprior, and defining Cprior via a differential precision operator is efficient and scalable for (4) and (5).

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 20 / 58

slide-30
SLIDE 30

Motivation Bayesian Inference Stokes Solvers References

GAUSSIAN PRIORS IN INFINITE-DIMENSIONS

To ensure correctness as h → 0, we consider an infinite-dimensional parameter space M. We treat M as an abstract Wiener space: Cprior = EE∗; E embeds a Cameron-Martin space. E converts a cylinder-set measure on C-M space to true measure on M. Cprior must be a Hilbert-Schmidt operator, ergo compact. C−1

prior is only densely-defined: a probability density is not defined, only

a probability measure. The prior-preconditioned Hessian is CpriorH(m) = I + CpriorHmisfit(m), a compact perturbation of identity. Therefore the condition number of CpriorH(m) is uniformly bounded as h → 0.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 21 / 58

slide-31
SLIDE 31

Motivation Bayesian Inference Stokes Solvers References

PRIOR-PRECONDITIONED HESSIAN SPECTRA, ANTARCTICA

1000 2000 3000 4000 10

−2

10 10

2

10

4

10

6

10

8

number eigenvalue 409,545 parameters 1,190,403 parameters

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 22 / 58

slide-32
SLIDE 32

Motivation Bayesian Inference Stokes Solvers References

TWO WAYS TO THINK ABOUT POSTERIOR SAMPLING

  • 1. PRIOR COVARIANCE FACTORIZATION

Let Cprior = LL∗; drop mMAP arguments from H, Hmisfit.

1 Factorize:

C−1

post = C−1 prior + Hmisfit = L−∗(I + L∗HmisfitL)L−1. 2 EVD:

L∗HmisfitL =: V ΛV ∗ [V ∗V = I].

3 Sherman-Morrison-Woodbury:

Cpost = L(I + V ΛV ∗)−1L∗ = L(I − V (V ∗V + Λ−1)−1V ∗)L∗ = LV (I + Λ)−1V ∗L∗.

4 Sample:

y = mMAP + LV (I + Λ)−1/2z, z ∼ N(0, I).

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 23 / 58

slide-33
SLIDE 33

Motivation Bayesian Inference Stokes Solvers References

TWO WAYS TO THINK ABOUT POSTERIOR SAMPLING

  • 2. PRIOR COVARIANCE NORM

1 G-EVD:

Hmisfit =: ˜ V Λ ˜ V ∗ [ ˜ V ∗Cprior ˜ V = I].

2 Sherman-Morrison-Woodbury: let W = Cprior ˜

V (and thus W ˜ V ∗ = I), Cpost = (C−1

prior + Hmisfit)−1

= (I + Cprior ˜ V Λ ˜ V ∗)−1Cprior = (I − W(I + Λ−1)−1 ˜ V ∗)Cprior = W(I − Λ(I + Λ)−1)W ∗ = W(I + Λ)−1/2 ˜ V ∗Cprior ˜ V (I + Λ)−1/2W ∗.

3 Sample:

y = mMAP + W(I + Λ)−1/2 ˜ V ∗z, z ∼ N(0, Cprior).

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 24 / 58

slide-34
SLIDE 34

Motivation Bayesian Inference Stokes Solvers References

TWO WAYS TO THINK ABOUT POSTERIOR SAMPLING

Why prefer one over the other? Philosophically different:

(1) observations change the Cameron-Martin space; (2) observations change the embedding of the Cameron-Martin space.

Practically different:

(1) assumes an invertible factorization of Cprior; (2) treats sampling from Cprior as a black box: allows, e.g.,

rectangular symmetric factorizations of Cprior; sampling Cprior not based on a matrix-vector product [Parker and Fox, 2012].

You may object: WW ∗ = Cprior, we end up computing an invertible factorization of Cprior anyway!

That’s if we do the complete eigen-decomposition . . .

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 25 / 58

slide-35
SLIDE 35

Motivation Bayesian Inference Stokes Solvers References

LOW-RANK APPROXIMATION

By the known compactness of CpriorHmisfit, limh→0 tr(Λ) < ∞. For a given tolerance ǫ, there is an r such that the leading r eigenvalues Λr are sufficient, lim suph→0 tr(Λ) − tr(Λr) < ǫ.

r is how much we are able to learn from the observations.

Using randomized methods [Halko et al., 2011], we can approximate ˜ V Λ ˜ V ∗ ≈ ˜ VrΛr ˜ V ∗

r

in O(r) + O(1) Hmisfit-vector products (depends on eigenvalue separation). We modify our sampling operation: y = mMAP + (W(I + Λ)−1/2 ˜ V ∗ − I + I)z, = mMAP + (W(I + Λ)−1/2 ˜ V ∗ − W ˜ V ∗ + I)z, ≈ mMAP + (Wr((Ir + Λr)−1/2 − Ir) ˜ V ∗

r + I)z.

Only computes as much factorization (Wr) as necessary.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 26 / 58

slide-36
SLIDE 36

Motivation Bayesian Inference Stokes Solvers References

CPRIOR AND CPOST SAMPLES, ANTARCTICA

Prior: Posterior:

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 27 / 58

slide-37
SLIDE 37

Motivation Bayesian Inference Stokes Solvers References

DIFFERENTIAL PRECISION OPERATORS

Cprior needs three scalable, efficient operations:

apply C−1

prior,

apply Cprior, and sample N(0, Cprior).

Encoding C−1

prior as an elliptic differential operator (e.g., FEM

discretization) is appealing:

C−1

prior is a sparse matrix;

Known optimal solvers (e.g. multigrid) to apply Cprior; Familiar similarity to Tikhonov regularization for PDE-constrained

  • ptimization;

C−1

prior has a sparse, rectangular factorization RR∗ via quadrature points,

so samples can be quickly generated with CpriorR; One catch: Cprior is only Hilbert-Schmidt in L2(Ω) for sufficiently smooth operators: e.g., bi-Laplacian in 2D/3D.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 28 / 58

slide-38
SLIDE 38

Motivation Bayesian Inference Stokes Solvers References

BAYESIAN INFERENCE SUMMARY

We are able to scalably approximate the posterior distribution of a Bayesian inverse problem for spatially variable parameter fields with O(106) parameters and forward models defined by large-scale nonlinear PDEs (O(108) dofs). The key components are: Scalable efficient solvers for forward, linearized, and adjoint PDE models (next section) for J (m), DJ (m) and H(m) calculations. A scalable PDE-constrained optimization solver to compute mMAP. A randomized, approximate, generalized eigenvalue computation to find the parameter modes most informed (w.r.t. the prior covariance) by the observations. A prior precision (C−1

prior) representation with fast operations for

precision-vector and covariance-vector products and sampling.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 29 / 58

slide-39
SLIDE 39

Motivation Bayesian Inference Stokes Solvers References

DIRECTIONS

H-matrix covariance approximation for highly-data-informed ˆ πpost(m). Non-Gaussian Proposal Distributions:

−0.5 0.5 1 −0.5 0.5 1 x y

Create maps that deform Gaussian-generated samples to match [Moselhy and Marzouk, 2012]. Compose with low-rank approximations to select directions of deformation.

Multilevel Uncertainty Quantification:

Multilevel Monte Carlo is well-developed theory for forward uncertainty propagation Multilevel Markov Chain Monte Carlo for inference is relatively new [Ketelsen et al., 2013]

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 30 / 58

slide-40
SLIDE 40

Motivation Bayesian Inference Stokes Solvers References

ICE SHEETS AND MANTLE: ymir

We use one FEM framework for both ice sheets and mantle Similarities: both sheer-thinning viscous flows with temperature-dependent viscosity Differences: ice sheet flow is gravity/geometry driven; mantle flow is convection driven; boundary conditions; degree of viscosity variation. Ymir: Norse ice giant formed from rime melted by fire

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 31 / 58

slide-41
SLIDE 41

Motivation Bayesian Inference Stokes Solvers References

ICE SHEETS AND MANTLE: ymir

We use one FEM framework for both ice sheets and mantle Similarities: both sheer-thinning viscous flows with temperature-dependent viscosity Differences: ice sheet flow is gravity/geometry driven; mantle flow is convection driven; boundary conditions; degree of viscosity variation. Ymir: Norse ice giant formed from rime melted by fire

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 31 / 58

slide-42
SLIDE 42

Motivation Bayesian Inference Stokes Solvers References

NAÏVITÉ

ymir follows rhea, a mantle convection code [Burstedde et al., 2013]: Q1 × Q1 stabilized (Bochev-Dohrmann) FEM discretization Picard iteration for nonlinear solver MINRES iteration for Krylov solver hypre/ML algebraic multigrid for elliptic (1,1)-block “Great! Let’s do this for ice, with high-order FEM and Newton. I’m sure the linear solver will still work.” –Me

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 32 / 58

slide-43
SLIDE 43

Motivation Bayesian Inference Stokes Solvers References

NAÏVITÉ

ymir follows rhea, a mantle convection code [Burstedde et al., 2013]: Q1 × Q1 stabilized (Bochev-Dohrmann) FEM discretization Picard iteration for nonlinear solver MINRES iteration for Krylov solver hypre/ML algebraic multigrid for elliptic (1,1)-block “Great! Let’s do this for ice, with high-order FEM and Newton. I’m sure the linear solver will still work.” –Me (circa 2010)

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 32 / 58

slide-44
SLIDE 44

Motivation Bayesian Inference Stokes Solvers References

HIGH-ORDER MIXED SPACES

Inf-sup stable discretizations with discontinuous pressure spaces that satisfy incompressibility at the element level. [Qk]3 Primary space: Fast matrix-free matvec Pdisc

k−1 Incompressibility constraint

space: inf-sup constant bounded

  • ver h and K

balanced approximation

  • rders

inf-sup degrades with aspect ratio Both families inf-sup stable on octree-refined hexahedra [Heuveline and Schieweck, 2007].

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 33 / 58

slide-45
SLIDE 45

Motivation Bayesian Inference Stokes Solvers References

HIGH-ORDER MIXED SPACES

Inf-sup stable discretizations with discontinuous pressure spaces that satisfy incompressibility at the element level. [Qk]3 Primary space: Fast matrix-free matvec Qdisc

k−2 Incompressibility constraint

space: inf-sup constant bounded

  • ver h, O(k−1)

inf-sup constant unaffected by aspect ratio [Toselli and Schwab, 2003] lose one order of approximation Both families inf-sup stable on octree-refined hexahedra [Heuveline and Schieweck, 2007].

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 33 / 58

slide-46
SLIDE 46

Motivation Bayesian Inference Stokes Solvers References

SOLVING THE NONLINEAR STOKES EQUATIONS

Armijo-Newton-Eisenstat-Walker-Saad-Schur-??? Saad: FGMRES(k) Krylov solver allows use of variable Schur: Use block upper-triangular preconditioner for Stokes operator: A = F B∗ B

  • ,

P = ˜ F B∗ ˜ S

  • ,

where ˜ F preconditions elliptic (1,1)-block, ˜ S preconditions Schur complement S = −BF −1B∗.

What should the block preconditioners ˜ F and ˜ S be?

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 34 / 58

slide-47
SLIDE 47

Motivation Bayesian Inference Stokes Solvers References

ICE SHEET PRECONDITIONING CHOICES

˜ F:

Antarctica has a complicated geometry: coarsest unstructured mesh has at least O(105) high-order elements. Without geometric coarsening, geometric multigrid would require very good coarse solve. (1,1)-block is elliptic with rigid-body modes in near null space (in the absence of boundary conditions): Smoothed aggregation algebraic multigrid is good for this type of problem (implemented by ML and PETSc’s PCGAMG). Design space:

Which smoother to use on each level of the hierarchy? How aggressively to aggregate degrees of freedom?

˜ S:

For constant / smoothly varying viscosity, S ∼ −µ−1M. ˜ S = −µ−1 ˆ M: lumped mass matrix.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 35 / 58

slide-48
SLIDE 48

Motivation Bayesian Inference Stokes Solvers References

ANISOTROPY IS IMPORTANT TO SMOOTHER EFFECTIVENESS

20 40 60 80 100 120 140 160 180 10−16 10−14 10−12 10−10 10−8 10−6 10−4 10−2 100

GMRES(30) iterations ||rk|//||r0|| Smoother: symmetric Gauss-Seidel 1:1 10:1 100:1

Pointwise smoothers do not smooth residuals in all directions of anisotropic problems.

Typical width:height aspect ratios for ice sheet discretizations are on the order of 100:1.

Example: µ = β = 1, Q3 preconditioned by Q1, Gauss-Seidel

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 36 / 58

slide-49
SLIDE 49

Motivation Bayesian Inference Stokes Solvers References

ANISOTROPY IS IMPORTANT TO SMOOTHER EFFECTIVENESS

10 20 30 40 50 60 70 10−16 10−14 10−12 10−10 10−8 10−6 10−4 10−2 100

GMRES(30) iterations ||rk||/||r0|| Smoother: ILU(0) 1:1 10:1 100:1

Pointwise smoothers do not smooth residuals in all directions of anisotropic problems.

Typical width:height aspect ratios for ice sheet discretizations are on the order of 100:1.

Example: µ = β = 1, Q3 preconditioned by Q1, ILU(0), ordered nodes

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 36 / 58

slide-50
SLIDE 50

Motivation Bayesian Inference Stokes Solvers References

INCOMPLETE FACTORIZATION SMOOTHING

SUMMARY OF GEOMETRIC MULTIGRID THEORY

(2d + 1)-point stencil Laplacian: if the dofs are ordered first in the strong direction, GMG smoothed by lexically ordered ILU on each level converges independent of ε. The ordering is key: tridiagonal solves on tightly coupled subsystems. 1 ε

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 37 / 58

slide-51
SLIDE 51

Motivation Bayesian Inference Stokes Solvers References

INCOMPLETE FACTORIZATION SMOOTHING

IN PRACTICE WITH STANDARD SA-AMG

u · n = 0; (I − n · n)(σ + βu) = 0 Convergence factor, standard SA-AMG (randomized MIS) ℓ β = 10 1−2 1−4 1−8 0.14 0.14 0.57 0.63 1 0.17 0.27 0.75 0.78 2 0.20 0.51 0.82 0.83

(1,1)-block of Stokes system, ℓ levels of refinement, column-ordered dofs on the fine grid

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 38 / 58

slide-52
SLIDE 52

Motivation Bayesian Inference Stokes Solvers References

SMOOTHER/AGGREGATE INCOMPATIBILITY

Strong β masks inefficient smoothing: without it, smoother/aggregate incompatibility is exposed. Column structure is not preserved. Discretization does not neatly separate into tightly coupled groups. “Smoothed” aggregates actually introduce a lot of high-frequency error. ILU smoother is no longer effective on coarse grids.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 39 / 58

slide-53
SLIDE 53

Motivation Bayesian Inference Stokes Solvers References

SMOOTHER/AGGREGATE INCOMPATIBILITY

Strong β masks inefficient smoothing: without it, smoother/aggregate incompatibility is exposed. Column structure is not preserved. Discretization does not neatly separate into tightly coupled groups. “Smoothed” aggregates actually introduce a lot of high-frequency error. ILU smoother is no longer effective on coarse grids.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 39 / 58

slide-54
SLIDE 54

Motivation Bayesian Inference Stokes Solvers References

SMOOTHER/AGGREGATE INCOMPATIBILITY

Strong β masks inefficient smoothing: without it, smoother/aggregate incompatibility is exposed. Column structure is not preserved. Discretization does not neatly separate into tightly coupled groups. “Smoothed” aggregates actually introduce a lot of high-frequency error. ILU smoother is no longer effective on coarse grids.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 39 / 58

slide-55
SLIDE 55

Motivation Bayesian Inference Stokes Solvers References

SMOOTHER/AGGREGATE INCOMPATIBILITY

Strong β masks inefficient smoothing: without it, smoother/aggregate incompatibility is exposed. Column structure is not preserved. Discretization does not neatly separate into tightly coupled groups. “Smoothed” aggregates actually introduce a lot of high-frequency error. ILU smoother is no longer effective on coarse grids.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 39 / 58

slide-56
SLIDE 56

Motivation Bayesian Inference Stokes Solvers References

CONSTRUCTING COLUMN-PRESERVING AGGREGATES: DofColumns

Plugin1 for PETSc’s GAMG preconditioner: Describe column structure of matrix A with MatSetDofColumns(A,...) Add -pc_gamg_type dofcol to command line

Creates aggregate columns, partitions columns to make node aggregates, preserving column structure

1bitbucket.org/tisaac/DofColumns

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 40 / 58

slide-57
SLIDE 57

Motivation Bayesian Inference Stokes Solvers References

INCOMPLETE FACTORIZATION SMOOTHING

IN PRACTICE WITH COLUMN PRESERVING SA-AMG

Convergence factor, standard SA-AMG (3D randomized MIS) ℓ β = 10 1−2 1−4 1−8 0.14 0.14 0.57 0.63 1 0.17 0.27 0.75 0.78 2 0.20 0.51 0.82 0.83 Convergence factor, DofColumns SA-AMG (2D randomized MIS) ℓ β = 10 1−2 1−4 1−8 0.14 0.14 0.30 0.35 1 0.17 0.20 0.39 0.47 2 0.19 0.25 0.48 0.54

(1,1)-block system, ℓ levels of refinement, column-ordered dofs on the fine grid

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 41 / 58

slide-58
SLIDE 58

Motivation Bayesian Inference Stokes Solvers References

INCOMPLETE FACTORIZATION SMOOTHING

DIRECT COMPARISON WITH SSOR ON ANTARCTIC ICE SHEET

SSOR, standard SA-AMG IC(0), column-preserving SA-AMG #Ndof P #N #K #N #K 38M 1,024 8 147 7 85 268M 8,192 9 243 8 98 Comparison for solve to 10−12 of nonlinear Stokes problem, described in [Isaac et al., 2015b], on the Antarctic ice sheet, on the same mesh.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 42 / 58

slide-59
SLIDE 59

Motivation Bayesian Inference Stokes Solvers References

FULL NONLINEAR SOLVER FOR ANTARCTICA

P Solve eff. #N #K Op eff. PC eff. Setup eff. ℓ = 0 uniform refinement, Ndof =7M 128 67.7 1.00 7 66 0.0275 1.00 0.981 1.00 7.361 1.00 256 36.9 0.91 7 67 0.0146 0.94 0.577 0.93 4.255 0.86 512 20.3 0.83 7 65 0.0080 0.86 0.299 0.82 2.722 0.67 ℓ = 1 uniform refinement, Ndof =51M 1,024 78.0 0.86 11 75 0.0276 0.99 0.987 0.99 8.657 0.85 2,048 44.2 0.76 10 75 0.0152 0.90 0.561 0.87 5.810 0.63 4,096 30.2 0.56 10 75 0.0087 0.79 0.383 0.64 5.406 0.34 ℓ = 2 uniform refinement, Ndof =383M 8,192 108.0 0.62 13 91 0.0295 0.93 1.13 0.87 15.58 0.47 16,384 74.7 0.45 10 89 0.0168 0.82 0.80 0.61 17.22 0.21 (from [Isaac et al., 2015c].)

Note the (AMG) preconditioner setup is the least scalable component: Requires sparse matrix triple product for Galerkin projection of course grid operators.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 43 / 58

slide-60
SLIDE 60

Motivation Bayesian Inference Stokes Solvers References

MANTLE CONVECTION: ATTRIBUTION

The following slides summarize “An Extreme-Scale Implicit Solver for Complex PDEs” (SC15 Gordon Bell Prize): Johann Rudi contributed most of the geometric multigrid work, both for the (1,1)-block and the Schur-complement preconditioner. I contributed the highly-scalable [Isaac et al., 2015a] FEM framework common with the Antarctic ice sheet model, the AMG coarse solve, and components of the Schur-complement preconditioner. the IBM Zürich research team contributed kernel and BG/Q-specific

  • ptimizations.
  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 44 / 58

slide-61
SLIDE 61

Motivation Bayesian Inference Stokes Solvers References

MANTLE CONVECTION PRECONDITIONING CHOICES

PARALLEL ADAPTIVE HIGH-ORDER SPECTRAL-GEOMETRIC MULTIGRID

We use hybrid spectral-geometric-algebraic multigrid to precondition ˜ F:

The multigrid hierarchy of nested meshes is generated from an adaptively refined octree-based mesh via spectral-geometric coarsening Parallel repartitioning of coarser meshes for load-balancing High-order L2-projection of fields onto coarser levels Re-discretization of PDEs at coarser geometric multigrid levels

HMG hierarchy pressure space spectral p-coarsening geometric h-coarsening algebraic coars.

  • discont. modal
  • cont. nodal

high-order F.E. trilinear F.E. decreasing #cores #cores < 1000 small MPI communicator single core HMG V-cycle p-MG h-MG AMG direct modal to nodal proj. high-order L2-projection linear L2-projection linear projection

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 45 / 58

slide-62
SLIDE 62

Motivation Bayesian Inference Stokes Solvers References

MANTLE CONVECTION PRECONDITIONING CHOICES

PARALLEL ADAPTIVE HIGH-ORDER SPECTRAL-GEOMETRIC MULTIGRID

We use hybrid spectral-geometric-algebraic multigrid to precondition ˜ F:

Chebyshev accelerated Jacobi smoother (PETSc) with tensorized matrix-free high-order stiffness apply; assembly of high-order diagonal only Restriction & interpolation operators are adjoints of each other in L2-sense. No collective communication needed in spectral-geometric MG cycles Coarse grid solver: AMG (PETSc’s GAMG), invoked on small core counts

HMG hierarchy pressure space spectral p-coarsening geometric h-coarsening algebraic coars.

  • discont. modal
  • cont. nodal

high-order F.E. trilinear F.E. decreasing #cores #cores < 1000 small MPI communicator single core HMG V-cycle p-MG h-MG AMG direct modal to nodal proj. high-order L2-projection linear L2-projection linear projection

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 45 / 58

slide-63
SLIDE 63

Motivation Bayesian Inference Stokes Solvers References

˜ S: BFBT/LEAST SQUARES COMMUTATOR (LSC) PRECONDITIONER

With high-viscosity contrasts (6 orders of magnitude) and narrow weak-zones, −µ−1M ceases to be a good approximation of S = −BF −1B∗. Basic BFBT idea [Elman]: use pseudo inverses for B and B∗, (−BF −1B∗)−1 ≈ −B+FB+, [B∗+ = (B∗B)−1B] Improved by a diagonal symmetric rescaling of equations by D−1 = diag(F)−1 [May and Moresi, 2008]:

(−BF −1B∗)−1 ≈ −B+D−1/2FD−1/2B∗+, [B+ = (B∗D−1B)−1BD−1/2] ≈ −(B∗D−1B)−1BD−1FD−1B∗(B∗D−1B)−1

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 46 / 58

slide-64
SLIDE 64

Motivation Bayesian Inference Stokes Solvers References

REDISCRETIZED LSC

Standard approach: use AMG to approximately apply (B∗D−1B)−1, because it cannot be expressed as an FEM computation. However, if q is smooth, q∗Bv = (q, ∇ · v)Ω = −(∇q, v)Ω + (boundary terms). So B∗D−1B ≈ ∇ · (D−1∇).

We leverage our GMG framework to approximately invert this BVP to approximate (B∗D−1B)−1.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 47 / 58

slide-65
SLIDE 65

Motivation Bayesian Inference Stokes Solvers References

MANTLE CONVECTION PRECONDITIONING COMPARISON

GMRES iteration

100 200 300

Residual reduction

10 -6 10 -4 10 -2 10 0

GMRES restart

GMG_BFBT_4.0km GMG_BFBT_1.5km GMG_BFBT_0.7km AMG_mass_4.0km AMG_mass_1.5km AMG_mass_0.7km

Comparison between “inverse viscosity” mass matrix Schur complement approximation and our improved BFBT preconditioner as weak zone narrows from 4 km to 0.7 km. Mass-matrix Schur complement approximation fails for problems with strong viscosity variations BFBT approximation is much more complicated than mass matrix, but pays off in 4–5 orders of magnitude improvement in convergence

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 48 / 58

slide-66
SLIDE 66

Motivation Bayesian Inference Stokes Solvers References

SOLVER ROBUSTNESS W.R.T. PLATE BOUNDARY THICKNESS

Plate boundary thickness (km) DOF (at solution) # Newton iterations total # GMRES iterations 15 1.00B 25 3915 10 1.63B 27 4645 5 4.78B 29 6112

1 multigrid V-cycle with 3+3 smoothing for both ˜ A and (BD−1B⊤) discretization fixed at Q2 × Pdisc

1

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 49 / 58

slide-67
SLIDE 67

Motivation Bayesian Inference Stokes Solvers References

ALGORITHMIC SCALABILITY WITH DECREASING h

Max level of Finest resolution DOF Newton total GMRES refinement [m] [×106] iterations iterations 10 2443 0.96 14 1408 11 1222 2.67 18 1160 12 611 5.58 21 1185 13 305 11.82 21 1368 14 153 36.35 27 1527

1 multigrid V-cycle with 3+3 smoothing for both ˜ A and (BD−1B⊤) discretization fixed at Q2 × Pdisc

1

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 50 / 58

slide-68
SLIDE 68

Motivation Bayesian Inference Stokes Solvers References

IMPLEMENTATION OPTIMIZATIONS FOR BG/Q

(A) Before optimizations (B) Reduction of blocking MPI communication (C) Minimization of integer

  • perations & cache misses

(D) Optimization of element-local derivatives; SIMD vectorization

Optimization phase A B C D E F G H GFlops/s per node 2.5 5 7.5 10 Normalized time [s] 10 0 10 1 10 2 10 3 10 4

GFlops/s Time

(E) OpenMP threading of matrix-free apply loops (e.g. multigrid smoothing, intergrid projection) (F) MPI communication reduction, overlapping with computations, OpenMP threading in intergrid operators (G) Finite element kernel optimizations (e.g. increase of flop-byte ratio, consecutive memory access, pipelining) (H) Low-level optimizations (e.g. boundary condition enforcement, interpolation of hanging finite element nodes)

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 51 / 58

slide-69
SLIDE 69

Motivation Bayesian Inference Stokes Solvers References

WEAK SCALING ON SEQUOIA: SOLVER AND FULL CODE

Number of Blue Gene/Q cores 16,384 32,768 65,536 131,072 262,144 524,288 1,048,576 Billions DOF / s per GMRES it. 10 -1 10 0 10 1 10 2 1,572,864

0.98 0.99 1.03 1.03 1.03 0.98 0.97 0.98 0.99 1.03 1.03 1.03 0.98 0.96 Solver scalability Full code scalability Ideal scalability

Vulcan Sequoia

Performance normalized by time and number of GMRES iterations Numbers indicate parallel efficiency w.r.t. ideal speedup (baseline = 1 rack) Red indicates linear solver (10 iterations) only Green indicates projected total runtime (includes measured I/O and setup time and estimate of total solver iterations to convergence) Largest problem size has 602 billion DOF on 96 racks

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 52 / 58

slide-70
SLIDE 70

Motivation Bayesian Inference Stokes Solvers References

STRONG SCALING ON SEQUOIA: SOLVER AND FULL CODE

Number of Blue Gene/Q cores 16,384 32,768 65,536 131,072 262,144 524,288 1,048,576 Speedup 1 2 4 8 16 32 64 96 1,572,864

1.03 1.00 0.88 0.76 0.61 0.43 0.33 1.03 1.00 0.88 0.76 0.61 0.43 0.32 Solver speedup Full code speedup Ideal speedup

Vulcan Sequoia

Performance normalized by time and number of GMRES iterations Numbers indicate parallel efficiency w.r.t. ideal speedup (baseline = 1 rack) Red indicates linear solver (10 iterations) only Green indicates projected total runtime (includes measured I/O and setup time and estimate of total solver iterations to convergence) Problem size fixed at 8.3 billion DOF

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 53 / 58

slide-71
SLIDE 71

Motivation Bayesian Inference Stokes Solvers References

STOKES SOLVERS SUMMARY

All though they are both nonlinear Stokes equations with similar rheologies, the solver demands of ice sheet dynamics and mantle convection are quite different: The primary challenge of ice sheets dynamics is efficiently dealing with the anisotropy of the discretization, which must be accounted for when designing smoothers for multigrid V-cycles. The primary challenges of mantle convection are viscosity contrast, which requires robust Schur-complement preconditioning, and shear scale, which must be addressed with both algorithmic scalability and parallel scalability.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 54 / 58

slide-72
SLIDE 72

Motivation Bayesian Inference Stokes Solvers References

FURTHER READING I

Carsten Burstedde, Georg Stadler, Laura Alisic, Lucas C. Wilcox, Eh Tan, Michael Gurnis, and Omar Ghattas. Large-scale adaptive mantle convection simulation. Geophysical Journal International, 192(3): 889–906, 2013. doi: 10.1093/gji/ggs070. Timothy T Creyts and Christian G Schoof. Drainage through subglacial water sheets. Journal of Geophysical Research, 114(F4):F04008, 2009. Nathan Halko, Per Gunnar Martinsson, and Joel A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217–288, 2011. Vincent Heuveline and Friedhelm Schieweck. On the inf-sup condition for higher order mixed fem on meshes with hanging nodes. ESAIM: Mathematical Modelling and Numerical Analysis, 41(01):1–20, 2007.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 55 / 58

slide-73
SLIDE 73

Motivation Bayesian Inference Stokes Solvers References

FURTHER READING II

Tobin Isaac, Carsten Burstedde, Lucas C. Wilcox, and Omar Ghattas. Recursive algorithms for distributed forests of octrees. SIAM Journal on Scientific Computing, 37(5), September 2015a. doi: doi:10.1137/140970963. http://dx.doi.org/10.1137/140970963. Tobin Isaac, Noemi Petra, Georg Stadler, and Omar Ghattas. Scalable and efficient algorithms for the propagation of uncertainty from data through inference to prediction for large-scale problems, with application to flow

  • f the Antarctic ice sheet. Journal of Computational Physics, 296:

348–368, September 2015b. doi: doi:10.1016/j.jcp.2015.04.047. http://dx.doi.org/10.1016/j.jcp.2015.04.047. Tobin Isaac, Georg Stadler, and Omar Ghattas. Solution of nonlinear Stokes equations discretized by high-order finite elements on nonconforming and anisotropic meshes, with application to ice sheet

  • dynamics. SIAM Journal on Scientific Computing (to appear), 2015c.

http://arxiv.org/abs/1406.6573.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 56 / 58

slide-74
SLIDE 74

Motivation Bayesian Inference Stokes Solvers References

FURTHER READING III

C Ketelsen, R Scheichl, and AL Teckentrup. A hierarchical multilevel markov chain monte carlo algorithm with applications to uncertainty quantification in subsurface flow. arXiv preprint arXiv:1303.7343, 2013. Dave A. May and Louis Moresi. Preconditioned iterative methods for Stokes flow problems arising in computational geodynamics. Physics of the Earth and Planetary Interiors, 171:33–47, 2008. Tarek A. El Moselhy and Yousssef M. Marzouk. Bayesian inference with

  • ptimal maps. Journal of Computational Physics, 231(23):7815–7850,

2012. Albert Parker and Colin Fox. Sampling Gaussian distributions in Krylov spaces with conjugate gradients. SIAM Journal on Scientific Computing, 34(3):B312–B334, 2012. E Rignot, J Mouginot, and B Scheuchl. Ice flow of the Antarctic ice sheet. Science, 333(6048):1427–1430, 2011a.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 57 / 58

slide-75
SLIDE 75

Motivation Bayesian Inference Stokes Solvers References

FURTHER READING IV

E Rignot, J Mouginot, and B Scheuchl. MEaSUREs InSAR-Based Antarctica Ice Velocity Map. 2011b.

  • A. Toselli and C. Schwab. Mixed hp-finite element approximations on

geometric edge and boundary layer meshes in three dimensions. Numerische Mathematik, 94(4):771–801, 2003.

  • T. Isaac (U. Chicago)

Scalable Forward & Inverse Solvers January 11, 2016 58 / 58