Abstract. In this paper we demonstrate a new technique for deriving discrete adjoint and tangent linear models of finite element models. 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 code for the derived models. The approach requires only trivial changes to a large class of forward models, including complicated time-dependent nonlinear models. The adjoint model automatically employs optimal checkpointing schemes to mitigate storage requirements for nonlinear models, without any user management or intervention. Furthermore, both the tangent linear and adjoint models naturally work in parallel, without any need to differentiate through calls to MPI or to parse OpenMP directives. The generality, applicability and efficiency of the approach are demonstrated with examples from a wide range of scientific applications. . While deriving the adjoint model associated with a linear stationary forward model is straightforward, the development and implementation of adjoint models for nonlinear or time-dependent forward models is notoriously difficult, for several reasons. First, each nonlinear operator of the forward model must be differentiated, which can be difficult for complex models. Second, the control flow of the adjoint model runs backwards, from the final time to the initial time, and requires access to the solution variables computed during the forward run if the forward problem is nonlinear. Since it is generally impractical for physically relevant simulations to store all variables during the forward run, the adjoint model developer must implement some checkpointing scheme that balances recomputation and storage [17]. The control flow of such a checkpointing scheme must alternate between the solution of forward variables and adjoint variables, and is thus highly

Oceanic tides have the potential to yield a vast amount of renewable energy. Tidal stream generators are one of the key technologies for extracting and harnessing this potential. In order to extract an economically useful amount of power, hundreds of tidal turbines must typically be deployed in an array. This naturally leads to the question of how these turbines should be configured to extract the maximum possible power: the positioning and the individual tuning of the turbines could significantly influence the extracted power, and hence is of major economic interest. However, manual optimisation is difficult due to legal site constraints, nonlinear interactions of the turbine wakes, and the cubic dependence of the power on the flow speed. The novel contribution of this paper is the formulation of this problem as an optimisation problem constrained by a physical model, which is then solved using an efficient gradient-based optimisation algorithm. In each optimisation iteration, a two-dimensional finite element shallow water model predicts the flow and the performance of the current array configuration. The gradient of the power extracted with respect to the turbine positions and their tuning parameters is then computed in a fraction of the time taken for a flow solution by solving the associated adjoint equations. These equations propagate causality backwards through the computation, from the power extracted back to the turbine positions and the tuning parameters. This yields the gradient at a cost almost independent of the number of turbines, which is crucial for any practical application. The utility of the approach is demonstrated by optimising turbine arrays in four idealised scenarios and a more realistic case with up to 256 turbines in the Inner Sound of the Pentland Firth, Scotland.

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 meshrefinement, 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 * Corresponding author Email addresses: a.beskos@ucl.ac.uk (Alexandros Beskos), M.Girolami@warwick.ac.uk (Mark Girolami), slan@caltech.edu (Shiwei Lan), patrick.farrell@maths.ox.ac.uk (Patrick E. Farrell), A.M.Stuart@warwick.ac.uk (Andrew M. Stuart)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. posals which are appropriate for non-linear likelihood-informed subspaces. A simpler step is then developed for a complementary subspace obtained by truncating the Karhunen-Loève expansion of the Gaussian prior. Other such separation methods used in the non-geometric context (likelihood informed subspace [LIS 11] or the active subspace [AS 12]) could potentially be brought into our 50 setting, though this requires further research. Lastly, we apply the geometric methods together with other main MCMC algorithms on three challenging inverse problems and contrast their efficiency. Two elliptic inverse problems, involving a groundwater flow and a thermal fin, aim to infer the coefficients of the elliptic PDEs (representing the permeability of a porous medium and the 55 heat conductivity of a material respectively) from data taken at given locations of the forward solver. The third inverse problem involves an incompressible Navier-Stokes equation, with the objective to infer the inflow vel...

No abstract

Nonlinear systems of partial differential equations (PDEs) may permit several distinct solutions. The typical current approach to finding distinct solutions is to start Newton's 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.

No abstract

SUMMARYThe variational approach to fracture is effective for simulating the nucleation and propagation of complex crack patterns but is computationally demanding. The model is a strongly nonlinear non-convex variational inequality that demands the resolution of small length scales. The current standard algorithm for its solution, alternate minimization, is robust but converges slowly and demands the solution of large, ill-conditioned linear subproblems. In this paper, we propose several advances in the numerical solution of this model that improve its computational efficiency. We reformulate alternate minimization as a nonlinear Gauss-Seidel iteration and employ over-relaxation to accelerate its convergence; we compose this accelerated alternate minimization with Newton's method, to further reduce the time to solution, and we formulate efficient preconditioners for the solution of the linear subproblems arising in both alternate minimization and in Newton's method. We investigate the improvements in efficiency on several examples from the literature; the new solver is five to six times faster on a majority of the test cases considered.

In Benzi & Olshanskii (SIAM J. Sci. Comput., 28(6) (2006)) a preconditioner of augmented Lagrangian type was presented for the two-dimensional stationary incompressible Navier-Stokes equations that exhibits convergence almost independent of Reynolds number. The algorithm relies on a highly specialized multigrid method involving a custom prolongation operator and for robustness requires the use of piecewise constant finite elements for the pressure. However, the prolongation operator and velocity element used do not directly extend to three dimensions: the local solves necessary in the prolongation operator do not satisfy the inf-sup condition. In this work we generalize the preconditioner to three dimensions, proposing alternative finite elements for the velocity and prolongation operators for which the preconditioner works robustly. The solver is effective at high Reynolds number: on a three-dimensional lid-driven cavity problem with approximately one billion degrees of freedom, the average number of Krylov iterations per Newton step varies from 4.5 at Re = 10 to 3 at Re = 1000 and 5 at Re = 5000. , to thank M. A. Olshanskii for supplying the code for the Oseen solver described in [7], and to thank M. G. Knepley for assistance with PETSc.