PaperlayerView exact revision on arXiv

Source-linked AI summary

Firedrake: automating the finite element method by composing abstractions

Florian Rathgeber, David A. Ham, Lawrence Mitchell, Michael Lange, Fabio Luporini, Andrew T. T. McRae, Gheorghe-Teodor Bercea, Graham R. Markall, Paul H. J. Kelly

arXiv:1501.01809v3cs.MSmath.NA

TL;DR

Finite element software must reconcile high-level mathematical descriptions with fine-grained, performance-critical parallel implementations. Firedrake composes abstractions, especially PyOP2, to separate local discretisation from execution while retaining a Python-based finite element interface. Its assembly is consistently faster than DOLFIN, with COFFEE delivering up to a fourfold speed increase over FFC quadrature optimisations.

  • Problem

    Finite element implementations must express innermost mesh operations where loop order, data layout, and fine-grained parallelism are performance-critical.

  • Method

    Firedrake composes abstractions and uses PyOP2 to separate local finite element discretisation from parallel execution over the mesh.

  • Results

    Up to a fourfold speed increase over FFC quadrature optimisations is reported for COFFEE, with assembly consistently faster than DOLFIN.

  • Takeaways & Limitations

    The abstraction composition separates contributions to finite elements and parallel execution while supporting automated performance optimisations.

  • Takeaways & Limitations

    The performance evaluation is intended to give a broad impression rather than provide a comprehensive evaluation of a single problem.

Abstract

from arXiv · show

Firedrake is a new tool for automating the numerical solution of partial differential equations. Firedrake adopts the domain-specific language for the finite element method of the FEniCS project, but with a pure Python runtime-only implementation centred on the composition of several existing and new abstractions for particular aspects of scientific computing. The result is a more complete separation of concerns which eases the incorporation of separate contributions from computer scientists, numerical analysts and application specialists. These contributions may add functionality, or improve performance. Firedrake benefits from automatically applying new optimisations. This includes factorising mixed function spaces, transforming and vectorising inner loops, and intrinsically supporting block matrix operations. Importantly, Firedrake presents a simple public API for escaping the UFL abstraction. This allows users to implement common operations that fall outside pure variational formulations, such as flux-limiters.

1. INTRODUCTION

Firedrake addresses the multidisciplinary complexity of finite element software by composing abstractions that separate mathematical formulation, implementation, and execution. Its PyOP2 layer separates local discretisation from parallel mesh execution, while the resulting codebase remains compact.

  • Finite element software combines expertise in PDEs, discretisation, meshes, solvers, parallel algorithms, vectorisation, and loop optimisation.
  • Abstraction and composition let specialists collaborate without understanding every implementation layer in full detail.
  • Firedrake introduces PyOP2 to separate local mathematical operators from their parallel execution over the mesh.
  • This separation supports contributions from numerical analysts developing finite elements and computer scientists developing execution strategies.
  • Around 5000 executable lines comprise core Firedrake, while PyOP2 contains fewer than 9000 executable lines.The paper presents this compactness as benefiting maintainability and extensibility.
  • The paper evaluates Firedrake’s performance and capability through an extensive computational verification.

2. MATHEMATICAL AND SOFTWARE ABSTRACTION OF THE FINITE ELEMENT METHOD

Finite element methods have a high-level mathematical description, but efficient implementations expose performance-critical low-level details. Existing library and DSL approaches preserve parts of that abstraction while facing the fine-grained variability of PDE discretisation.

  • Finite element algorithms can often be specified by a weak-form PDE, boundary conditions, and discrete function spaces.
  • Hand-coded low-level implementations replace symbolic equations and integrals with array loops, communication calls, threading, and vectorisation directives.
  • Low-level implementations commit discretisation and equation choices early, mix mathematical and implementation concerns, and risk bugs when changed.
  • Object-oriented scientific libraries preserve higher-level mathematical objects while hiding primitive floating-point operations inside methods.
  • Library interfaces work well for coarse-grained operations, but PDE discretisation varies at innermost mesh loops where layout and vectorisation affect performance.
  • Dune-FEM and Deal.II therefore require application developers to provide low-level C++ integral implementations over quadrature points.
  • DSLs retain symbolic variational specifications while allowing compilers to choose optimised implementations and automate mathematical reasoning.

3. EXPLOITING COMPOSABLE ABSTRACTIONS IN FIREDRAKE

Firedrake composes domain-specific abstractions into a Python finite element toolchain, with PyOP2 providing a uniform interface for mesh iteration and parallel execution. The design retains mathematical problem specification while enabling runtime code generation and optimisations.

  • PyOP2 uniformly specifies mesh iterations, decoupling local computation from execution over the whole domain.
  • Firedrake composes separate abstract processes and packages into a largely seamless finite element abstraction.
  • Firedrake adopts UFL and a mostly compatible DOLFIN Python interface, while adding language extensions and omitting some DOLFIN features.
  • A finite element problem requires more than PDE and element specifications: it also needs a mesh, field values, forcing functions, and solve ordering.
  • PyOP2 dynamically generates runtime code by inspecting Python objects and data structures, and supplies finite-element sparse matrices.

4. PYOP2

PyOP2 represents unstructured-mesh computation through sets, maps, data objects, and parallel loops. Access descriptors expose enough information for automated scheduling, communication, and contention avoidance while preserving a uniform kernel abstraction.

  • Mesh-based numerical algorithms apply local computational kernels independently across mesh entities.
  • PyOP2 models meshes as graphs of entity sets and maps connecting source entities to target data.
  • PyOP2 does not assign semantic meaning to meshes or function spaces; Firedrake interprets its sets and maps.
  • Constant-arity maps fix loop bounds, enabling vectorisation and related optimisations, while excluding some irregular mappings.
  • Dats represent discretised vectors, Mats sparse matrices, and Globals data independent of individual set members.
  • A Parloop applies a kernel over an iteration set using access descriptors and indirection maps for its data arguments.
  • PyOP2 permits unordered parallel execution across threads, vector lanes, or processes and uses descriptors to manage indirect-access write contention.
  • Deterministic colouring provides bit-reproducible results when the same number of processors is used.The paper leaves sustainability under increasingly fine-grained hardware parallelism unresolved.

4.5. Kernel optimisation in COFFEE

COFFEE optimises finite element assembly kernels represented as abstract syntax trees, targeting fewer floating-point operations and greater instruction-level parallelism through SIMD vectorisation.

  • Kernel representation: Kernels accept C code strings or abstract syntax trees, enabling COFFEE to optimise short loops around non-trivial mathematical expressions.The AST is used to generate C code and exposes the kernel structure to the optimiser.
  • Optimisation goals: COFFEE minimises floating-point operations and improves instruction-level parallelism through platform-specific AST transformations.Its transformations include invariant-subexpression hoisting, loop permutation, loop unrolling, and expression vectorisation.
  • Optimisation goals: SIMD vectorisation may require padding data and enforcing alignment constraints for the target architecture.

5. THE FIREDRAKE LAYER

The Firedrake layer composes UFL, FIAT, FFC, PETSc, and PyOP2 abstractions to turn high-level finite element specifications into efficient parallel computations. It also supports richer mesh, function-space, mixed-system, and non-variational workflows.

  • Layer composition: Firedrake marshals UFL, FIAT, FFC, PETSc, and PyOP2 to efficiently produce solutions from finite element problems specified in the FEniCS language.
  • Layer composition: Firedrake implements FEniCS-language operations by selecting PyOP2 objects and composing parallel loops, while PETSc supplies complementary objects.
  • Mesh abstraction: Coordinates are represented as first-class vector-valued fields, allowing field operations on geometry and supporting curved elements.A fully featured function space also provides a mechanism for isoparametric elements, available in a Firedrake branch.
  • Function spaces and functions: Function and function-space objects map naturally onto PyOP2 data types, delegating data storage and communication to PyOP2.Function spaces contain maps from mesh entities to degrees of freedom, while functions hold coefficient data.
  • Assembly: Variational assembly extracts iteration sets, maps, coordinate fields, coefficients, and output tensors from the problem’s function spaces.Bilinear forms produce PETSc matrices, while linear forms produce data objects associated with the test space.
  • Modified FFC: Modified FFC produces unscheduled ASTs that Firedrake converts into PyOP2 kernels and parallel loops, preserving PyOP2’s separation between kernel specification and optimisation.The Firedrake-specific FFC is effectively a fork, although UFLACS was expected to support a unified compiler infrastructure.
  • Escaping the abstraction: Firedrake enables graceful escape from UFL by allowing custom C or AST kernels for mesh-local operations outside pure variational formulations.These kernels can be explicitly executed over the mesh through parallel loops.
  • Mixed function spaces: Mixed forms are symbolically split into blocks, generating simpler kernels and hierarchical PETSc nested matrices that support block solver techniques such as Schur complements.Smaller submatrices are expected to accelerate sparse insertion by reducing the number of non-zero columns searched per row.

6. EXPERIMENTS

The experiments assess Firedrake across several finite element problems and implementation characteristics, while explicitly presenting broad performance impressions rather than a comprehensive evaluation.

  • Problems and capabilities: Experiments cover stationary Poisson, nonlinear time-dependent Cahn–Hilliard, and explicitly time-stepped linear wave equations.They examine regular and mixed finite element assembly among other implementation characteristics.
  • Evaluation scope: The evaluation is designed to provide an impression of Firedrake’s broad performance characteristics, not a comprehensive performance evaluation.The authors note that a comprehensive evaluation of even one problem could occupy an entire paper.
  • Comparison: DOLFIN is the comparison package because it is Firedrake’s closest analogue and supports nearly identical test cases from similar code.The setup uses developer-recommended compilation flags and DOLFIN-based examples to reduce risks from inexpert use.
  • Reproducibility: Benchmark source code and driver scripts are available in the firedrake-bench repository, with the experimental version archived on Zenodo.

6.1. Experimental setup

Experiments ran on ARCHER with specified compilers, libraries, optimisations, mesh reorderings, and measurement controls. The protocol used cache-warming dry runs and reported the minimum of three consecutive timings.

  • Hardware: Experiments used ARCHER, a Cray XC30 supercomputer with two 12-core Intel Xeon processors and 64GB memory per node.The machine used an Aries Dragonfly interconnect and NUMA-linked processor memory regions.
  • Software: Firedrake and PETSc used GNU compilers 4.9.2 with Cray MPICH2 7.1.1 and asynchronous progress enabled for parallel runs.
  • Compilation: Firedrake generated code used -O3 -fno-tree-vectorize, whereas DOLFIN used -O3 -ffast-math -march=native.
  • Compilation: Intel and Cray compilers were unavailable to PyOP2’s runtime compilation because of technical limitations accessing the licence server.
  • Runtime configuration: Both configurations used quadrature representation and mesh reordering, with Firedrake additionally applying COFFEE loop-invariant motion, alignment, and padding optimisations.
  • Measurement protocol: Runs used exclusive node access and process pinning, followed a cache-warming dry run, and reported the minimum of three consecutive measurements.

6.2. Poisson

The Poisson benchmark evaluates Firedrake’s assembly and solver performance against DOLFIN across polynomial degrees and core counts. Firedrake assembles efficiently and maintains strong matrix-assembly efficiency, while solver costs dominate and scaling eventually flattens.

  • Problem setup: The benchmark assembles bilinear and linear forms into a sparse matrix and vector, then solves the linear system with preconditioned Krylov methods.It uses CG with HYPRE BoomerAMG on unit-cube meshes of varying resolution and polynomial degree.
  • Strong scaling: Strong-scaling runtimes compare Firedrake and DOLFIN for assembly and solves on up to 1536 cores and approximately 0.5M–14M degrees of freedom.The comparison covers first- through third-order basis functions.
  • Strong scaling: Solve time dominates all cases, especially at higher order and near the strong-scaling limit, where scaling flattens around 10k DOFs per core.Firedrake assembles left- and right-hand sides faster in almost all cases, but right-hand-side overhead appears sooner because its sequential baseline is faster.
  • Strong scaling: Firedrake matrix assembly shows the highest strong-scaling efficiency across degree one, two, and three basis functions.Right-hand-side efficiency tails off relative to DOLFIN, while solver efficiencies are almost identical with a slight Firedrake advantage at third order.
  • Weak scaling: With 50k DOFs per core, Firedrake has better assembly efficiency beyond one node, although DOLFIN remains faster overall for right-hand-side assembly.Within a node, shared-resource contention and Firedrake’s faster sequential baseline reduce its assembly efficiency.

6.3. Linear Wave Equation

The linear wave benchmark uses explicit time stepping with pointwise updates enabled by mass lumping, avoiding linear-system solves. Firedrake’s direct update scales close to its overhead model, while the indirect update is limited by communication at small core counts and weak scaling remains strong.

  • Method: An explicit symplectic scheme offsets p and φ by half a timestep, and mass lumping converts the p update’s mass-matrix inversion into pointwise multiplication.The φ update is already pointwise because its equation has no spatial derivatives.
  • Method: The benchmark uses Firedrake’s expression compiler for p and φ updates without solving a linear system or invoking a PETSc solver.Symbolic expressions generate pointwise calculations when assignment and accumulation operations are called.
  • Scaling scope: The wave benchmark’s strong scaling is limited by measured non-parallelisable overhead, while weak-scaling tests cover up to 24 intra-node and 384 inter-node cores.Timings for the weak-scaling figure cover 100 timesteps.
  • Strong scaling: The direct-loop φ update follows its projected strong-scaling curve almost perfectly, whereas the indirect-loop p update trails because halo-exchange communication begins at 3 cores.Caching assembled expressions keeps sequential overheads low.
  • Weak scaling: For 84k DOFs per core, both updates retain similarly high intra-node weak-scaling efficiency, declining to about 80% at 24 cores.Across nodes, scaling is almost perfect for both p and φ updates over 24–384 cores.

6.4. Cahn-Hilliard Equation

The Cahn–Hilliard experiment reformulates a fourth-order nonlinear equation as two coupled second-order equations and solves the resulting mixed system with a fieldsplit-preconditioned GMRES method.

  • Problem formulation: The Cahn–Hilliard problem models phase separation in a binary fluid and contains first-order time derivatives with second- and fourth-order spatial derivatives.Introducing the chemical potential μ restates it as two coupled second-order equations.
  • Problem formulation: The variational problem uses unknown concentration c and chemical potential μ in a shared suitable function space V.The formulation is time-discretized with Crank–Nicolson.
  • Experimental setup: The experiment uses a fully unstructured unit-square mesh, first-order Lagrange basis functions, λ = 0.01, M = 1, and dt = 5 · 10^-6.The nonlinear function is f = 100c^2(1 − c^2).
  • Implementation: Firedrake sets the random initial concentration through a custom Kernel executed in a parallel loop, with equivalent lower-level PyOP2 code provided for comparison.The setup seeds randomness by MPI rank and writes the initial value through a kernel.
  • Solver: The mixed system is solved with GMRES and a fieldsplit preconditioner using a lower Schur complement factorisation.An approximate A^-1 uses one HYPRE BoomerAMG V-cycle, while the inverse Schur complement is approximated with a custom PETSc matrix preconditioner.

6.5. Performance discussion

Firedrake’s performance advantages arise from PyOP2 and mixed-space optimisations, while fixed per-field overheads limit strong scaling and PETSc MATNEST improves solver behavior.

  • Scaling performance: In the Cahn–Hilliard strong-scaling experiment, Firedrake is about twice as fast for assembly and about two orders of magnitude faster for initial-condition evaluation.Both systems scale close to linearly for assembly down to 10k DOFs per core.
  • Solver performance: Firedrake’s solver scaling gains significantly from about 80k DOFs per core because PETSc MATNEST avoids costly sub-block extraction copies for fieldsplit preconditioning.The same MATNEST approach avoids the memory allocation and deallocation burden affecting DOLFIN’s monolithic preconditioner.
  • Assembly: Firedrake assembly is consistently much faster than DOLFIN because PyOP2 kernels can be completely inlined, unlike DOLFIN’s repeated C++ virtual function calls.COFFEE optimisations provide up to a fourfold speed increase over FFC quadrature optimisations, with form splitting especially benefiting mixed spaces.
  • Scaling performance: Weak scaling of pure Firedrake code beyond one node is uniformly excellent, although within-node resource contention reduces efficiency.In strong scaling, fixed overheads of some hundreds of microseconds per field cause loss of optimal scaling at higher DOF counts.
  • Scaling performance: Firedrake is about twice as fast for weak-scaling assembly and solve, and almost two orders of magnitude faster for initial-condition evaluation.Inter-node assembly scaling is almost perfect and even superlinear for Firedrake, while initial-condition efficiency stabilises just below 70%.

7. CURRENT LIMITATIONS AND FUTURE EXTENSIONS

The paper documents scope boundaries for Firedrake’s production-quality backends and identifies hybrid parallelism and p-refined finite elements as future extensions.

  • Accelerators and threads: Only the MPI CPU backend with COFFEE-assisted instruction-level vector parallelism is presented as production quality and suitable for full documentation.Other backends support assembly and solving only within the limits of their available solver libraries.
  • Accelerators and threads: Hybrid parallel approaches combining message passing with shared memory are identified as a future direction for increasingly fine-grained CPU and accelerator hardware.The limitation reflects the more restricted hybrid and GPU linear-solver libraries compared with PETSc’s MPI-only functionality.
  • hp-adaptive finite element methods: Supporting p-refined finite elements requires lifting PyOP2’s fixed-arity map restriction.The paper proposes variable-arity maps or container maps holding maps of different arity and corresponding kernels.
  • hp-adaptive finite element methods: Container maps could support both p-refinement and mixed-geometry meshes, but variable-arity loops would interfere with some COFFEE low-level optimisations.The map storage format would also need to record each source element’s arity.
  • Related extensions: The paper relates Firedrake to dolfin-adjoint through the shared mathematical abstraction captured by the FEniCS language.A short Python wrapper extended dolfin-adjoint to support Firedrake solvers written in that language.
Loading 1501.01809v3…