Solution of Nonlinear Stokes Equations Discretized By High-Order Finite Elements on Nonconforming and Anisotropic Meshes, with Application to Ice Sheet Dynamics

Abstract

Abstract. Motivated by the need for efficient and accurate simulation of the dynamics of the polar ice sheets, we design high-order finite element discretizations and scalable solvers for the solution of nonlinear incompressible Stokes equations. In particular, we focus on power-law, shear thinning rheologies commonly used in modeling ice dynamics and other geophysical flows. We use nonconforming hexahedral meshes and the conforming inf-sup stable finite element velocity-pressure pairings, where k ≥ 2 is the polynomial order of the velocity space. To solve the nonlinear equations, we propose a Newton-Krylov method with a block upper triangular preconditioner for the linearized Stokes systems. The diagonal blocks of this preconditioner are sparse approximations of the (1,1)-block and of its Schur complement. The (1,1)-block is approximated using linear finite elements based on the nodes of the high-order discretization, and the application of its inverse is approximated using algebraic multigrid with an incomplete factorization smoother. This preconditioner is designed to be efficient on anisotropic meshes, which are necessary to match the high aspect ratio domains typical for ice sheets. As part of this work, we develop and make available extensions to two libraries-a hybrid meshing scheme for the p4est parallel adaptive mesh refinement library, and a modified smoothed aggregation scheme for PETSc-to improve their support for solving PDEs in high aspect ratio domains. In a comprehensive numerical study, we find that our solver yields fast convergence that is independent of the element aspect ratio, the occurrence of nonconforming interfaces, and of the mesh refinement, and that depends only weakly on the polynomial finite element order. We simulate the ice flow in a realistic description of the Antarctic ice sheet derived from field data, and study the parallel scalability of our solver for problems with up to 383 million unknowns.Key words. Viscous incompressible flow, nonlinear Stokes equations, shear-thinning, high-order finite elements, preconditioning, multigrid, Newton-Krylov method, ice sheet modeling, Antarctic ice sheet.1. Introduction. We design high-order finite element discretizations and scalable solvers for incompressible nonlinear Stokes equations describing creeping flows of power-law rheology fluids. Applications include ice sheet dynamics [31], mantle convection [53], magma dynamics [44] and other problems involving non-Newtonian fluids [26]. Among the main challenges for the solution of these problems are the presence of local features that emerge from the nonlinear constitutive relation, the strongly varying and anisotropic coefficients arising upon linearization, the incompressibility condition leading to indefinite matrix problems, complex geometry and boundary conditions, a wide range of length scales that may require highly-adapted meshes with high aspect ratios, and large problem sizes that necessitate parallel solution on large supercomputers. Our approach to cope with these challenges ...

<i>Albany/FELIX</i>: a parallel, scalable and robust, finite element, first-order Stokes approximation ice sheet solver built for advanced analysis

“…For this process to be computationally tractable during both the inverse (parameter estimation and uncertainty assignment) and forward propagation steps, it is crucial to have robust, efficient, and scalable solves on HPC platforms (Isaac et al, 2014). This, in turn, requires advanced dynamical core capabilities, such as access to model derivatives (e.g., the Jacobian matrix), and advanced algorithms for the solution of the nonlinear and linear equations.…”

“…Oscillatory components are effectively reduced through a simple iterative procedure, while smooth components are tackled using auxiliary lower-resolution versions of the problem. Different geometric multigrid methods have been successfully applied to the linear systems arising from ice sheet modeling simulations, e.g., Brown et al (2013); Cornford et al (2013); Isaac et al (2014).…”

“…At the discrete level, these aspect ratios give rise to matrices where entries representing vertical coupling are generally much larger than entries representing horizontal coupling. Anisotropic phenomena within ice sheets and fairly different types of multigrid methods have been considered in recent prior works (Brown et al, 2013;Isaac et al, 2014;Jouvet and Graser, 2013).…”

Extreme-Scale Multigrid Components within PETSc

“…high aspect ratio elements, see Chapter 7 [28]). Such an approach was advocated in [6,12], where ICC(0) with a column-oriented (e.g. perpendicular to the anisotropy) ordering of the unknowns provides an exact column solve and thus is a highly efficient smoother.…”

Optimal numerical solvers for transient simulations of ice flow using the Ice Sheet System Model (ISSM versions 4.2.5 and 4.11)

“…The finite-element discretization of the full Stokes model leads to a well-studied saddle point problem, which represents an active area of research in geophysics (e.g., Benzi et al, 2005;Elman et al, 2014). While recent work (e.g., Isaac et al, 2015) has shown promising results, stable iterative full Stokes solvers are not readily available and, in general, are significantly disruptive to integrate in terms of their code base, which is the reason we will not be considering them in this study.…”