Robert C. Kirby

Robert C. Kirby, Anders Logg

As a key step towards a complete automation of the finite element method, we present a new algorithm for automatic and efficient evaluation of multilinear variational forms. The algorithm has been implemented in the form of a compiler, the FEniCS Form Compiler FFC. We present benchmark results for a series of standard variational forms, including the incompressible Navier-Stokes equations and linear elasticity. The speedup compared to the standard quadrature-based approach is impressive; in some cases the speedup is as large as a factor 1000.

Marie Rognes, Robert C. Kirby, Anders Logg

In this paper, we discuss how to efficiently evaluate and assemble general finite element variational forms on H(div) and H(curl). The proposed strategy relies on a decomposition of the element tensor into a precomputable reference tensor and a mesh-dependent geometry tensor. Two key points must then be considered: the appropriate mapping of basis functions from a reference element, and the orientation of geometrical entities. To address these issues, we extend here a previously presented representation theorem for affinely mapped elements to Piola-mapped elements. We also discuss a simple numbering strategy that removes the need to contend with directions of facet normals and tangents. The result is an automated, efficient, and easy-to-use implementation that allows a user to specify finite element variational forms on H(div) and H(curl) in close to mathematical notation.

Robert C. Kirby, Matthew Knepley, Anders Logg et al.

Assembling stiffness matrices represents a significant cost in many finite element computations. We address the question of optimizing the evaluation of these matrices. By finding redundant computations, we are able to significantly reduce the cost of building local stiffness matrices for the Laplace operator and for the trilinear form for Navier-Stokes. For the Laplace operator in two space dimensions, we have developed a heuristic graph algorithm that searches for such redundancies and generates code for computing the local stiffness matrices. Up to cubics, we are able to build the stiffness matrix on any triangle in less than one multiply-add pair per entry. Up to sixth degree, we can do it in less than about two. Preliminary low-degree results for Poisson and Navier-Stokes operators in three dimensions are also promising.

Robert C. Kirby, Anders Logg

We investigate the compilation of general multilinear variational forms over affines simplices and prove a representation theorem for the representation of the element tensor (element stiffness matrix) as the contraction of a constant reference tensor and a geometry tensor that accounts for geometry and variable coefficients. Based on this representation theorem, we design an algorithm for efficient pretabulation of the reference tensor. The new algorithm has been implemented in the FEniCS Form Compiler (FFC) and improves on a previous loop-based implementation by several orders of magnitude, thus shortening compile-times and development cycles for users of FFC.

Robert C. Kirby, Anders Logg, L. Ridgway Scott et al.

We present a topological framework for finding low-flop algorithms for evaluating element stiffness matrices associated with multilinear forms for finite element methods posed over straight-sided affine domains. This framework relies on phrasing the computation on each element as the contraction of each collection of reference element tensors with an element-specific geometric tensor. We then present a new concept of complexity-reducing relations that serve as distance relations between these reference element tensors. This notion sets up a graph-theoretic context in which we may find an optimized algorithm by computing a minimum spanning tree. We present experimental results for some common multilinear forms showing significant reductions in operation count and also discuss some efficient algorithms for building the graph we use for the optimization.

Robert C. Kirby, Anders Logg

We examine the effect of using complexity-reducing relations to generate optimized code for the evaluation of finite element variational forms. The optimizations are implemented in a prototype code named FErari, which has been integrated as an optimizing backend to the FEniCS Form Compiler, FFC. In some cases, FErari provides very little speedup, while in other cases, we obtain reduced local operation counts of a factor of as much as 7.9 and speedups for the assembly of the global sparse matrix of as much as a factor of 2.8.

Miklós Homolya, Robert C. Kirby, David A. Ham

Code generation based software platforms, such as Firedrake, have become popular tools for developing complicated finite element discretisations of partial differential equations. We extended the code generation infrastructure in Firedrake with optimisations that can exploit the structure inherent to some finite elements. This includes sum factorisation on cuboid cells for continuous, discontinuous, H(div) and H(curl) conforming elements. Our experiments confirm optimal algorithmic complexity for high-order finite element assembly. This is achieved through several novel contributions: the introduction of a more powerful interface between the form compiler and the library providing the finite elements; a more abstract, smarter library of finite elements called FInAT that explicitly communicates the structure of elements; and form compiler algorithms to automatically exploit this exposed structure.

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.

Pablo D. Brubeck, Patrick E. Farrell, Robert C. Kirby et al.

We present new finite elements for solving the Riesz maps of the de Rham complex on triangular and tetrahedral meshes at high order. The finite elements discretize the same spaces as usual, but with different basis functions, so that the resulting matrices have desirable properties. These properties mean that we can solve the Riesz maps to a given accuracy in a $p$-robust number of iterations with $\mathcal{O}(p^6)$ flops in three dimensions, rather than the naïve $\mathcal{O}(p^9)$ flops. The degrees of freedom build upon an idea of Demkowicz et al., and consist of integral moments on an equilateral reference simplex with respect to a numerically computed polynomial basis that is orthogonal in two different inner products. As a result, the interior-interface and interior-interior couplings are provably weak, and we devise a preconditioning strategy by neglecting them. The combination of this approach with a space decomposition method on vertex and edge star patches allows us to efficiently solve the canonical Riesz maps at high order. We apply this to solving the Hodge Laplacians of the de Rham complex with novel augmented Lagrangian preconditioners.

Robert C. Kirby

The use of a reference element on which a finite element basis is constructed once and mapped to each cell in a mesh greatly expedites the structure and efficiency of finite element codes. However, many famous finite elements such as Hermite, Morley, Argyris, and Bell, do not possess the kind of equivalence needed to work with a reference element in the standard way. This paper gives a generalizated approach to mapping bases for such finite elements by means of studying relationships between the finite element nodes under push-forward.

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.

Robert C. Kirby

We consider the discontinuous Galerkin method for hyperbolic conservation laws, with some particular attention to the linear acoustic equation, using Bernstein polynomials as local bases. Adapting existing techniques leads to optimal-complexity computation of the element and boundary flux terms. The element mass matrix, however, requires special care. In particular, we give an explicit formula for its eigenvalues and exact characterization of the eigenspaces in terms of the Bernstein representation of orthogonal polynomials. We also show a fast algorithm for solving linear systems involving the element mass matrix to preserve the overall complexity of the DG method. Finally, we present numerical results investigating the accuracy of the mass inversion algorithms and the scaling of total run-time for the function evaluation needed in DG time-stepping.