Search

Accurate transport algorithms are crucial for computational fluid dynamics and more accurate and efficient schemes are always in development. One dimensional limiting is commonly employed to suppress nonphysical oscillations. However, the application of such limiters can reduce accuracy. It is important to identify the weakest set of sufficient conditions required on the limiter as to allow the development of successful numerical algorithms. The main goal of this paper is to identify new less restrictive sufficient conditions for flux form in-compressible advection to remain monotonic. We identify additional necessary conditions for incompressible flux form advection to be monotonic, demonstrating that the Spekreijse limiter region is not sufficient for incompressible flux form advection to remain monotonic. Then a convex combination argument is used to derive new sufficient conditions that are less restrictive than the Sweby region for a discrete maximum principle. This allows the introduction of two new more general limiter regions suitable for flux form incompressible advection.

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., Comment: 45 pages, 9 figures

When estimating quantities and fields that are difficult to measure directly, such as the fluidity of ice, from point data sources, such as satellite altimetry, it is important to solve a numerical inverse problem that is formulated with Bayesian consistency. Otherwise, the resultant probability density function for the difficult to measure quantity or field will not be appropriately clustered around the truth. In particular, the inverse problem should be formulated by evaluating the numerical solution at the true point locations for direct comparison with the point data source. If the data are first fitted to a gridded or meshed field on the computational grid or mesh, and the inverse problem formulated by comparing the numerical solution to the fitted field, the benefits of additional point data values below the grid density will be lost. We demonstrate, with examples in the fields of groundwater hydrology and glaciology, that a consistent formulation can increase the accuracy of results and aid discourse between modellers and observationalists. To do this, we bring point data into the finite element method ecosystem as discontinuous fields on meshes of disconnected vertices. Point evaluation can then be formulated as a finite element interpolation operation (dual-evaluation). This new abstraction is well-suited to automation, including automatic differentiation. We demonstrate this through implementation in Firedrake, which generates highly optimised code for solving Partial Differential Equations (PDEs) with the finite element method. Our solution integrates with dolfin-adjoint/pyadjoint, allowing PDE-constrained optimisation problems, such as data assimilation, to be solved through forward and adjoint mode automatic differentiation., Comment: This version: Added missing affiliation

This article surveys research on the application of compatible finite element methods to large scale atmosphere and ocean simulation. Compatible finite element methods extend Arakawa's C-grid finite difference scheme to the finite element world. They are constructed from a discrete de Rham complex, which is a sequence of finite element spaces which are linked by the operators of differential calculus. The use of discrete de Rham complexes to solve partial differential equations is well established, but in this article we focus on the specifics of dynamical cores for simulating weather, oceans and climate. The most important consequence of the discrete de Rham complex is the Hodge-Helmholtz decomposition, which has been used to exclude the possibility of several types of spurious oscillations from linear equations of geophysical flow. This means that compatible finite element spaces provide a useful framework for building dynamical cores. In this article we introduce the main concepts of compatible finite element spaces, and discuss their wave propagation properties. We survey some methods for discretising the transport terms that arise in dynamical core equation systems, and provide some example discretisations, briefly discussing their iterative solution. Then we focus on the recent use of compatible finite element spaces in designing structure preserving methods, surveying variational discretisations, Poisson bracket discretisations, and consistent vorticity transport., Comment: correction of some typos

Compatible finite element discretisations for the atmospheric equations of motion have recently attracted considerable interest. Semi-implicit timestepping methods require the repeated solution of a large saddle-point system of linear equations. Preconditioning this system is challenging since the velocity mass matrix is non-diagonal, leading to a dense Schur complement. Hybridisable discretisations overcome this issue: weakly enforcing continuity of the velocity field with Lagrange multipliers leads to a sparse system of equations, which has a similar structure to the pressure Schur complement in traditional approaches. We describe how the hybridised sparse system can be preconditioned with a non-nested two-level preconditioner. To solve the coarse system, we use the multigrid pressure solver that is employed in the approximate Schur complement method previously proposed by the some of the authors. Our approach significantly reduces the number of solver iterations. The method shows excellent performance and scales to large numbers of cores in the Met Office next-generation climate- and weather prediction model LFRic., Comment: 24 pages, 13 figures, 5 tables; accepted for publication in Quarterly Journal of the Royal Meteorological Society

We derive a linearized rotating shallow water system modeling tides, which can be discretized by mixed finite elements. Unlike previous models, this model allows for multiple layers stratified by density. Like the single-layer case~\cite{kirby2021preconditioning} a weighted-norm preconditioner gives a (nearly) parameter-robust method for solving the resulting linear system at each time step, but the all-to-all coupling between the layers in the model poses a significant challenge to efficiency. Neglecting the inter-layer coupling gives a preconditioner that degrades rapidly as the number of layers increases. By a careful analysis of the matrix that couples the layers, we derive a robust method that requires solving a reformulated system that only involves coupling between adjacent layers. Numerical results obtained using Firedrake confirm the theory., Comment: 23 pages, 7 figures

We present a new implementation of the geometric method of Cullen & Purser (1984) for solving the semi-geostrophic Eady slice equations which model large scale atmospheric flows and frontogenesis. The geometric method is a Lagrangian discretisation, where the PDE is approximated by a particle system. An important property of the discretisation is that it is energy conserving. We restate the geometric method in the language of semi-discrete optimal transport theory and exploit this to develop a fast implementation that combines the latest results from numerical optimal transport theory with a novel adaptive time-stepping scheme. Our results enable a controlled comparison between the Eady-Boussinesq vertical slice equations and their semi-geostrophic approximation. We provide further evidence that weak solutions of the Eady-Boussinesq vertical slice equations converge to weak solutions of the semi-geostrophic Eady slice equations as the Rossby number tends to zero., Comment: 43 pages, 9 figures

Adaptive meshes have the potential to improve the accuracy and efficiency of atmospheric modelling by increasing resolution where it is most needed. Mesh re-distribution, or r-adaptivity, adapts by moving the mesh without changing the connectivity. This avoids some of the challenges with h-adaptivity (adding and removing points): the solution does not need to be mapped between meshes, which can be expensive and introduces errors, and there are no load balancing problems on parallel computers. A long standing problem with both forms of adaptivity has been changes in volume of the domain as resolution changes at an uneven boundary. We propose a solution to exact local conservation and maintenance of uniform fields while the mesh changes volume as it moves over orography. This is solved by introducing a volume adjustment parameter which tracks the true cell volumes without using expensive conservative mapping. A finite volume solution of the advection equation over orography on moving meshes is described and results are presented demonstrating improved accuracy for cost using moving meshes. Exact local conservation and maintenance of uniform fields is demonstrated and the corrected mesh volume is preserved. We use optimal transport to generate meshes which are guaranteed not to tangle and are equidistributed with respect to a monitor function. This leads to a Monge-Amp\`{e}re equation which is solved with a Newton solver. The superiority of the Newton solver over other techniques is demonstrated in the appendix. However the Newton solver is only efficient if it is applied to the left hand side of the Monge-Amp\`{e}re equation with fixed point iterations for the right hand side.

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., Comment: 19 pages, 8 figures

We describe a proof-of-concept development and application of a phase averaging technique to the nonlinear rotating shallow water equations on the sphere, discretised using compatible finite element methods. Phase averaging consists of averaging the nonlinearity over phase shifts in the exponential of the linear wave operator. Phase averaging aims to capture the slow dynamics in a solution that is smoother in time (in transformed variables) so that larger timesteps may be taken. We overcome the two key technical challenges that stand in the way of studying the phase averaging and advancing its implementation: 1) we have developed a stable matrix exponential specific to finite elements and 2) we have developed a parallel finite averaging proceedure. Following Peddle et al (2019), we consider finite width phase averaging windows, since the equations have a finite timescale separation. In our numerical implementation, the averaging integral is replaced by a Riemann sum, where each term can be evaluated in parallel. This creates an opportunity for parallelism in the timestepping method, which we use here to compute our solutions. Here, we focus on the stability and accuracy of the numerical solution. We confirm there is an optimal averaging window, in agreement with theory. Critically, we observe that the combined time discretisation and averaging error is much smaller than the time discretisation error in a semi-implicit method applied to the same spatial discretisation. An evaluation of the parallel aspects will follow in later work., Comment: Author Beth Wingate added for submitted version

We introduce a higher order phase averaging method for nonlinear oscillatory systems. Phase averaging is a technique to filter fast motions from the dynamics whilst still accounting for their effect on the slow dynamics. Phase averaging is useful for deriving reduced models that can be solved numerically with more efficiency, since larger timesteps can be taken. Recently, Haut and Wingate (2014) introduced the idea of computing finite window numerical phase averages in parallel as the basis for a coarse propagator for a parallel-in-time algorithm. In this contribution, we provide a framework for higher order phase averages that aims to better approximate the unaveraged system whilst still filtering fast motions. Whilst the basic phase average assumes that the solution independent of changes of phase, the higher order method expands the phase dependency in a basis which the equations are projected onto. In this new framework, the original numerical phase averaging formulation arises as the lowest order version of this expansion. Our new projection onto functions that are $k$th degree polynomials in the phase gives rise to higher order corrections to the phase averaging formulation. We illustrate the properties of this method on an ODE that describes the dynamics of a swinging spring due to Lynch (2002). Although idealized, this model shows an interesting analogy to geophysical flows as it exhibits a slow dynamics that arises through the resonance between fast oscillations. On this example, we show convergence to the non-averaged (exact) solution with increasing approximation order also for finite averaging windows. At zeroth order, our method coincides with a standard phase average, but at higher order it is more accurate in the sense that solutions of the phase averaged model track the solutions of the unaveraged equations more accurately.

We present an energy conserving space discretisation based on a Poisson bracket that can be used to derive the dry compressible Euler as well as thermal shallow water equations. It is formulated using the compatible finite element method, and extends the incorporation of upwinding for the shallow water equations as described in Wimmer, Cotter, and Bauer (2019). While the former is restricted to DG upwinding, an energy conserving SUPG scheme for the (partially) continuous Galerkin thermal field space is newly introduced here. The energy conserving property is validated by coupling the Poisson bracket based spatial discretisation to an energy conserving time discretisation. Further, the discretisation is demonstrated to lead to an improved temperature field development with respect to stability when upwinding is included. An approximately energy conserving full discretisation with a smaller computational cost is also presented., Comment: 32 pages, 12 figures; main differences to first version: simplified formulation for auxiliary SUPG operator, improved numerical results section

We introduce an efficient split finite element (FE) discretization of a y-independent (slice) model of the rotating shallow water equations. The study of this slice model provides insight towards developing schemes for the full 2D case. Using the split Hamiltonian FE framework (Bauer, Behrens and Cotter, 2019), we result in structure-preserving discretizations that are split into topological prognostic and metric-dependent closure equations. This splitting also accounts for the schemes' properties: the Poisson bracket is responsible for conserving energy (Hamiltonian) as well as mass, potential vorticity and enstrophy (Casimirs), independently from the realizations of the metric closure equations. The latter, in turn, determine accuracy, stability, convergence and discrete dispersion properties. We exploit this splitting to introduce structure-preserving approximations of the mass matrices in the metric equations avoiding to solve linear systems. We obtain a fully structure-preserving scheme with increased efficiency by a factor of two., Comment: 8 pages, conference paper

A famous feature of the Camassa-Holm equation is its admission of peaked soliton solutions known as peakons. We investigate this equation under the influence of stochastic transport. Noting that peakons are weak solutions of the equation, we present a finite element discretisation for it, which we use to explore the formation of peakons. Our simulations using this discretisation reveal that peakons can still form in the presence of stochastic perturbations. Peakons can emerge both through wave breaking, as the slope turns vertical, and without wave breaking as the inflection points of the velocity profile rise to reach the summit.

We present new discretisation of the moist compressible Euler equations, using the compatible finite element framework identified in Cotter and Shipton (2012). The discretisation strategy is described and details of the parametrisations of moist processes are presented. A procedure for establishing hydrostatic balance for moist atmospheres is introduced, and the model's performance is demonstrated through several test cases, two of which are new.

We present a new compatible finite element advection scheme for the compressible Euler equations. Unlike the discretisations described in Cotter and Kuzmin (2016) and Shipton et al (2018), the discretisation uses the lowest-order family of compatible finite element spaces, but still retains second-order numerical accuracy. This scheme obtains this second-order accuracy by first `recovering' the function in higher-order spaces, before using the discontinuous Galerkin advection schemes of Cotter and Kuzmin (2016). As well as describing the scheme, we also present its stability properties and a strategy for ensuring boundedness. We then demonstrate its properties through some numerical tests, before presenting its use within a model solving the compressible Euler equations.

We consider the variational formulation for vertical slice models introduced in Cotter and Holm (Proc Roy Soc, 2013). These models have a Kelvin circulation theorem that holds on all materially-transported closed loops, not just those loops on isosurfaces of potential temperature. Potential vorticity conservation can be derived directly from this circulation theorem. In this paper, we show that this property is due to these models having a relabelling symmetry for every single diffeomorphism of the vertical slice that preserves the density, not just those diffeomorphisms that preserve the potential temperature. This is developed using the methodology of Cotter and Holm (Foundations of Computational Mathematics, 2012).

Within the finite element community, discontinuous Galerkin (DG) and mixed finite element methods have become increasingly popular in simulating geophysical flows. However, robust and efficient solvers for the resulting saddle-point and elliptic systems arising from these discretizations continue to be an on-going challenge. One possible approach for addressing this issue is to employ a method known as hybridization, where the discrete equations are transformed such that classic static condensation and local post-processing methods can be employed. However, it is challenging to implement hybridization as performant parallel code within complex models, whilst maintaining separation of concerns between applications scientists and software experts. In this paper, we introduce a domain-specific abstraction within the Firedrake finite element library that permits the rapid execution of these hybridization techniques within a code-generating framework. The resulting framework composes naturally with Firedrake's solver environment, allowing for the implementation of hybridization and static condensation as runtime-configurable preconditioners via the Python interface to PETSc, petsc4py. We provide examples derived from second order elliptic problems and geophysical fluid dynamics. In addition, we demonstrate that hybridization shows great promise for improving the performance of solvers for mixed finite element discretizations of equations related to large-scale geophysical flows., Comment: Revisions for submission to GMD

We present a numerical investigation of stochastic transport in ideal fluids. According to Holm (Proc Roy Soc, 2015) and Cotter et al. (2017), the principles of transformation theory and multi-time homogenisation, respectively, imply a physically meaningful, data-driven approach for decomposing the fluid transport velocity into its drift and stochastic parts, for a certain class of fluid flows. In the current paper, we develop new methodology to implement this velocity decomposition and then numerically integrate the resulting stochastic partial differential equation using a finite element discretisation for incompressible 2D Euler fluid flows. The new methodology tested here is found to be suitable for coarse graining in this case. Specifically, we perform uncertainty quantification tests of the velocity decomposition of Cotter et al. (2017), by comparing ensembles of coarse-grid realisations of solutions of the resulting stochastic partial differential equation with the "true solutions" of the deterministic fluid partial differential equation, computed on a refined grid. The time discretization used for approximating the solution of the stochastic partial differential equation is shown to be consistent. We include comprehensive numerical tests that confirm the non-Gaussianity of the stream function, velocity and vorticity fields in the case of incompressible 2D Euler fluid flows., Comment: 41 pages, 26 figures Minor changes -- updated figures to improve readability. Corrected typos. Shifted Remark 7 to just after Assumption A1. Added Remark 8

We describe an energy-enstrophy conserving discretisation for the rotating shallow water equations with slip boundary conditions. This relaxes the assumption of boundary-free domains (periodic solutions or the surface of a sphere, for example) in the energy-enstrophy conserving formulation of McRae and Cotter (2014). This discretisation requires extra prognostic vorticity variables on the boundary in addition to the prognostic velocity and layer depth variables. The energy-enstrophy conservation properties hold for any appropriate set of compatible finite element spaces defined on arbitrary meshes with arbitrary boundaries. We demonstrate the conservation properties of the scheme with numerical solutions on a rotating hemisphere., Comment: Revision in response to reviewers at JCP

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..., Comment: Updated following reviewer comments

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., Comment: updated in peer review

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.

In {\em{Holm}, Proc. Roy. Soc. A 471 (2015)} stochastic fluid equations were derived by employing a variational principle with an assumed stochastic Lagrangian particle dynamics. Here we show that the same stochastic Lagrangian dynamics naturally arises in a multi-scale decomposition of the deterministic Lagrangian flow map into a slow large-scale mean and a rapidly fluctuating small scale map. We employ homogenization theory to derive effective slow stochastic particle dynamics for the resolved mean part, thereby justifying stochastic fluid partial equations in the Eulerian formulation. To justify the application of rigorous homogenization theory, we assume mildly chaotic fast small-scale dynamics, as well as a centering condition. The latter requires that the mean of the fluctuating deviations is small, when pulled back to the mean flow., Comment: 11 pages, geometric mechanics, stochastic fluid models, stochastic processes, multi-scale fluid dynamics, symmetry reduced variational principles, homogenisation

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\`ere equation for a scalar mesh potential. We adapt an existing finite element scheme for the standard Monge-Amp\`ere 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\`ere equation due to the implicit dependence of the monitor function on the resulting mesh. We also derive the equivalent Monge-Amp\`ere-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., Comment: Updated following reviews, 36 pages

A vertical slice model is developed for the Euler-Boussinesq equations with a constant temperature gradient in the direction normal to the slice (the Eady-Boussinesq model). The model is a solution of the full three-dimensional equations with no variation normal to the slice, which is an idealized problem used to study the formation and subsequent evolution of weather fronts. A compatible finite element method is used to discretise the governing equations. To extend the Charney-Phillips grid staggering in the compatible finite element framework, we use the same node locations for buoyancy as the vertical part of velocity and apply a transport scheme for a partially continuous finite element space. For the time discretisation, we solve the semi-implicit equations together with an explicit strong-stability-preserving Runge-Kutta scheme to all of the advection terms. The model reproduces several quasi-periodic lifecycles of fronts despite the presence of strong discontinuities. An asymptotic limit analysis based on the semi-geostrophic theory shows that the model solutions are converging to a solution in cross-front geostrophic balance. The results are consistent with the previous results using finite difference methods, indicating that the compatible finite element method is performing as well as finite difference methods for this test problem. We observe dissipation of kinetic energy of the cross-front velocity in the model due to the lack of resolution at the fronts, even though the energy loss is not likely to account for the large gap on the strength of the fronts between the model result and the semi-geostrophic limit solution., Comment: 27 pages, 13 figures. Revised version, accepted for publication in JCP

We analyse the multiscale properties of energy-conserving upwind-stabilised finite element discretisations of the two-dimensional incompressible Euler equations. We focus our attention on two particular methods: the Lie derivative discretisation introduced in Natale and Cotter (2016a) and the Streamline Upwind/Petrov-Galerkin (SUPG) discretisation of the vorticity advection equation. Such discretisations provide control on enstrophy by modelling different types of scale interactions. We quantify the performance of the schemes in reproducing the non-local energy backscatter that characterises two-dimensional turbulent flows.

We propose a finite element discretisation approach for the incompressible Euler equations which mimics their geometric structure and their variational derivation. In particular, we derive a finite element method that arises from a nonholonomic variational principle and an appropriately defined Lagrangian, where finite element H(div) vector fields are identified with advection operators; this is the first successful extension of the structure-preserving discretisation of Pavlov et al. (2009) to the finite element setting. The resulting algorithm coincides with the energy-conserving scheme presented in Guzm\'an et al. (2016). Through the variational derivation, we discover that it also satisfies a discrete analogous of Kelvin's circulation theorem. Further, we propose an upwind-stabilised version of the scheme which dissipates enstrophy whilst preserving energy conservation and the discrete Kelvin's theorem. We prove error estimates for this version of the scheme, and we study its behaviour through numerical tests., Comment: Revision to version published on IMAJNA

Compatible finite elements provide a framework for preserving important structures in equations of geophysical fluid dynamics, and are becoming important in their use for building atmosphere and ocean models. We survey the application of compatible finite element spaces to geophysical fluid dynamics, including the application to the nonlinear rotating shallow water equations, and the three-dimensional compressible Euler equations. We summarise analytic results about dispersion relations and conservation properties, and present new results on approximation properties in three dimensions on the sphere, and on hydrostatic balance properties., Comment: Revision to version published on DSCS

The selective frequency damping (SFD) method is an alternative to classical Newton's method to obtain unstable steady-state solutions of dynamical systems. However this method has two main limitations: it does not converge for arbitrary control parameters; and when it does converge, the time necessary to reach the steady-state solution may be very long. In this paper we present an adaptive algorithm to address these two issues. We show that by evaluating the dominant eigenvalue of a "partially converged" steady flow, we can select a control coefficient and a filter width that ensure an optimum convergence of the SFD method. We apply this adaptive method to several classical test cases of computational fluid dynamics and we show that a steady-state solution can be obtained without any a priori knowledge of the flow stability properties.

We describe and implement a symbolic algebra for scalar and vector-valued finite elements, enabling the computer generation of elements with tensor product structure on quadrilateral, hexahedral and triangular prismatic cells. The algebra is implemented as an extension to the domain-specific language UFL, the Unified Form Language. This allows users to construct many finite element spaces beyond those supported by existing software packages. We have made corresponding extensions to FIAT, the FInite element Automatic Tabulator, to enable numerical tabulation of such spaces. This tabulation is consequently used during the automatic generation of low-level code that carries out local assembly operations, within the wider context of solving finite element problems posed over such function spaces. We have done this work within the code-generation pipeline of the software package Firedrake; we make use of the full Firedrake package to present numerical examples., Comment: Submitted to SISC special issue on CSE Software. Updated version, following reviewer comments

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.

A new numerical method is presented for solving the rotating shallow water equations on a rotating sphere using quasi-uniform polygonal meshes. The method uses special families of finite element function spaces to mimic key mathematical properties of the continuous equations and thereby capture several desirable physical properties related to balance and conservation. The method relies on two novel features. The first is the use of compound finite elements to provide suitable finite element spaces on general polygonal meshes. The second is the use of dual finite element spaces on the dual of the original mesh, along with suitably defined discrete Hodge star operators to map between the primal and dual meshes, enabling the use of a finite volume scheme on the dual mesh to compute potential vorticity fluxes. The resulting method has the same mimetic properties as a finite volume method presented previously, but is more accurate on a number of standard test cases.

This paper presents a family of spatial discretisations of the nonlinear rotating shallow-water equations that conserve both energy and potential enstrophy. These are based on two-dimensional mixed finite element methods and hence, unlike some finite difference methods, do not require an orthogonal grid. Numerical verification of the aforementioned properties is also provided., Comment: 25 pages, 14 figures, submitted to QJRMS; parts of sections 3 and 4 rewritten for clarification, assorted typo fixes

The traditional approach to investigating the stability of a physical system is to linearise the equations about a steady base solution, and to examine the eigenvalues of the linearised operator. Over the past several decades, it has been recognised that this approach only determines the asymptotic stability of the system, and neglects the possibility of transient perturbation growth arising due to the nonnormality of the system. This observation motivated the development of a more powerful generalised stability theory (GST), which focusses instead on the singular value decomposition of the linearised propagator of the system. While GST has had significant successes in understanding the stability of phenomena in geophysical fluid dynamics, its more widespread applicability has been hampered by the fact that computing the SVD requires both the tangent linear operator and its adjoint: deriving the tangent linear and adjoint models is usually a considerable challenge, and manually embedding them inside an eigensolver is laborious. In this paper, we present a framework for the automation of generalised stability theory, which overcomes these difficulties. Given a compact high-level symbolic representation of a finite element discretisation implemented in the FEniCS system, efficient C++ code is automatically generated to assemble the forward, tangent linear and adjoint models; these models are then used to calculate the optimally growing perturbations to the forward model, and their growth rates. By automating the stability computations, we hope to make these powerful tools a more routine part of computational analysis. The efficiency and generality of the framework is demonstrated with applications drawn from geophysical fluid dynamics, phase separation and quantum mechanics., Comment: Accepted in SISC

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.