Statistical Inversion for Basal Parameters for the Antarctic Ice - PowerPoint PPT Presentation

Statistical Inversion for Basal Parameters for the Antarctic Ice Sheet Tobin Isaac, Noémi Petra, Georg Stadler, Omar Ghattas 1 Institute for Computational Engineering & Sciences The University of Texas at Austin SIAM UQ 2014 Savannah, Georgia April 3, 2014 1 Departments of Geological Sciences and Mechanical Engineering Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 1 / 28

1 Forward Model: Ice Sheet Dynamics 2 Bayesian Inversion via methods of Constrained Optimization Example: Pine Island Glacier 3 Example: Antarctic Ice Sheet 4 These slides can be found at www.ices.utexas.edu/~tisaac/slides. Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 2 / 28

Motivation Antarctic ice sheet inversion for the basal friction parameter field using InSAR surface velocity measurements Observed surface flow velocity [Rignot et al., 2011] IPCC Fourth Assessment Report [Meehl et al., 2007] [T]he uncertainty in the projections of the land ice contributions [to sea level rise] is dominated by the various uncertainties in the land ice models themselves ... rather than in the temperature projections.’’ Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 3 / 28

Equations Balance of linear momentum, mass, and energy ε = 1 2 ( ∇ u + ∇ u T )] − ∇ · [ η ( θ, u ) ˙ ε − I p ] = ρ g , [˙ ∇ · u = 0 , � ∂θ � ε 2 ) ρ c ∂ t + u · ∇ θ − ∇ · ( K ∇ θ ) = 2 η tr (˙ Constitutive relations �� − 1 � � − Q 1 − n n ε II = 1 ε 2 )] ˙ [˙ 2 tr (˙ η ( θ, u ) = ε 2 n A 0 exp II R θ Basal boundary conditions T ( σ n ) + β ( x ) Tu = 0 Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 4 / 28

Discretization High-order finite element discretization: Q p × Q p − 2 complex geometry high-order accurate 3 4 LBB stable 0 element-wise mass conservative 1 2 inf-sup constant independent of h and weakly dependent on p inf-sup constant unaffected by high aspect ratio ice sheet geometry amenable to scalable parallel AMR schemes better efficiency than low-order (increased flops per memop), if there is an efficient solver. Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 5 / 28

Solution Nonlinear solver is inexact Newton-Krylov method: Newton’s method: linearized Stokes equations have 4th-order tensor coefficient, as do true adjoint equations. Inexact Krylov tolerances [Eisenstat and Walker, 1996]: minimizes Krylov iterations in Newton steps not believed to be in asymptotic regime. Linear solver: trilinear submesh preconditioner for high-order elements [Brown, 2010] Smoothed Aggregation AMG: tolerates heterogeneities/anisotropies of viscous operator better than geometric multigrid ASM( i )/ILU( j ) smoother for high-aspect-ratio geometries Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 6 / 28

Antarctic simulation: scaling Newton-Krylov Solver for Antarctic Ice Sheet Stokes Equations 10 − 1 p = 3 [18] p = 4 [23] 10 − 4 p = 5 [23] � r k � / � r 0 � 10 − 7 10 − 10 10 − 13 0 50 100 150 200 Krylov iterations Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 7 / 28

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 bC bC bC Antarctic simulation: scaling Newton-Krylov Solver for Antarctic Ice Sheet Stokes Equations 10 0 7 35 16 5 10 − 2 33 7 40 42 4 5 3 3 4 30 17 7 4 10 − 4 4 14 33 38 6 21 || r k || / || r 0 || 5 93 10 − 6 54 8 3 13 10 − 8 8 196 75 10 − 10 82 28M dofs, 1024 cores 143 122M dofs, 4096 cores 10 − 12 122M dofs, 16384 cores 356 104 10 − 14 0 3 6 9 12 15 Newton iterations (Krylov iterations labeled) Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 7 / 28

Inventory We have Scalable iterative solution of nonlinear forward problem Scalable iterative solution of Newton linearization Scalable iterative solution of adjoint equations Time to solution is more-or-less linear in the tolerance Algebraic multilevel representation of linear operator We do not currently have Full multilevel approximation scheme of the nonlinear problem Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 8 / 28

Bayesian-inversion-in-a-slide Let f be the map from the parameter m (controlling the basal coefficient β ) to the observations through the forward problem: f ( m ) = B ( u ) , A ( u , m ) = 0 We assume additive Gaussian noise in the measurements: d = f ( m ) + e , e ∼ N ( 0 , Γ noise ) Thus: 2 ( f ( m ) − d ) T Γ − 1 � − 1 � π like ( d | m ) = exp noise ( f ( m ) − d ) If the prior is Gaussian with mean m prior and covariance Γ prior , then we obtain for the posterior pdf: � � − 1 2 � f ( m ) − d � 2 noise − 1 2 � m − m prior � 2 π post ( m ) ∝ exp Γ − 1 Γ − 1 prior Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 9 / 28

Posterior statistic 0: MAP point The maximum a posteriori point is def m MAP = arg max π post ( m ) m def = arg min − log π post ( m ) [ = J ( m )] m 2 � B ( u ) − d � 2 2 � m − m prior � 2 1 noise + 1 = arg min Γ − 1 Γ − 1 m , u prior such that A ( u , m ) = 0 . 2 Γ − 1 This is classical PDE-constrained optimization: 1 prior characterizes the regularization. Γ − 1 prior is (should be) given a priori by expert knowledge, not chosen a posteriori when the solution has desirable regularity (although if ‘‘expert prior knowledge’’ is ‘‘I know it when I see it’’, there is little difference) Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 10 / 28

Scalable MAP point computation via RSINK The KKT first-order necessary conditions ( J ′ ( m ) = 0 ) of a solution m are A ( u , m ) = 0 [state] , A T u ( u , m ) p = B ( u ) − d [adjoint] , m ( u , m ) p + Γ − 1 A T prior m = 0 [gradient] , where u and p are state and adjoint variables. We solve this using R educed S pace I nexact N ewton- K rylov: Reduced space: search for a minimum of J ( m ) ; given m k , u k and p k satisfy state and adjoint equations. Inexact Newton-Krylov: like inexact Newton-Krylov used to solve state equations J ′′ ( m k ) δ k = − J ′ ( m k ) m k + 1 = m k + α δ k , where α is determined by line search. Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 11 / 28

Necessary components of scalable, efficient RSINK implementation Scalable, efficient cost evaluation ⇔ Scalable, efficient nonlinear forward solver Scalable, efficient gradient computation Scalable, efficient (Gauss-)Newton Hessian application def = J ′′ ( m ) . H ( m ) The Hessian, or an approximation, is useful for UQ: please attend Noémi’s talk in MS7 three days ago for more details. Scalable, efficient (Gauss-)Newton Hessian preconditioner Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 12 / 28

Adjoint-methods-in-a-slide If u k and p k solve state and adjoint equations ( u k computed from m k , p k computed form u k and m k ), then J ′ ( m ) = p T k Γ − 1 k A m ( u k , m k ) + m T prior All terms in equations have finite element expressions (not differentiating matrices). Marginal cost of J ′ ( m ) , given J ( m ) : one linear PDE solve. Linearizing the first-order KKT system of equations results in set of PDEs for computing Hessian action Computing H ( m k ) ˜ m costs two (‘‘incremental’’) linear PDE solves with the same operators as Newton step / adjoint. Another opportunity to exploit inexactness: loose tolerance on incrementals approximates H ( m k ) , but quadratic/superlinear convergence of inexact Newton only requires accurate Hessian near the solution [Hicken and Alonso, 2013]. Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 13 / 28

Prior preconditioning the Hessian misfit + Γ − 1 We precondition the full Hessian H = H prior with Γ prior , so the preconditioned system is 1 / 2 1 / 2 H prior + I ∼ Γ prior H prior + I . misfit Γ misfit Γ The justification from PDE-constrained optimization: The precision Γ − 1 prior is a differential operator. The prior Γ prior is compact (low effective rank), so H misfit Γ prior is compact. Identity + compact ⇒ good ( h -independent) Krylov method convergence. The implication for statistical inversion: post of Gaussianized posterior distribution is H ( m MAP ) − 1 . Covariance Γ The effective numerical rank of H misfit Γ prior is the dimension of the subspace of parameters informed by the data. The more information the data contains, the less compact is H misfit Γ prior Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 14 / 28