Massively parallel solvers for elliptic PDEs in Numerical Weather- and Climate Prediction

TL;DR

This paper demonstrates that advanced parallel solvers like Krylov subspace methods and multigrid algorithms can efficiently solve large-scale elliptic PDEs in weather and climate models, enabling high-resolution simulations within operational time constraints.

Contribution

It provides a comprehensive evaluation and optimization of scalable elliptic PDE solvers, including Krylov and multigrid methods, for large atmospheric models on supercomputers.

Findings

Multigrid solver shows robust performance across parameters.

Krylov and multigrid methods scale efficiently to 65,536 cores.

Elliptic solvers can meet operational time scales in high-resolution models.

Abstract

The demand for substantial increases in the spatial resolution of global weather- and climate- prediction models makes it necessary to use numerically efficient and highly scalable algorithms to solve the equations of large scale atmospheric fluid dynamics. For stability and efficiency reasons several of the operational forecasting centres, in particular the Met Office and the ECMWF in the UK, use semi-implicit semi-Lagrangian time stepping in the dynamical core of the model. The additional burden with this approach is that a three dimensional elliptic partial differential equation (PDE) for the pressure correction has to be solved at every model time step and this often constitutes a significant proportion of the time spent in the dynamical core. To run within tight operational time scales the solver has to be parallelised and there seems to be a (perceived) misconception that elliptic solvers do not scale to large processor counts and hence implicit time stepping can not be used in very high resolution global models. After reviewing several methods for solving the elliptic PDE for the pressure correction and their application in atmospheric models we demonstrate the performance and very good scalability of Krylov subspace solvers and multigrid algorithms for a representative model equation with more than $10^{10}$ unknowns on 65536 cores on HECToR, the UK's national supercomputer. For this we tested and optimised solvers from two existing numerical libraries (DUNE and hypre) and implemented both a Conjugate Gradient solver and a geometric multigrid algorithm based on a tensor-product approach which exploits the strong vertical anisotropy of the discretised equation. We study both weak and strong scalability and compare the absolute solution times for all methods; in contrast to one-level methods the multigrid solver is robust with respect to parameter variations.