Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton-Krylov method

  • Authors:

  • Jean-François Lemieux;Bruno Tremblay;Jan Sedláček;Paul Tupper;Stephen Thomas;David Huard;Jean-Pierre Auclair

  • Affiliations:

  • Department of Atmospheric and Oceanic Sciences, McGill University, 805 Sherbrooke Street West, Montréal, QC, Canada H3A 2K6;Department of Atmospheric and Oceanic Sciences, McGill University, 805 Sherbrooke Street West, Montréal, QC, Canada H3A 2K6;Institute for Atmospheric and Climate Science, ETH Zurich, Universitaetstrasse 16, Zurich CH-8092, Switzerland;Department of Mathematics, Simon Fraser University, 8888 University Drive, Burnaby, BC, Canada V5A 1S6;Acceleware Corp., 1600 37th Street SW, Calgary, AB, Canada T3C 3P1;Department of Atmospheric and Oceanic Sciences, McGill University, 805 Sherbrooke Street West, Montréal, QC, Canada H3A 2K6;Department of Oceanography, Dalhousie University, 1355 Oxford Street, Halifax, NS, Canada B3H 4J1

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

Quantified Score

Hi-index 31.48

Visualization

Abstract

We have implemented the Jacobian-free Newton-Krylov (JFNK) method to solve the sea ice momentum equation with a viscous-plastic (VP) formulation. The JFNK method has many advantages: the system matrix (the Jacobian) does not need to be formed and stored, the method is parallelizable and the convergence can be nearly quadratic in the vicinity of the solution. The convergence rate of our JFNK implementation is characterized by two phases: an initial phase with slow convergence and a fast phase for which the residual norm decreases significantly from one Newton iteration to the next. Because of this fast phase, the computational gain of the JFNK method over the standard solver used in existing VP models increases with the required drop in the residual norm (termination criterion). The JFNK method is between 3 and 6.6 times faster (depending on the spatial resolution and termination criterion) than the standard solver using a preconditioned generalized minimum residual method. Resolutions tested in this study are 80, 40, 20 and 10km. For a large required drop in the residual norm, both JFNK and standard solvers sometimes do not converge. The failure rate for both solvers increases as the grid is refined but stays relatively small (less than 2.3% of failures). With increasing spatial resolution, the velocity gradients (sea ice deformations) get more and more important. Nonlinear solvers such as the JFNK method tend to have difficulties when there are such sharp structures in the solution. This lack of robustness of both solvers is however a debatable problem as it mostly occurs for large required drops in the residual norm. Furthermore, when it occurs, it usually affects only a few grid cells, i.e., the residual is small for all the velocity components except in very localized regions. Globalization approaches for the JFNK solver, such as the line search method, have not yet proven to be successful. Further investigation is needed.