Fully implicit timestepping methods for the rotating shallow water equations

TL;DR

The paper addresses time integration for geophysical fluid dynamics by evaluating fully implicit Runge-Kutta (IRK) timestepping for the rotating shallow water equations on the sphere. It integrates a compatible finite element discretisation with a monolithic Newton-Krylov solver preconditioned by overlapping additive Schwarz methods using vertex-star patches, and compares Gauss-Legendre and Radau IIA collocation IRKs (1–3 stages) against the ARK2 IMEX benchmark. Results on the Rossby-Haurwitz test show that Gauss-Legendre IRKs can achieve faster wallclock times than ARK2 at equivalent accuracy for large stable timesteps, with Radau IIA showing similar trends albeit with different accuracy profiles; both exhibit some loss of timestep robustness at very large $\Delta t$. The findings demonstrate the viability of IRK timestepping for 2D atmosphere/ocean-scale problems and point to opportunities for performance tuning and extension to 3D with more scalable solvers.

Abstract

Fully implicit timestepping methods have several potential advantages for atmosphere/ocean simulation. First, being unconditionally stable, they degrade more gracefully as the Courant number increases, typically requiring more solver iterations rather than suddenly blowing up. Second, particular choices of implicit timestepping methods can extend energy conservation properties of spatial discretisations to the fully discrete method. Third, these methods avoid issues related to splitting errors that can occur in some situations, and avoid the complexities of splitting methods. Fully implicit timestepping methods have had limited application in geophysical fluid dynamics due to challenges of finding suitable iterative solvers, since the coupled treatment of advection prevents the standard elimination techniques. However, overlapping Additive Schwarz methods, provide a robust, scalable iterative approach for solving the monolithic coupled system for all fields and Runge-Kutta stages. In this study we investigate this approach applied to the rotating shallow water equations, facilitated by the Irksome package which provides automated code generation for implicit Runge-Kutta methods. We compare various schemes in terms of accuracy and efficiency using an implicit/explicit splitting method, namely the ARK2 scheme of Giraldo et al (2013), as a benchmark. This provides an initial look at whether implicit Runge-Kutta methods can be viable for atmosphere and ocean simulation.

Figures (4)

• Figure 1: A diagram showing the "vertex star" patch associated with the vertex at the middle of the diagram. Dots indicate $D$ degrees of freedom and thin lines indicate $u$ degrees of freedom. Of the latter, lines crossing cell edges indicate normal components of $u$, which are continuous across cell edges for our choice of finite element space. Lines parallel to cell edges indicate tangential components of $u$, which are discontinuous across cell edges. The shaded region indicates the patch. All degrees of freedom outside the shaded region (including normal but not tangential components on the patch boundary) are excluded from the patch.
• Figure 2: Top row: initial vorticity (left column) and depth (right column) fields for Williamson test case 6. Bottom row: both fields at day 1. The shown fields are for mesh level 5, which has 20480 cells, 61440 height degrees of freedom and 92160 velocity degrees of freedom.
• Figure 3: Relative errors of Gauss-Legendre 1-3 and ARK2 IMEX methods vs total wallclock time. The numbers on the curves indicate the time step sizes $\Delta t$ used for the corresponding solution. Solid lines indicate values for the elevation field $\eta=D-H$ and dashed lines for the velocity fields $u$. The solid and dashed vertical lines indicate the spatial discretisation error between mesh levels L6 and L7 of $\eta$ and $u$, respectively.
• Figure 4: Relative errors of Radau IIA 1-3 and ARK2 IMEX methods vs total wallclock time. The numbers on the curves indicate the time step sizes $\Delta t$ used for the corresponding solution. Solid lines indicate values for the elevation field $\eta=D-H$ and dashed lines for the velocity fields $u$. The solid and dashed vertical lines indicate the spatial discretisation error between mesh levels L6 and L7 of $\eta$ and $u$, respectively.