Massively parallel solvers for elliptic partial differential equations in numerical weather and climate prediction:scalability of elliptic solvers in NWP

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.Comment: 24 pages, 7 figures, 7 table

Matrix-free GPU implementation of a preconditioned conjugate gradient solver for anisotropic elliptic PDEs

Many problems in geophysical and atmospheric modelling require the fast solution of elliptic partial differential equations (PDEs) in "flat" three dimensional geometries. In particular, an anisotropic elliptic PDE for the pressure correction has to be solved at every time step in the dynamical core of many numerical weather prediction models, and equations of a very similar structure arise in global ocean models, subsurface flow simulations and gas and oil reservoir modelling. The elliptic solve is often the bottleneck of the forecast, and an algorithmically optimal method has to be used and implemented efficiently. Graphics Processing Units have been shown to be highly efficient for a wide range of applications in scientific computing, and recently iterative solvers have been parallelised on these architectures. We describe the GPU implementation and optimisation of a Preconditioned Conjugate Gradient (PCG) algorithm for the solution of a three dimensional anisotropic elliptic PDE for the pressure correction in NWP. Our implementation exploits the strong vertical anisotropy of the elliptic operator in the construction of a suitable preconditioner. As the algorithm is memory bound, performance can be improved significantly by reducing the amount of global memory access. We achieve this by using a matrix-free implementation which does not require explicit storage of the matrix and instead recalculates the local stencil. Global memory access can also be reduced by rewriting the algorithm using loop fusion and we show that this further reduces the runtime on the GPU. We demonstrate the performance of our matrix-free GPU code by comparing it to a sequential CPU implementation and to a matrix-explicit GPU code which uses existing libraries. The absolute performance of the algorithm for different problem sizes is quantified in terms of floating point throughput and global memory bandwidth.Comment: 18 pages, 7 figure

On the impact of three dimensional radiative transfer on cloud evolution

The goal of this study is to gain insight into cloud-radiative feedback mechanisms and what role three-dimensional radiative transfer effects play in the evolution of convective clouds. The usually employed one-dimensional radiative transfer solvers neglect any horizontal energy transfer and thereby introduce considerable errors in surface and atmospheric heating rates. While fully three-dimensional radiative transfer solvers exist, they are several orders of magnitude too slow. In conclusion, so far, there is no straightforward solution that would solve the task at hand — namely, compute accurate three-dimensional radiative heating rates in the atmosphere — fast enough to be coupled interactively to a cloud resolving model. This thesis presents a new method — the TenStream solver — that provides a fast yet accurate approximation for three-dimensional heating rates. The TenStream is furthermore integrated into the University of California, Los Angeles large-eddy simulation (UCLA-LES) cloud-resolving model. This setup allows to study the effects of three-dimensional radiative heating on the evolution of clouds. The TenStream method extends the well-known one-dimensional two-stream theory to 10 streams. The new solver significantly reduces the root mean square error for atmospheric heating and surface heating rates when compared to traditionally employed one-dimensional solvers. In the case of a cumulus cloud field and the solar zenith angle being 60 ◦ , the error is reduced from 178 % to 31 %. Parallel scalability was a primary concern developing the TenStream solver. This thesis documents the overall performance of the solver as well as the technical challenges of migrating from 1-D schemes to 3-D schemes. To understand the performance characteristics of the TenStream solver, weak as well as strong-scaling experiments are conducted. In this context, two matrix preconditioner are investigated: geometric algebraic multigrid preconditioning (GAMG) and block Jacobi incomplete LU (ILU) factorization and it is found that algebraic multigrid preconditioning performs well for complex scenes and highly parallelized simulations. The TenStream solver is tested on several state of the art super-computers for up to 4096 cores and shows a parallel scaling efficiency of 80 % to 90 %. The central part of this thesis examines the influence of three-dimensional radiative transfer effects on the development of convective cumulus clouds. The influence is tested on short time scales of a single convective warm-bubble and over a longer period of time and a reasonably large domain for shallow cumulus clouds. The directionality of the direct solar beam introduces an asymmetry in the atmospheric heating of the convective motion and tilts the updraft. While a cloud’s shadow is always directly beneath itself in a one-dimensional radiative transfer solver. In contrast, the TenStream solver correctly displaces the shadowy region according to the sun’s zenith angle. The constant supply of warm and moist air due to the local heating in the updraft region beneath the cloud, prolongs the cloud’s lifetime by a factor of two and generally increases cloud development. The influence of three-dimensional heating on the evolution of clouds shows to be persistent even in the presence of a horizontal wind. The results presented here motivate further research in the field of cloud- radiative feedbacks and their role in weather and climate prediction simulations

A generalized curvilinear solver for spherical shell Rayleigh-B\'enard convection

A three-dimensional finite-difference solver has been developed and implemented for Boussinesq convection in a spherical shell. The solver transforms any complex curvilinear domain into an equivalent Cartesian domain using Jacobi transformation and solves the governing equations in the latter. This feature enables the solver to account for the effects of the non-spherical shape of the convective regions of planets and stars. Apart from parallelization using MPI, implicit treatment of the viscous terms using a pipeline alternating direction implicit scheme and HYPRE multigrid accelerator for pressure correction makes the solver efficient for high-fidelity direct numerical simulations. We have performed simulations of Rayleigh-B\'enard convection at three Rayleigh numbers $Ra=10^{5}, 10^{7}$ and $10^{8}$ while keeping the Prandtl number fixed at unity ($Pr=1$). The average radial temperature profile and the Nusselt number match very well, both qualitatively and quantitatively, with the existing literature. Closure of the turbulent kinetic energy budget, apart from the relative magnitude of the grid spacing compared to the local Kolmogorov scales, assures sufficient spatial resolution

Solution of partial differential equations on vector and parallel computers

The present status of numerical methods for partial differential equations on vector and parallel computers was reviewed. The relevant aspects of these computers are discussed and a brief review of their development is included, with particular attention paid to those characteristics that influence algorithm selection. Both direct and iterative methods are given for elliptic equations as well as explicit and implicit methods for initial boundary value problems. The intent is to point out attractive methods as well as areas where this class of computer architecture cannot be fully utilized because of either hardware restrictions or the lack of adequate algorithms. Application areas utilizing these computers are briefly discussed

Research in Applied Mathematics, Fluid Mechanics and Computer Science

This report summarizes research conducted at the Institute for Computer Applications in Science and Engineering in applied mathematics, fluid mechanics, and computer science during the period October 1, 1998 through March 31, 1999

A Bayesian conjugate gradient method (with Discussion)

A fundamental task in numerical computation is the solution of large linear systems. The conjugate gradient method is an iterative method which offers rapid convergence to the solution, particularly when an effective preconditioner is employed. However, for more challenging systems a substantial error can be present even after many iterations have been performed. The estimates obtained in this case are of little value unless further information can be provided about the numerical error. In this paper we propose a novel statistical model for this numerical error set in a Bayesian framework. Our approach is a strict generalisation of the conjugate gradient method, which is recovered as the posterior mean for a particular choice of prior. The estimates obtained are analysed with Krylov subspace methods and a contraction result for the posterior is presented. The method is then analysed in a simulation study as well as being applied to a challenging problem in medical imaging