Patrick E. Farrell

AUTOMATED DERIVATION OF THE ADJOINT OF HIGH-LEVEL TRANSIENT FINITE ELEMENT PROGRAMS

In this paper we demonstrate a new technique for deriving discrete adjoint and tangent linear models of a finite element model. The technique is significantly more efficient and automatic than standard algorithmic differentiation techniques. The approach relies on a high-level symbolic representation of the forward problem. In contrast to developing a model directly in Fortran or C++, high-level systems allow the developer to express the variational problems to be solved in near-mathematical notation. As such, these systems have a key advantage: since the mathematical structure of the problem is preserved, they are more amenable to automated analysis and manipulation. The framework introduced here is implemented in a freely available software package named dolfin-adjoint, based on the FEniCS Project. Our approach to automated adjoint derivation relies on run-time annotation of the temporal structure of the model and employs the FEniCS finite element form compiler to automatically generate the low-level co...

Anisotropic mesh adaptivity for multi-scale ocean modelling

Research into the use of unstructured mesh methods in oceanography has been growing steadily over the past decade. The advantages of this approach for domain representation and non-uniform resolution are clear. However, a number of issues remain, in particular those related to the computational cost of models produced using unstructured mesh methods compared with their structured mesh counterparts. Mesh adaptivity represents an important means to improve the competitiveness of unstructured mesh models, where high resolution is only used when and where necessary. In this paper, an optimization-based approach to mesh adaptivity is described where emphasis is placed on capturing anisotropic solution characteristics. Comparisons are made between the results obtained with uniform isotropic resolution, isotropic adaptive resolution and fully anisotropic adaptive resolution.

Modelling of fluid–solid interactions using an adaptive mesh fluid model coupled with a combined finite–discrete element model

Fluid–structure interactions are modelled by coupling the finite element fluid/ocean model ‘Fluidity-ICOM’ with a combined finite–discrete element solid model ‘Y3D’. Because separate meshes are used for the fluids and solids, the present method is flexible in terms of discretisation schemes used for each material. Also, it can tackle multiple solids impacting on one another, without having ill-posed problems in the resolution of the fluid’s equations. Importantly, the proposed approach ensures that Newton’s third law is satisfied at the discrete level. This is done by first computing the action–reaction force on a supermesh, i.e. a function superspace of the fluid and solid meshes, and then projecting it to both meshes to use it as a source term in the fluid and solid equations. This paper demonstrates the properties of spatial conservation and accuracy of the method for a sphere immersed in a fluid, with prescribed fluid and solid velocities. While spatial conservation is shown to be independent of the mesh resolutions, accuracy requires fine resolutions in both fluid and solid meshes. It is further highlighted that unstructured meshes adapted to the solid concentration field reduce the numerical errors, in comparison with uniformly structured meshes with the same number of elements. The method is verified on flow past a falling sphere. Its potential for ocean applications is further shown through the simulation of vortex-induced vibrations of two cylinders and the flow past two flexible fibres.

Parallel anisotropic mesh adaptivity with dynamic load balancing for cardiac electrophysiology

Abstract Simulations in cardiac electrophysiology generally use very fine meshes and small time steps to resolve highly localized wavefronts. This expense motivates the use of mesh adaptivity, which has been demonstrated to reduce the overall computational load. However, even with mesh adaptivity performing such simulations on a single processor is infeasible. Therefore, the adaptivity algorithm must be parallelised. Rather than modifying the sequential adaptive algorithm, the parallel mesh adaptivity method introduced in this paper focuses on dynamic load balancing in response to the local refinement and coarsening of the mesh. In essence, the mesh partition boundary is perturbed away from mesh regions of high relative error, while also balancing the computational load across processes. The parallel scaling of the method when applied to physiologically realistic heart meshes is shown to be good as long as there are enough mesh nodes to distribute over the available parallel processes. It is shown that the new method is dominated by the cost of the sequential adaptive mesh procedure and that the parallel overhead of inter-process data migration represents only a small fraction of the overall cost.

Hybrid OpenMP/MPI Anisotropic Mesh Smoothing

Abstract Mesh smoothing is an important algorithm for the improvement of element quality in unstructured mesh finite element methods. A new optimisation based mesh smoothing algorithm is presented for anisotropic mesh adaptivity. It is shown that this smoothing kernel is very effective at raising the minimum local quality of the mesh. A number of strategies are employed to reduce the algorithms cost while maintaining its effectiveness in improving overall mesh quality. The method is parallelised using hybrid OpenMP/MPI programming methods, and graph colouring to identify independent sets. Different approaches are explored to achieve good scaling performance within a shared memory compute node.

DEFLATION TECHNIQUES FOR FINDING DISTINCT SOLUTIONS OF NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS

Nonlinear systems of partial differential equations (PDEs) may permit several distinct solutions. The typical current approach to finding distinct solutions is to start Newtons method with many different initial guesses, hoping to find starting points that lie in different basins of attraction. In this paper, we present an infinite-dimensional deflation algorithm for systematically modifying the residual of a nonlinear PDE problem to eliminate known solutions from consideration. This enables the Newton--Kantorovitch iteration to converge to several different solutions, even starting from the same initial guess. The deflated Jacobian is dense, but an efficient preconditioning strategy is devised, and the number of Krylov iterations is observed not to grow as solutions are deflated. The power of the approach is demonstrated on several problems from special functions, phase separation, differential geometry and fluid mechanics that permit distinct solutions.

Geometric MCMC for infinite-dimensional inverse problems

Bayesian inverse problems often involve sampling posterior distributions on infinite-dimensional function spaces. Traditional Markov chain Monte Carlo (MCMC) algorithms are characterized by deteriorating mixing times upon mesh-refinement, when the finite-dimensional approximations become more accurate. Such methods are typically forced to reduce step-sizes as the discretization gets finer, and thus are expensive as a function of dimension. Recently, a new class of MCMC methods with mesh-independent convergence times has emerged. However, few of them take into account the geometry of the posterior informed by the data. At the same time, recently developed geometric MCMC algorithms have been found to be powerful in exploring complicated distributions that deviate significantly from elliptic Gaussian laws, but are in general computationally intractable for models defined in infinite dimensions. In this work, we combine geometric methods on a finite-dimensional subspace with mesh-independent infinite-dimensional approaches. Our objective is to speed up MCMC mixing times, without significantly increasing the computational cost per step (for instance, in comparison with the vanilla preconditioned Crank–Nicolson (pCN) method). This is achieved by using ideas from geometric MCMC to probe the complex structure of an intrinsic finite-dimensional subspace where most data information concentrates, while retaining robust mixing times as the dimension grows by using pCN-like methods in the complementary subspace. The resulting algorithms are demonstrated in the context of three challenging inverse problems arising in subsurface flow, heat conduction and incompressible flow control. The algorithms exhibit up to two orders of magnitude improvement in sampling efficiency when compared with the pCN method.

A Framework for the Automation of Generalized Stability Theory

The traditional approach to investigating the stability of a physical system is to linearize the equations about a steady base solution, and to examine the eigenvalues of the linearized operator. Over the past several decades, it has been recognized 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 generalized stability theory (GST), which focuses instead on the singular value decomposition (SVD) of the linearized 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 ...

Directional integration on unstructured meshes via supermesh construction

Unstructured meshes are in widespread use throughout computational physics, but calculating diagnostics of simulations on such meshes can be challenging. For example, in geophysical fluid dynamics, it is frequently desirable to compute directional integrals such as vertical integrals and zonal averages; however, it is difficult to compute these on meshes with no inherent spatial structure. This is widely regarded as an obstacle to the adoption of unstructured mesh numerical modelling in this field. In this paper, we describe an algorithm by which one can exactly compute such directional integrals on arbitrarily unstructured meshes. This is achieved via the solution of a problem of computational geometry, constructing the supermesh of two meshes. We demonstrate the utility of this approach by applying it to a classical geophysical fluid dynamics system: the thermally driven rotating annulus. This addresses an important objection to the more widespread use of unstructured mesh modelling.

Geostrophic balance preserving interpolation in mesh adaptive linearised shallow-water ocean modelling

Abstract The accurate representation of geostrophic balance is an essential requirement for numerical modelling of geophysical flows. Significant effort is often put into the selection of accurate or optimal balance representation by the discretisation of the fundamental equations. The issue of accurate balance representation is particularly challenging when applying dynamic mesh adaptivity, where there is potential for additional imbalance injection when interpolating to new, optimised meshes. In the context of shallow-water modelling, we present a new method for preservation of geostrophic balance when applying dynamic mesh adaptivity. This approach is based upon interpolation of the Helmholtz decomposition of the Coriolis acceleration. We apply this in combination with a discretisation for which states in geostrophic balance are exactly steady solutions of the linearised equations on an f-plane; this method guarantees that a balanced and steady flow on a donor mesh remains balanced and steady after interpolation onto an arbitrary target mesh, to within machine precision. We further demonstrate the utility of this interpolant for states close to geostrophic balance, and show that it prevents pollution of the resulting solutions by imbalanced perturbations introduced by the interpolation.