Colin J. Cotter

Andrew T. T. McRae, Colin J. Cotter, Chris J. Budd

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-Ampère equation for a scalar mesh potential. We adapt an existing finite element scheme for the standard Monge-Ampère 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-Ampère equation due to the implicit dependence of the monitor function on the resulting mesh. We also derive the equivalent Monge-Ampère-like equation 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.

C. J. Cotter, J. Thuburn

We describe discretisations of the shallow water equations on the sphere using the framework of finite element exterior calculus, which are extensions of the mimetic finite difference framework presented in Ringler, Thuburn, Klemp, and Skamarock (Journal of Computational Physics, 2010). The exterior calculus notation provides a guide to which finite element spaces should be used for which physical variables, and unifies a number of desirable properties. We present two formulations: a ``primal'' formulation in which the finite element spaces are defined on a single mesh, and a ``primal-dual'' formulation in which finite element spaces on a dual mesh are also used. Both formulations have velocity and layer depth as prognostic variables, but the exterior calculus framework leads to a conserved diagnostic potential vorticity. In both formulations we show how to construct discretisations that have mass-consistent (constant potential vorticity stays constant), stable and oscillation-free potential vorticity advection.

Alastair Gregory, Colin Cotter, Sebastian Reich

This paper extends the Multilevel Monte Carlo variance reduction technique to nonlinear filtering. In particular, Multilevel Monte Carlo is applied to a certain variant of the particle filter, the Ensemble Transform Particle Filter. A key aspect is the use of optimal transport methods to re-establish correlation between coarse and fine ensembles after resampling; this controls the variance of the estimator. Numerical examples present a proof of concept of the effectiveness of the proposed method, demonstrating significant computational cost reductions (relative to the single-level ETPF counterpart) in the propagation of ensembles.

C. J. Cotter, D. A. Ham, C. C. Pain et al.

We introduce a new mixed discontinuous/continuous Galerkin finite element for solving the 2- and 3-dimensional wave equations and equations of incompressible flow. The element, which we refer to as P1dg-P2, uses discontinuous piecewise linear functions for velocity and continuous piecewise quadratic functions for pressure. The aim of introducing the mixed formulation is to produce a new flexible element choice for triangular and tetrahedral meshes which satisfies the LBB stability condition and hence has no spurious zero-energy modes. We illustrate this property with numerical integrations of the wave equation in two dimensions, an analysis of the resultant discrete Laplace operator in two and three dimensions, and a normal mode analysis of the semi-discrete wave equation in one dimension.

Colin J. Cotter, Sebastian Reich

This survey paper is written with the intention of giving a mathematical introduction to filtering techniques for intermittent data assimilation, and to survey some recent advances in the field. The paper is divided into three parts. The first part introduces Bayesian statistics and its application to statistical inference and estimation. Basic aspects of Markov processes, as they typically arise from scientific models in the form of stochastic differential and/or difference equations, are covered in the second part. The third and final part describes the filtering approach to estimation of model states by assimilation of observational data into scientific models. While most of the material is of survey type, very recent advances in the field of nonlinear data assimilation covered in this paper include a discussion of Bayesian inference in the context of optimal transportation and coupling of random variables, as well as a discussion of recent advances in ensemble transform filters. References and sources for further reading material will be listed at the end of each section.

Colin Cotter, Simon Cotter, Paul Russell

Markov chain Monte Carlo methods are a powerful and commonly used family of numerical methods for sampling from complex probability distributions. As applications of these methods increase in size and complexity, the need for efficient methods increases. In this paper, we present a particle ensemble algorithm. At each iteration, an importance sampling proposal distribution is formed using an ensemble of particles. A stratified sample is taken from this distribution and weighted under the posterior, a state-of-the-art ensemble transport resampling method is then used to create an evenly weighted sample ready for the next iteration. We demonstrate that this ensemble transport adaptive importance sampling (ETAIS) method outperforms MCMC methods with equivalent proposal distributions for low dimensional problems, and in fact shows better than linear improvements in convergence rates with respect to the number of ensemble members. We also introduce a new resampling strategy, multinomial transformation (MT), which while not as accurate as the ensemble transport resampler, is substantially less costly for large ensemble sizes, and can then be used in conjunction with ETAIS for complex problems. We also focus on how algorithmic parameters regarding the mixture proposal can be quickly tuned to optimise performance. In particular, we demonstrate this methodology's superior sampling for multimodal problems, such as those arising from inference for mixture models, and for problems with expensive likelihoods requiring the solution of a differential equation, for which speed-ups of orders of magnitude are demonstrated. Likelihood evaluations of the ensemble could be computed in a distributed manner, suggesting that this methodology is a good candidate for parallel Bayesian computations.

Colin J. Cotter, Robert C. Kirby

We study mixed finite element methods for the linearized rotating shallow water equations with linear drag and forcing terms. By means of a strong energy estimate for an equivalent second-order formulation for the linearized momentum, we prove long-time stability of the system without energy accumulation -- the geotryptic state. A priori error estimates for the linearized momentum and free surface elevation are given in $L^2$ as well as for the time derivative and divergence of the linearized momentum. Numerical results confirm the theoretical results regarding both energy damping and convergence rates.

Andreas Bock, Colin J. Cotter, Robert C. Kirby

We study parameterisation-independent closed planar curve matching as a Bayesian inverse problem. The motion of the curve is modelled via a curve on the diffeomorphism group acting on the ambient space, leading to a large deformation diffeomorphic metric mapping (LDDMM) functional penalising the kinetic energy of the deformation. We solve Hamilton's equations for the curve matching problem using the Wu-Xu element [S. Wu, J. Xu, Nonconforming finite element spaces for $2m^\text{th}$ order partial differential equations on $\mathbb{R}^n$ simplicial grids when $m=n+1$, Mathematics of Computation 88 (316) (2019) 531-551] which provides mesh-independent Lipschitz constants for the forward motion of the curve, and solve the inverse problem for the momentum using Bayesian inversion. Since this element is not affine-equivalent we provide a pullback theory which expedites the implementation and efficiency of the forward map. We adopt ensemble Kalman inversion using a negative Sobolev norm mismatch penalty to measure the discrepancy between the target and the ensemble mean shape. We provide several numerical examples to validate the approach.

Thomas M. Bendall, Colin J. Cotter

A framework of variational principles for stochastic fluid dynamics was presented by Holm (2015), and these stochastic equations were also derived by Cotter et al. (2017). We present a conforming finite element discretisation for the stochastic quasi-geostrophic equation that was derived from this framework. The discretisation preserves the first two moments of potential vorticity, i.e. the mean potential vorticity and the enstrophy. Following the work of Dubinkina and Frank (2007), who investigated the statistical mechanics of discretisations of the deterministic quasi-geostrophic equation, we investigate the statistical mechanics of our discretisation of the stochastic quasi-geostrophic equation. We compare the statistical properties of our discretisation with the Gibbs distribution under assumption of these conserved quantities, finding that there is agreement between the statistics under a wide range of set-ups.

Andreas Bock, Colin J. Cotter

We study the problem of diffeomorphometric geodesic landmark matching where the objective is to find a diffeomorphism that via its group action maps between two sets of landmarks. It is well-known that the motion of the landmarks, and thereby the diffeomorphism, can be encoded by an initial momentum leading to a formulation where the landmark matching problem can be solved as an optimisation problem over such momenta. The novelty of our work lies in the application of a derivative-free Bayesian inverse method for learning the optimal momentum encoding the diffeomorphic mapping between the template and the target. The method we apply is the ensemble Kalman filter, an extension of the Kalman filter to nonlinear observation operators. We describe an efficient implementation of the algorithm and show several numerical results for various target shapes.

Andreas Bock, Alexis Arnaudon, Colin Cotter

We present a framework for shape matching in computational anatomy allowing users control of the degree to which the matching is diffeomorphic. This control is given as a function defined over the image and parameterises the template deformation. By modelling localised template deformation we have a mathematical description of growth only in specified parts of an image. The location can either be specified from prior knowledge of the growth location or learned from data. For simplicity, we consider landmark matching and infer the distribution of a finite dimensional parameterisation of the control via Markov chain Monte Carlo. Preliminary numerical results are shown and future paths of investigation are laid out. Well-posedness of this new problem is studied together with an analysis of the associated geodesic equations.

Alastair Gregory, Colin Cotter

This paper presents a seamless algorithm for the application of the multilevel Monte Carlo (MLMC) method to the ensemble transform particle filter (ETPF). The algorithm uses a combination of optimal coupling transformations between coarse and fine ensembles in difference estimators within a multilevel framework, to minimise estimator variance. It differs from that of Gregory et al. (2016) in that strong coupling between the coarse and fine ensembles is seamlessly maintained during all stages of the assimilation algorithm, instead of using independent transformations to equal weights followed by recoupling with an assignment problem. This modification is found to lead to an increased rate in variance decay between coarse and fine ensembles with level in the hierarchy, a key component of MLMC. This offers the potential for greater computational cost reductions. This is shown, alongside evidence of asymptotic consistency, in numerical examples.

Colin J. Cotter, P. Jameson Graber, Robert C. Kirby

We study mixed finite element methods for the rotating shallow water equations with linearized momentum terms but nonlinear drag. By means of an equivalent second-order formulation, we prove long-time stability of the system without energy accumulation. We also give rates of damping in unforced systems and various continuous dependence results on initial conditions and forcing terms. \emph{A priori} error estimates for the momentum and free surface elevation are given in $L^2$ as well as for the time derivative and divergence of the momentum. Numerical results confirm the theoretical results regarding both energy damping and convergence rates.

C. J. Cotter, D. A. Ham, A. T. T. McRae et al.

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.

C. J. Cotter

Inertia-gravity mode and Rossby mode dispersion properties are examined for discretisations of the linearized rotating shallow-water equations using the $P1_{DG}$-$P2$ finite element pair on arbitrary triangulations in planar geometry. A discrete Helmholtz decomposition of the functions in the velocity space based on potentials taken from the pressure space is used to provide a complete description of the numerical wave propagation for the discretised equations. In the $f$-plane case, this decomposition is used to obtain decoupled equations for the geostrophic modes, the inertia-gravity modes, and the inertial oscillations. As has been noticed previously, the geostrophic modes are steady. The Helmholtz decomposition is used to show that the resulting inertia-gravity wave equation is third-order accurate in space. In general the \pdgp finite element pair is second-order accurate, so this leads to very accurate wave propagation. It is further shown that the only spurious modes supported by this discretisation are spurious inertial oscillations which have frequency $f$, and which do not propagate. The Helmholtz decomposition also allows a simple derivation of the quasi-geostrophic limit of the discretised $P1_{DG}$-$P2$ equations in the $β$-plane case, resulting in a Rossby wave equation which is also third-order accurate.

C. J. Cotter

We introduce a family of mixed finite element pairs for use on geodesic grids and with adaptive mesh refinement for numerical weather prediction and ocean modelling. We prove that when these finite element pairs are applied to the linear rotating shallow water equations, the geostrophically balanced states are exactly steady, which means that the numerical schemes do not introduce any spurious inertia-gravity waves; this makes these finite element pairs in some sense optimal for numerical weather prediction and ocean modelling applications. We further prove that these finite element pairs satisfy an inf-sup condition which means that they are free of spurious pressure modes which would pollute the numerical solution over the timescales required for large-scale geophysical applications. We then discuss the extension to incompressible Euler-Boussinesq equations with rotation, and show that for the linearised equations the balanced states are again exactly steady on arbitrary unstructured meshes. We also show that the discrete pressure Poisson equation resulting from these discretisations satisfies an optimal stencil property. All these properties make the discretisations in this family excellent candidates for numerical weather prediction and large-scale ocean modelling applications when unstructured grids are required.

Colin J Cotter, Darryl D Holm

We focus on the spatial discretization produced by the Variational Particle-Mesh (VPM) method for a prototype fluid equation the known as the EPDiff equation}, which is short for Euler-Poincaré equation associated with the diffeomorphism group (of $\mathbb{R}^d$, or of a $d$-dimensional manifold $Ω$). The EPDiff equation admits measure valued solutions, whose dynamics are determined by the momentum maps for the left and right actions of the diffeomorphisms on embedded subspaces of $\mathbb{R}^d$. The discrete VPM analogs of those dynamics are studied here. Our main results are: (i) a variational formulation for the VPM method, expressed in terms of a constrained variational principle principle for the Lagrangian particles, whose velocities are restricted to a distribution $D_{\VPM}$ which is a finite-dimensional subspace of the Lie algebra of vector fields on $Ω$; (ii) a corresponding constrained variational principle on the fixed Eulerian grid which gives a discrete version of the Euler-Poincaré equation; and (iii) discrete versions of the momentum maps for the left and right actions of diffeomorphisms on the space of solutions.

Colin Cotter

Given a fluid equation with reduced Lagrangian $l$ which is a functional of velocity $\MM{u}$ and advected density $D$ given in Eulerian coordinates, we give a general method for semidiscretising the equations to give a canonical Hamiltonian system; this system may then be integrated in time using a symplectic integrator. The method is Lagrangian, with the variables being a set of Lagrangian particle positions and their associated momenta. The canonical equations obtained yield a discrete form of Euler-Poincaré equations for $l$ when projected onto the grid, with a new form of discrete calculus to represent the gradient and divergence operators. Practical symplectic time integrators are suggested for a large family of equations which include the shallow-water equations, the EP-Diff equations and the 3D compressible Euler equations, and we illustrate the technique by showing results from a numerical experiment for the EP-Diff equations.