Implementation of the Jacobian-free Newton-Krylov method for solving the first-order ice sheet momentum balance

  • Authors:

  • Jean-FrançOis Lemieux;Stephen F. Price;Katherine J. Evans;Dana Knoll;Andrew G. Salinger;David M. Holland;Antony J. Payne

  • Affiliations:

  • Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY 10012-1185, USA;Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, USA;Oak Ridge National Laboratory, 1 Bethel Valley Road, P.O. Box 2008, Oak Ridge, TN 37831, USA;Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, USA;Sandia National Laboratories, New Mexico P.O. Box 5800, Albuquerque, NM 87185, USA;Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY 10012-1185, USA;School of Geophysical Sciences, University of Bristol, Bristol, UK

  • Venue:
  • Journal of Computational Physics
  • Year:
  • 2011

Quantified Score

Hi-index 31.45

Visualization

Abstract

We have implemented the Jacobian-free Newton-Krylov (JFNK) method for solving the first-order ice sheet momentum equation in order to improve the numerical performance of the Glimmer-Community Ice Sheet Model (Glimmer-CISM), the land ice component of the Community Earth System Model (CESM). Our JFNK implementation is based on significant re-use of existing code. For example, our physics-based preconditioner uses the original Picard linear solver in Glimmer-CISM. For several test cases spanning a range of geometries and boundary conditions, our JFNK implementation is 1.8-3.6 times more efficient than the standard Picard solver in Glimmer-CISM. Importantly, this computational gain of JFNK over the Picard solver increases when refining the grid. Global convergence of the JFNK solver has been significantly improved by rescaling the equation for the basal boundary condition and through the use of an inexact Newton method. While a diverse set of test cases show that our JFNK implementation is usually robust, for some problems it may fail to converge with increasing resolution (as does the Picard solver). Globalization through parameter continuation did not remedy this problem and future work to improve robustness will explore a combination of Picard and JFNK and the use of homotopy methods.