Discretization and Solvers for the Stokes Equations of Ice Sheet - - PowerPoint PPT Presentation

Discretization and Solvers for the Stokes Equations of Ice Sheet Dynamicss

Tobin Isaac

Institute for Computational Engineering and Sciences The University of Texas at Austin

SIAM GS, Padua, June 18, 2013

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 1 / 30

Outline

Motivation Discretization Solver Results

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 2 / 30

Motivation

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 3 / 30

Ice sheet dynamics

Ice sheets exhibit creeping, shear-thinning, incompressible flow. Continental ice sheet flow has

◮ narrow margins between slow and fast flow, ◮ variable viscosity, ◮ complex boundary conditions, ◮ a domain with a large width:height aspect ratio.

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 4 / 30

Our simulation approach

We use PDE-constrained optimization techniques to infer parameters1, and to infer statistical uncertainty in parameters2. These techniques require that we solve nonlinear Stokes equations:

◮ accurately:

◮ model state equations in 3D, ◮ satisfy incompressibility without stabilization, ◮ resolve flow features at appropriate length scales;

◮ efficiently:

◮ adaptive mesh refinement (AMR), ◮ Newton-Krylov method to solve nonlinear equations, ◮ scalable solvers for Newton step;

◮ robustly:

◮ parameter regimes may vary widely during optimization.

We want to bring this all together at continental scale.

• 1N. Petra, H. Zhu, G. Stadler, T. J. R. Hughes, O. Ghattas, Journal of Glaciology 58, 889 (2012)
• 2N. Petra, J. Martin, G. Stadler, O. Ghattas, In preparation (2013)
• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 5 / 30

Discretization

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 6 / 30

Meshing

The meshing pipeline:

◮ describe geometry from NetCDF data provided by glaciologists3, ◮ create coarse hexahedral mesh (postprocessed Triangle4 mesh), ◮ adapt.

• 3A. M. Le Brocq, A. J. Payne, A. Vieli, Earth System Science Data 2, 247 (2010)
• 4J. R. Shewchuk, Proceedings of the Twelfth Annual Symposium on Computational Geometry (Association for

Computing Machinery, 1996), pp. 141–150

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 7 / 30

Adaptive Mesh Refinement with p4est

0.1 1 10 100 12 24 48 96 192 384 768 1536 3072 6144 Seconds Number of CPU cores Perfect Scaling Old New 1 2 3 4 5 6 12 96 768 6144 49152 112128 Seconds per (million elements / core) Number of CPU cores Old New

Parallel AMR provided by p4est library5:

◮ forest-of-octrees allow arbitrary hierarchical refinement of complicated

geometries,

◮ fully distributed fine mesh organized by space-filling curve, ◮ highly scalable algorithms6 7.

• 5C. Burstedde, p4est: Parallel AMR on forests of octrees (2010). http://www.p4est.org/
• 6C. Burstedde, et al., SC10: Proceedings of the International Conference for High Performance Computing,

Networking, Storage and Analysis (ACM/IEEE, 2010)

• 7T. Isaac, C. Burstedde, O. Ghattas, Proceedings of the 26th IEEE International Parallel & Distributed

Processing Symposium (IEEE, 2012)

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 8 / 30

High-order FEM

Primary space:

◮ Qk hexahedral elements ◮ Fast matrix-free matvec8

Incompressibility constraint space: discontinuous spaces that satisfy incompressibility on each element

◮ Qk−2:

◮ inf-sup constant bounded

• ver h, weak k-dependence

◮ inf-sup constant unaffected

by boundary-layer refinement 9

• 8J. Brown, Journal of Scientific Computing 45, 48 (2010)
• 9V. Heuveline, F. Schieweck, ESAIM: Mathematical Modelling and Numerical Analysis 41, 1 (2007)
• 9A. Toselli, C. Schwab, Numerische Mathematik 94, 771 (2003)
• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 9 / 30

Discretized, Linearized PDEs

Γbase Ω

We use an inexact Newton’s method to solve nonlinear Stokes and for

• ptimization: our linear solver must handle the linearized equations

(Newton step) and their adjoint equations. Both have the variational form {(Dv, µ′(u)D˜ u)Ω + (T

v, βT ˜

u)Γbase} −(∇ · v, ˜ p)Ω −(q, ∇ · ˜ u)Ω

• =

ru,p(v) ru,p(q)

• ∀(v, q).

◮ The effective viscosity 2µ′(u) is a 4th-order tensor and varies over

several orders of magnitude.

◮ The effective friction β also varies over several orders of magnitude.

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 10 / 30

Solver

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 11 / 30

Preconditioning the indefinite system

We precondition the Stokes matrix A with ˜ A: A = F BT B

• ;

˜ A = ˜ F BT − ˜ S

• .

◮ Converges in 2 iterations of ˜

F = F and ˜ S = BF −1BT.

◮ Choices for ˜

S:

◮ Scaled lumped mass matrix

1 |µ| ˆ

M,

◮ Least squares commutator (BBT)−1BFBT (BBT)−1,

(BBT)−1 computed approximately10.

• 10H. C. Elman, V. E. Howle, J. Shadid, D. J. Silvester, R. Tuminaro, SIAM Journal on Scientific Computing 30,

290 (2007)

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 12 / 30

Momentum block preconditioner

˜ A = ˜ F BT − ˜ S

• We would like the convergence rate to be independent of

◮ h, ◮ k, ◮ µ′, ◮ β, ◮ the mesh (nonconforming interfaces), ◮ φ = H/L (δ = L/H).

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 13 / 30

h-independence: algebraic multigrid

◮ A multilevel method for the momentum block is necessary for

h-independence.

◮ We use existing smoothed aggregation AMG libraries (Trilinos’s ML,

PETSc’s GAMG).

◮ There is considerable flexibility

◮ hierarchy must complement ◮ smoothing must complement ◮ discretization.

◮ Not µ′-independence, but more µ′-robust than geometric multigrid.

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 14 / 30

φ-independence (φ = H/L)

10 20 30 40 50 60

Krylov iterations

10−15 10−12 10−9 10−6 10−3 100

rk/r0

Gauss-Seidel smoothing φ = 1 φ = 10 φ = 100

◮ Pointwise smoothers (Jacobi, GS) are φ-dependent.

◮ 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 (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 15 / 30

φ-independence (φ = H/L)

10 20 30 40 50 60

Krylov iterations

10−15 10−12 10−9 10−6 10−3 100

rk/r0

Default parameters, k = 3 φ = 1 φ = 10 φ = 100

◮ Pointwise smoothers (Jacobi, GS) are φ-dependent.

◮ Typical width:height aspect ratios for ice sheet discretizations are on

the order of 100:1.

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

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 15 / 30

k-independence: approximating higher order-elements

Element matrix FK is

◮ dense (O(n2 e) entries), ◮ expensive (O(n7/3 e

) steps), (ne = (k + 1)3) We treat the nodes as trilinear nodes, create ˜ FK:

◮ sparse (O(ne) entries), ◮ cheap (O(ne) steps), ◮ spectrally equivalent (proven

in 2D)11.

• 11S. Kim, Electronic Transactions on Numerical Analysis 26, 228 (2007)
• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 16 / 30

k-independence: approximating higher order-elements

10 20 30 40 50 60

Krylov iterations

10−15 10−12 10−9 10−6 10−3 100

rk/r0

φ = 100, varying k k = 3 k = 6 k = 9

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 17 / 30

β-dependence

◮ For large β, the same convergence rate is seen for all φ ◮ For small β, the convergence is φ-dependent

10 20 30 40 50 60

Krylov iterations

10−15 10−12 10−9 10−6 10−3 100

rk/r0

β ≡ 10−10 φ = 1 φ = 10 φ = 100

◮ For fixed φ, the convergence rate is bounded below as β → 0

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 18 / 30

Mesh independence: nonconforming meshes

10 20 30 40 50 60

Krylov iterations

10−15 10−12 10−9 10−6 10−3 100

rk/r0

Default parameters, k = 3 φ = 1 φ = 10 φ = 100

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 19 / 30

Mesh independence: nonconforming meshes

10 20 30 40 50 60

Krylov iterations

10−15 10−12 10−9 10−6 10−3 100

rk/r0

Txyz mesh { ˜ RK}, φ = 1 { ˜ RK}, φ = 10 { ˜ RK}, φ = 100 {RK}, φ = 1 {RK}, φ = 10 {RK}, φ = 100

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 20 / 30

Results

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 21 / 30

ISMIP-HOM Benchmark C, stream variation

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 22 / 30

ISMIP-HOM Benchmark C, stream variation

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 23 / 30

ISMIP-HOM Benchmark C, stream variation

10 20 30 40 50 60 70 80 Krylov iterations 10−12 10−11 10−10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 rk/r0

ISMIP-HOM benchmark C, stream 3 4 5

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 24 / 30

Antarctic simulation: varying β

Three nonlinear solves:

◮ full Antarctic ice sheet (no shelves), ◮ vary weakest friction β over three orders of magnitude, ◮ Q3 × Q1 discretization, ◮ moderate refinement at grounding line (∼ (1.25 − 2.5km)2), ◮ “good mesh refinement”, ◮ ∼264K elements, ∼28M variables, ◮ ˜

FS,GL approximate momentum block,

◮ 1 µ ˆ

M approximate Schur complement,

◮ Processor-block additive Schwarz, ILU(1) smoother, ML generated

hierarchy,

◮ adaptive Newton step tolerance12.

• 12S. C. Eisenstat, H. F. Walker, SIAM Journal on Scientific Computing 17, 16 (1996)
• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 25 / 30

Antarctic simulation: varying β

Convergence plots for varying choices of β

2 4 6 8 10 12 10−12 10−10 10−8 10−6 10−4 10−2 100 102 104

Newton iterations (Krylov iterations labeled) ||rk||/||r0|| Newton-Krylov Solver for Antarctic Ice Sheet Stokes Equations

bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC

1011 Pa s / m, min. friction 109 Pa s / m 108 Pa s / m

bC bC bC

16 24 22 13 75 129 467 695 74 100 248 300 78 123 22 15 7 7 4 3 23 15 7 4 28 27 29

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 26 / 30

Antarctic simulation: algorithmic scaling

Simulations conducted on TACC’s Stampede machine:

◮ minimum β 109 Pa s/m, ◮ same solver parameters as before, ◮ two meshes:

◮ mesh from previous example (265K elements, 28M variables), ◮ solved on 1024 cores (∼27K variables/core), ◮ larger mesh: ◮ one additional level uniform refinement, ◮ finer grounding line refinement (∼ (0.6 − 1.2km)2), ◮ ∼1.2M elements, ∼122M variables, ◮ solved on 4096 cores (∼30K variables/core), ◮ solved on 16384 cores (∼7K variables/core).

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 27 / 30

Antarctic simulation: scaling

3 6 9 12 15 10−14 10−12 10−10 10−8 10−6 10−4 10−2 100

Newton iterations (Krylov iterations labeled) ||rk||/||r0|| Newton-Krylov Solver for Antarctic Ice Sheet Stokes Equations

bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC

28M dofs, 1024 cores 122M dofs, 4096 cores 122M dofs, 16384 cores

bC bC bC

30 4 7 14 21 54 8 3 8 82 104 33 5 7 17 33 38 6 5 13 75 143 35 5 7 16 40 42 4 3 3 4 4 93 196 356

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 28 / 30

Antarctic simulation: larger problem

Nonlinear solve on Antarctic ice sheet

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 29 / 30

Antarctic simulation: larger problem

Nonlinear solve on Antarctic ice sheet

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 29 / 30

Antarctic simulation: larger problem

Nonlinear solve on Antarctic ice sheet

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 29 / 30

Thank you!

• T. Isaac (ICES, UT Austin)

Stokes Solvers for Ice June 18, 2013 30 / 30