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 parameters 1 , and to infer statistical uncertainty in parameters 2 . 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 glaciologists 3 , ◮ create coarse hexahedral mesh (postprocessed Triangle 4 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 Perfect Scaling Old New Old New 100 6 5 Seconds per (million elements / core) 10 4 Seconds 3 1 2 1 0.1 0 12 24 48 96 192 384 768 1536 3072 6144 12 96 768 6144 49152 112128 Number of CPU cores Number of CPU cores Parallel AMR provided by p4est library 5 : ◮ forest-of-octrees allow arbitrary hierarchical refinement of complicated geometries, ◮ fully distributed fine mesh organized by space-filling curve, ◮ highly scalable algorithms 6 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 Incompressibility constraint space: discontinuous spaces that satisfy incompressibility on each element ◮ Q k − 2 : ◮ inf-sup constant bounded over h , weak k -dependence ◮ inf-sup constant unaffected Primary space: by boundary-layer refinement 9 ◮ Q k hexahedral elements ◮ Fast matrix-free matvec 8 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 optimization: our linear solver must handle the linearized equations (Newton step) and their adjoint equations. Both have the variational form � { ( D v , µ ′ ( u ) D˜ u ) Ω + (T � v , β T u ) Γ base } − ( ∇ · v , ˜ p ) Ω � � r u ,p ( v ) � � ˜ = ∀ ( v , q ) . − ( q, ∇ · ˜ u ) Ω r u ,p ( 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 : � ˜ B T B T � F � � F ˜ A = ; A = . − ˜ B 0 0 S ◮ Converges in 2 iterations of ˜ F = F and ˜ S = BF − 1 B T . ◮ Choices for ˜ S : | µ | ˆ ◮ Scaled lumped mass matrix 1 M , ◮ Least squares commutator ( BB T ) − 1 BFB T ( BB T ) − 1 , ( BB T ) − 1 computed approximately 10 . 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 � ˜ B T � F ˜ A = − ˜ 0 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 ) Gauss-Seidel smoothing 10 0 φ = 1 10 − 3 φ = 10 φ = 100 � r k � / � r 0 � 10 − 6 10 − 9 10 − 12 10 − 15 0 10 20 30 40 50 60 Krylov iterations ◮ Pointwise smoothers (Jacobi, GS) are φ -dependent. ◮ Typical width:height aspect ratios for ice sheet discretizations are on the order of 100:1. ◮ Example: µ = β = 1 , Q 3 preconditioned by Q 1 , Gauss-Seidel T. Isaac (ICES, UT Austin) Stokes Solvers for Ice June 18, 2013 15 / 30
φ -independence ( φ = H/L ) Default parameters, k = 3 10 0 φ = 1 10 − 3 φ = 10 φ = 100 � r k � / � r 0 � 10 − 6 10 − 9 10 − 12 10 − 15 0 10 20 30 40 50 60 Krylov iterations ◮ Pointwise smoothers (Jacobi, GS) are φ -dependent. ◮ Typical width:height aspect ratios for ice sheet discretizations are on the order of 100:1. ◮ Example: µ = β = 1 , Q 3 preconditioned by Q 1 , ILU(0) T. Isaac (ICES, UT Austin) Stokes Solvers for Ice June 18, 2013 15 / 30
k -independence: approximating higher order-elements We treat the nodes as trilinear Element matrix F K is nodes, create ˜ F K : ◮ sparse ( O ( n e ) entries), ◮ dense ( O ( n 2 e ) entries), ◮ cheap ( O ( n e ) steps), ◮ expensive ( O ( n 7 / 3 ) steps), e ◮ spectrally equivalent (proven ( n e = ( k + 1) 3 ) 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 φ = 100, varying k 10 0 k = 3 10 − 3 k = 6 k = 9 � r k � / � r 0 � 10 − 6 10 − 9 10 − 12 10 − 15 0 10 20 30 40 50 60 Krylov iterations 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 − 10 10 0 φ = 1 10 − 3 φ = 10 φ = 100 � r k � / � r 0 � 10 − 6 10 − 9 10 − 12 10 − 15 0 10 20 30 40 50 60 Krylov iterations ◮ 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 Default parameters, k = 3 10 0 φ = 1 10 − 3 φ = 10 φ = 100 � r k � / � r 0 � 10 − 6 10 − 9 10 − 12 10 − 15 0 10 20 30 40 50 60 Krylov iterations T. Isaac (ICES, UT Austin) Stokes Solvers for Ice June 18, 2013 19 / 30
Mesh independence: nonconforming meshes T xyz mesh 10 0 { ˜ R K } , φ = 1 10 − 3 { ˜ R K } , φ = 10 { ˜ � r k � / � r 0 � R K } , φ = 100 10 − 6 { R K } , φ = 1 10 − 9 { R K } , φ = 10 { R K } , φ = 100 10 − 12 10 − 15 0 10 20 30 40 50 60 Krylov iterations 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
Recommend
More recommend
