# Publications by Andrew McRae

## Compatible Finite Element Methods for Geophysical Flows Automation and Implementation Using Firedrake

Springer, 2019

This book introduces recently developed mixed finite element methods for large-scale geophysical flows that preserve essential numerical properties for accurate simulations.

## The scaling and skewness of optimally transported meshes on the sphere

Journal of Computational Physics Elsevier **375** (2018) 540-564

In the context of numerical solution of PDEs, dynamic mesh redistribution methods (r-adaptive methods) are an important procedure for increasing the resolution in regions of interest, without modifying the connectivity of the mesh. Key to the success of these methods is that the mesh should be sufficiently refined (locally) and flexible in order to resolve evolving solution features, but at the same time not introduce errors through skewness and lack of regularity. Some state-of-the-art methods are bottom-up in that they attempt to prescribe both the local cell size and the alignment to features of the solution. However, the resulting problem is overdetermined, necessitating a compromise between these conflicting requirements. An alternative approach, described in this paper, is to prescribe only the local cell size and augment this an optimal transport condition to provide global regularity. This leads to a robust and flexible algorithm for generating meshes fitted to an evolving solution, with minimal need for tuning parameters. Of particular interest for geophysical modelling are meshes constructed on the surface of the sphere. The purpose of this paper is to demonstrate that meshes generated on the sphere using this optimal transport approach have good a-priori regularity and that the meshes produced are naturally aligned to various simple features. It is further shown that the sphere's intrinsic curvature leads to more regular meshes than the plane. In addition to these general results, we provide a wide range of examples relevant to practical applications, to showcase the behaviour of optimally transported meshes on the sphere. These range from axisymmetric cases that can be solved analytically to more general examples that are tackled numerically. Evaluation of the singular values and singular vectors of the mesh transformation provides a quantitative measure of the mesh aniso...

## Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements

SIAM Journal on Scientific Computing Society for Industrial and Applied Mathematics **40** (2018) A1121-A1148

In moving mesh methods, the underlying mesh is dynamically adapted without changing the connectivity of the mesh. We specifically consider the generation of meshes which are adapted to a scalar monitor function through equidistribution. Together with an optimal transport condition, this leads to a Monge–Ampere equation ` for a scalar mesh potential. We adapt an existing finite element scheme for the standard Monge–Ampere ` equation to this mesh generation problem; this is a mixed finite element scheme, in which an extra discrete variable is introduced to represent the Hessian matrix of second derivatives. The problem we consider has additional nonlinearities over the basic Monge–Ampere equation due to the implicit dependence of the monitor func- ` tion on the resulting mesh. We also derive an equivalent Monge–Ampere-like equa- ` tion for generating meshes on the sphere. The finite element scheme is extended to the sphere, and we provide numerical examples. All numerical experiments are performed using the open-source finite element framework Firedrake.

## Firedrake: Automating the finite element method by composing abstractions

ACM Transactions on Mathematical Software **43** (2016)

Firedrake is a new tool for automating the numerical solution of partial differential equations. Firedrake adopts the domain-specific language for the finite element method of the FEniCS project, but with a pure Python runtime-only implementation centered on the composition of several existing and new abstractions for particular aspects of scientific computing. The result is a more complete separation of concerns that eases the incorporation of separate contributions from computer scientists, numerical analysts, and application specialists. These contributions may add functionality or improve performance. Firedrake benefits from automatically applying new optimizations. This includes factorizing mixed function spaces, transforming and vectorizing inner loops, and intrinsically supporting block matrix operations. Importantly, Firedrake presents a simple public API for escaping the UFL abstraction. This allows users to implement common operations that fall outside of pure variational formulations, such as flux limiters.

## A structure-exploiting numbering algorithm for finite elements on extruded meshes, and its performance evaluation in Firedrake

Geoscientific Model Development **9** (2016) 3803-3815

We present a generic algorithm for numbering and then efficiently iterating over the data values attached to an extruded mesh. An extruded mesh is formed by replicating an existing mesh, assumed to be unstructured, to form layers of prismatic cells. Applications of extruded meshes include, but are not limited to, the representation of three-dimensional high aspect ratio domains employed by geophysical finite element simulations. These meshes are structured in the extruded direction. The algorithm presented here exploits this structure to avoid the performance penalty traditionally associated with unstructured meshes. We evaluate the implementation of this algorithm in the Firedrake finite element system on a range of low compute intensity operations which constitute worst cases for data layout performance exploration. The experiments show that having structure along the extruded direction enables the cost of the indirect data accesses to be amortized after 10-20 layers as long as the underlying mesh is well ordered. We characterize the resulting spatial and temporal reuse in a representative set of both continuous-Galerkin and discontinuous-Galerkin discretizations. On meshes with realistic numbers of layers the performance achieved is between 70 and 90% of a theoretical hardware-specific limit.

## AUTOMATED GENERATION AND SYMBOLIC MANIPULATION OF TENSOR PRODUCT FINITE ELEMENTS

SIAM JOURNAL ON SCIENTIFIC COMPUTING **38** (2016) S25-S47

## Energy- and enstrophy-conserving schemes for the shallow-water equations, based on mimetic finite elements

QUARTERLY JOURNAL OF THE ROYAL METEOROLOGICAL SOCIETY **140** (2014) 2223-2234

## Automating the solution of PDEs on the sphere and other manifolds in FEniCS 1.2

GEOSCIENTIFIC MODEL DEVELOPMENT **6** (2013) 2099-2119

## Using reduced-precision arithmetic in the adjoint model of MITgcm

ArXiv (0)

In recent years, it has been convincingly shown that weather forecasting models can be run in single-precision arithmetic. Several models or components thereof have been tested with even lower precision than this. This previous work has largely focused on the main nonlinear `forward' model. A nonlinear model (in weather forecasting or otherwise) can have corresponding tangent linear and adjoint models, which are used in 4D variational data assimilation. The linearised models are plausibly far more sensitive to reductions in numerical precision since unbounded error growth can occur with no possibility of nonlinear saturation. In this paper, we present a geophysical experiment that makes use of an adjoint model to calculate sensitivities and perform optimisation. Using software emulation, we investigate the effect of degrading the numerical precision of the adjoint model. We find that reasonable results are obtained with as few as 10 significand bits, equal to the significand precision in the IEEE half-precision standard.

## On the shallow atmosphere approximation in finite element dynamical cores

ArXiv (0)

We provide an approach to implementing the shallow atmosphere approximation in three dimensional finite element discretisations for dynamical cores. The approach makes use of the fact that the shallow atmosphere approximation metric can be obtained by writing equations on a three-dimensional manifold embedded in $\mathbb{R}^4$ with a restriction of the Euclidean metric. We show that finite element discretisations constructed this way are equivalent to the use of a modified three dimensional mesh for the construction of metric terms. We demonstrate our approach via a convergence test for a prototypical elliptic problem.

## Compatible finite element methods for numerical weather prediction

ArXiv (0)

This article takes the form of a tutorial on the use of a particular class of mixed finite element methods, which can be thought of as the finite element extension of the C-grid staggered finite difference method. The class is often referred to as compatible finite elements, mimetic finite elements, discrete differential forms or finite element exterior calculus. We provide an elementary introduction in the case of the one-dimensional wave equation, before summarising recent results in applications to the rotating shallow water equations on the sphere, before taking an outlook towards applications in three-dimensional compressible dynamical cores.