Precision Reduction for Elliptic Solvers

Semi-implicit timestepping schemes for atmosphere and ocean models require elliptic solvers that are working efficiently on modern supercomputers. We investigate potential computational savings when using mixed precision in the elliptic solvers to improve their performance. If precision is reduced from double precision arithmetic (64 bits per real number) to single precision arithmetic (32 bits per real number) in parts of the solver for operational weather forecasts, significant savings can be expected. While single precision is already used successfully within some models, the community remains sceptical whether the use of single precision for these solvers is a viable option since rounding errors may impair the results. To learn how numerical precision is affecting the performance of the elliptic solvers, we run important parts of a representative elliptic solver at even lower precisions (half precision with 16 bits per real number) and perform a detailed study on the impact on solver convergence and solution quality. We study this in the context of a shallow-water model on the sphere.

Our semi-implicit model is discretised based on the MPDATA approach. The resulting elliptic problem is solved using the generalized conjugated residual method, accelerated with an ADI-type preconditioner. As for the computational grid, here we deliberately choose a lat-lon grid, with its well-known pole problem, to better study the behavior of the reduced-precision in the vicinity of grid singularities. The precision reduction itself is done on a software level, using an emulator. In all experiments, the dynamic range of representable numbers is kept at double precision level (11 exponent bits). The reduced-precision experiments are performed for standard dynamical core test-cases like the zonal flow past a mountain and the propagation of the Rossby-Haurwitz wave (wavenumber 4) as well as for simulations of a zonal flow past realistic Earth orography.

We show for our setup that key components of the elliptic solve can be performed at the level of half precision (16 bits), most prominently the preconditioning. This in principle holds true for a large range of model resolutions. The only limitation we find when resolution is increased, is that for the reduced-precision the solver's convergence rate significantly deteriorates at the poles, as the pole problem becomes more dominant. We show however that the solver performance can be almost entirely recovered by either employing higher levels of precision in a small vicinity of the poles, or alternatively by applying a suitable smoothing of the prognostic model variables near the poles.