Automatic symbolic computation for discontinuous Galerkin finite element methods

04/06/2018 ∙ by Paul Houston, et al. ∙ The University of Nottingham 0

The implementation of discontinuous Galerkin finite element methods (DGFEMs) represents a very challenging computational task, particularly for systems of coupled nonlinear PDEs, including multiphysics problems, whose parameters may consist of power series or functionals of the solution variables. Thereby, the exploitation of symbolic algebra to express a given DGFEM approximation of a PDE problem within a high level language, whose syntax closely resembles the mathematical definition, is an invaluable tool. Indeed, this then facilitates the automatic assembly of the resulting system of (nonlinear) equations, as well as the computation of Fréchet derivative(s) of the DGFEM scheme, needed, for example, within a Newton-type solver. However, even exploiting symbolic algebra, the discretisation of coupled systems of PDEs can still be extremely verbose and hard to debug. Thereby, in this article we develop a further layer of abstraction by designing a class structure for the automatic computation of DGFEM formulations. This work has been implemented within the FEniCS package, based on exploiting the Unified Form Language. Numerical examples are presented which highlight the simplicity of implementation of DGFEMs for the numerical approximation of a range of PDE problems.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The finite element method (FEM) represents an indispensable computational tool for the accurate, efficient, and rigorous numerical approximation of continuum models arising within a wide range of scientific and engineering application areas. Key reasons for the success of FEMs include their applicability to very general classes of partial differential equations (PDEs), simple treatment of complicated computational geometries and enforcement of boundary conditions, ease of adaptivity including both local mesh subdivision (

–refinement) and local polynomial enrichment (–refinement), and, from a mathematical point of view, the availability of tools for their rigorous error analysis. However, when compared with their finite difference counterparts, FEMs are typically regarded as being complicated

to implement. Indeed, assembly of the underlying matrix stemming from the FEM discretisation of a given (linear, for example) PDE problem typically involves mapping each element present in the computational mesh, defined in the global coordinate system, to a given reference or canonical element, where both the local FEM basis and corresponding quadrature, needed to approximate the underlying integral, is defined. In this manner, local elemental stiffness matrices and load vectors may be computed; these entries are then inserted into the global matrix and right-hand side vector, respectively, according to the elementwise local-to-global degree of freedom mapping, subject to the enforcement of inter-element continuity constraints and the imposition of boundary conditions. This general assembly strategy is almost universally employed within both open source and commercial software, e.g., FreeFem++ 

[23], DUNE [8], deal.II [4] and OpenFOAM [40], to name just a few. Thereby, in principle, assuming that a user can define both the element stiffness matrix and load vector, then a numerical approximation may be readily determined. However, most open source software packages do not provide easy-to-use user interfaces; indeed, typically the user must first understand the low level language in which the code is written, e.g., C++, Fortran, etc, and develop an understanding of the underlying datastructures and function/subroutine calls defined within their chosen package in order to be able to develop code specific to their own application. This can, of course, be a rather time-consuming exercise, and often requires one to use multiple software libraries, for example, when certain features are not available within the chosen package which the user needs to utilise.

In the nonlinear setting, assuming a Newton-type iteration is exploited to solve the underlying system of nonlinear equations stemming from the given FEM employed, the general strategy of assembly of the resulting linearised equations is similar to the linear case; in this setting the load vector represents the residual of the numerical scheme. However, in this case, the Fréchet derivative of the FEM must now be computed; we point out that, in the context of the numerical detection of bifurcation points, the number of derivatives that must be computed just to form the so-called extended system needed to accurately compute the bifurcation point, even before a Newton iteration is implemented, is dependent on the codimension of the singularity being sought, cf. [10, 11]. In general, the evaluation and implementation of both the FEM residual and, moreover, the Fréchet derivative(s) of the scheme is a difficult and time consuming task, which is inevitably prone to human error. This is particularly the case when exploiting FEMs for the numerical approximation of systems of coupled nonlinear PDEs, including multiphysics problems, whose parameters may consist of power series or functionals of the solution variables, cf. [28].

Thus far, our discussion has primarily focussed on the application of conforming FEMs, whereby, excluding Neumann, or weakly imposed boundary conditions, only element contributions need to be evaluated. The exploitation of more general FEMs, and in particular discontinuous Galerkin FEMs (DGFEMs) requires the implementation of inter-element flux terms, which involves combinations of both inner– and outer–traces of the FEM solution, relative to a given orientation of the element face. In recent years, DGFEMs have become an increasingly popular class of FEMs, most notably due to their local conservation properties, inherent numerical stability for convection–dominated diffusion problems, limited interelement communication, which is restricted only to neighbouring elements, and has important advantages for the implementation of boundary conditions and the parallel efficiency of the method, and finally the ease in which so–called hanging nodes can be treated, and the efficient implementation of –adaptivity. Indeed, tremendous progress has been made on both the analytical and computational aspects of DGFEMs; for a review of some of the key developments in the subject, we refer to the recent monographs [13, 14, 24, 35]. For a historical review of DGFEMs, we refer to the articles [3, 12], and the references cited therein.

The addition of inter-element face terms within DGFEMs further complicates the implementation of such schemes within standard FEM packages. Assuming for a moment that the underlying PDE problem is linear, we note that for a given interior face, shared by two neighbouring elements, there are, in general, four local face matrices, which stem from the different combinations of traces of the FEM solution from the two elements whose boundaries form the current face. Once these local face–wise matrices have been assembled, they can then be inserted into the global FEM matrix in an analogous manner to the treatment of the element stiffness matrix. On boundary faces, contributions to the load vector must also be computed. Again, in the nonlinear setting, the task of computing both the residual vector and Fréchet derivative(s) of the DGFEM scheme is a very technical and time consuming task.

The purpose of this article is to discuss the use of symbolic algebra to facilitate the assembly of FEM matrix problems, and in particular those arising from the application of DGFEMs, for the numerical approximation of general nonlinear systems of PDEs. The general approach is to develop a high level language syntax, which closely corresponds to the mathematical formulation of the underlying FEM. Thereby, through this layer of abstraction, the user needs to only specify the FEM residual in a concise and easy to read manner, whereby the evaluation of the Fréchet derivative(s) of the scheme are automatically computed symbolically and the resulting low level C++/Fortran code for element stiffness and face matrices and load/residual vectors are automatically generated by the so-called form compiler. In this way, any existing open source FEM package may be utilised, subject to the implementation of a suitable interface which directly calls the automatically generated snippets of code when assembling element and face matrices and load/residual vectors. Most importantly, we stress that once the user has selected a particular FEM package which is appropriate for their purposes, the interface to that software platform which links to the automatically generated code only needs to be written and debugged once; it may then subsequently be exploited to solve a potentially huge variety of PDE problems, with a plethora of FEM schemes. At this point it is pertinent to mention that some FEM packages do indeed include such a symbolic interface at the heart of their design; most notably we mention the excellent FEniCS package, cf. [1, 31], for example. The Unified Form Language (UFL) component of FEniCS provides an easy to use python interface which allows for the automatic FEM numerical approximation of systems of PDEs in a user-friendly manner. Other form compilers include SyFi [2], and Manycore Form Compiler [32], for example.

However, even exploiting a package such as UFL, as powerful as it is, the definition of DGFEMs in this framework is still rather verbose, particularly for nonlinear systems of coupled PDEs. With this in mind, we present a further layer of abstraction for the automatic computation of DGFEM formulations employing symbolic algebra. Indeed, in this article DGFEM utility functions for general PDE operators are developed, which significantly simplifies the specification of the DGFEM discretisation of a given problem. This work was originally inspired by the need to numerically model the formation of a hydrogen plasma in a microwave power assisted chemical vapour deposition reactor employed for the manufacture of synthetic diamond; this includes constituent equations for the background gas mass average velocity, gas temperature, electromagnetic field energy and plasma density, cf. [28, 36].

In [36] we originally developed easy-to-use DGFEM utility functions as part of our inhouse software package AptoPy [36], for application with our own FEM package AptoFEM [26]. AptoPy is written in Python and exploits the open source symbolic algebra package SymPy [37]. We stress that these choices are entirely user-dependent; indeed, in the past we have employed the symbolic algebra packages REDUCE [22], Mathematica, and Maple. In addition to AptoFEM, ENTWIFE has also been employed as the low level FEM package. The DGFEM utility functions written in AptoPy have been ported to UFL for use within the FEniCS package; full open source codes are available from With this in mind, for simplicity of presentation, throughout this article, we only show snippets of UFL code in order to highlight how the DGFEM utility functions may be exploited in practice; the corresponding AptoPy syntax is quite similar, cf. [36].

The outline of this article is as follows. In Section 2 we briefly outline the DGFEM discretisation of general systems of nonlinear conservation laws. Then, in Section 3 we propose a computational framework for the automatic generation of DGFEM schemes within a simple unified setting. On the basis of this work, in Section 4 we provide some examples to illustrate the flexibility and ease within which systems of PDEs may be numerically approximated using the software developed in this article. Finally, in Section 5 we provide a summary of the work undertaken in this article, as well as outlining potential future developments.

2 Discontinuous Galerkin Finite Element Methods

As a representative PDE example, in this section we outline the DGFEM discretisation for the following system of conservation laws:


where is an open bounded domain in , , with boundary . Here, , , and represent the convective and diffusive fluxes, respectively, which are assumed to be continuously differentiable, and is a given source function. For simplicity of presentation, we assume that (1) may be supplemented with the boundary conditions:


where , and and are two disjoint subsets, with nonempty and relatively open in . We stress that more general boundary conditions can also be considered; for example, in the case when (1) represents the compressible Euler or Navier-Stokes equations, we refer to, for example, [17] and [20], respectively.

For the purposes of discretisation, we rewrite (1) in the following equivalent form:

Here, the matrices


are the homogeneity tensors defined by

, . We write the homogeneity tensor product


such that .

For simplicity of presentation, we now proceed to discretise (1)–(3) based on employing the symmetric interior penalty (SIPG) formulation presented in [20]; however, we stress that other DGFEMs can easily be included within this general setting, cf. below. To this end, we partition into a mesh consisting of non-overlapping open element domains , such that . For each , we denote by the unit outward normal vector to the boundary . We assume that each is an image of a fixed reference element , that is, for all . On the reference element we define spaces of polynomials of degree as follows:

Thereby, with this notation, we define the following DGFEM finite element space

An interior face of is defined as the (non-empty) –dimensional interior of , where and are two adjacent elements of , not necessarily matching. A boundary face of is defined as a (non-empty) –dimensional facet of , , where is a boundary element of , such that . We denote by the union of all interior faces of . Let and be two adjacent elements of  , and an arbitrary point on the interior face . Furthermore, let and  be vector- and matrix-valued functions, respectively, that are smooth inside each element . By , we denote the traces of on taken from within the interior of , respectively. Then, the averages of and at are given by and , respectively. Similarly, the jump of at is given by , where we denote by the unit outward normal vector of , respectively. On , we set , and , where denotes the unit outward normal vector to .

The interior penalty DGFEM discretisation of (1)–(3) is given by: find such that


for all in . The subscript on the operator is used to denote the discrete counterpart of , defined elementwise. Here, denotes the (convective) numerical flux function; this may be chosen to be any two–point monotone Lipschitz function which is both consistent and conservative; typical choices include the (local) Lax–Friedrichs flux, the HLLE flux, the Roe flux, and the Vijayasundaram flux, cf. [30, 39].

The penalisation term arising in the DGFEM scheme (2) is given by

where is a (sufficiently large) positive constant, cf. [15]. Moreover, is defined as , if is in the interior of for two neighbouring elements in the mesh , and , if is in the interior of . Here, for a given (open) bounded set , , we write to denote the –dimensional measure of .

Finally, we define the boundary terms present in the form by

where . Here, the viscous boundary flux and the corresponding homogeneity tensor are defined by

The convective boundary flux is defined by

Finally, the boundary function is given according to the type of boundary condition imposed; in the current setting on and on . For further details regarding the imposition of more general boundary conditions, we refer to [20], for example.

3 Computational Framework for DGFEMs

In this section we present the general computational framework for the automatic generation of DGFEM (semi-) linear forms in a concise and easy to use manner. We begin by outlining the treatment of both convective and viscous components arising in conservation laws. In this setting, the DGFEM formulations can be constructed in a consistent manner. We exploit this by designing utility functions to automatically generate these symbolic DGFEM formulations. We then proceed by proposing a hierarchical framework for computing DGFEM formulations of PDE operators. This hierarchy takes advantage of the DGFEM utility functions, providing a modular framework for a suite of DGFEM formulations for various operators arising in PDE problems of engineering interest.

3.1 Automatic Treatment of Convective Terms

In this section we discuss the automatic symbolic representation of the DGFEM discretisation of the term involving the convective flux function . To this end, we recall that the convective numerical flux function , cf. (2), must be defined on the element interfaces. Defining a callable function F_c which specifies we construct the abstraction of the numerical flux function as shown in Listing 1. The methods interior and exterior will be called to construct the flux function on interior and exterior faces, respectively. The method setup is provided for the inheriting class to initialise any members prior to calls to interior and exterior. This design allows for the flux function to be different on the two types of faces present in the mesh.

class ConvectiveFlux:
    def __init__(self):
    def setup(self):
    def interior(self, F_c, u_p, u_m, n):
    def exterior(self, F_c, u_p, u_m, n):
Listing 1: The abstraction of the numerical flux function .

For simplicity, here we consider two implementations of the ConvectiveFlux class, based on employing the local Lax-Friedrichs and HLLE fluxes, though we stress that other fluxes, such as the Vijayasundaram or Roe flux, for example, may also be employed. Thereby, on the boundary of an element , we define the local-Lax Friedrichs flux and HLLE flux by


respectively. Here,

is the local dissipation parameter, which is selected to be the maximum of the (absolute value) of the eigenvalues of the Jacobi matrix

evaluated on the element face. More precisely, we have that

, where denotes the largest eigenvalue (in modulus) of . Additionally,

where denotes the smallest eigenvalue (in modulus) of .

Given F_c, which specifies , the numerical flux functions and can be automatically generated in order to yield the discretisation of the convective term present in the underlying PDE problem; we note that the constructors of both classes LocalLaxFriedrichs and HLLE require the symbolic representation of the eigenvalues of . As an example, we consider the application of the DGFEM employing the local Lax-Friedrichs flux for the numerical approximation of the linear advection equation


where , , and is some given forcing function. Here, the dissipation parameter can be shown to be . The UFL code required to generate the numerical flux function, together with the corresponding element boundary term arising in the DGFEM scheme for this problem is given in Listing 2. We note that the use of the class HLLE defining the HLLE numerical flux follows in an analogous manner.

def F_c(u): return b*u
H = LocalLaxFriedrichs(lambda u, n: dot(b, n))
H.setup(F_c, u(’+’), u(’-’), n(’+’))
conv_interior = H.interior(F_c, u(’+’), u(’-’), n(’+’))*(v(’+’) - v(’-’))*dS
Listing 2: Example of the automatic calculation and symbolic representation of the local Lax-Friedrichs flux for the linear advection equation shown in (8).

In many cases an analytical expression for the eigenvalues of the Jacobi matrix , , may be difficult to evaluate; thereby, packages such as SymPy [37] may be used to compute them symbolically. To this end, can be assembled using symbolic differentiation; the eigenvalues of this matrix may then be computed symbolically by exploiting the Berkowitz algorithm [7]. In Listing 3 we show an example of this method applied to the convective component of the incompressible Navier-Stokes equations, i.e., .

from sympy import *
dim = 2
u = Matrix([Symbol("u%d" % d, real=True) for d in range(dim)])
n = Matrix([Symbol("n%d" % d, real=True) for d in range(dim)])
F_c = u*u.T
B = zeros(dim, dim)
for d in range(dim):
  B += F_c[:, d].jacobian(u)*n[d]
Listing 3: Automatic symbolic algebra computation of flux Jacobian eigenvalues required by the local Lax-Friedrichs and HLLE fluxes applied to the incompressible Navier- Stokes convective flux component .

3.2 Automatic Treatment of Viscous Terms

In this section we develop utility functions which automatically generate the DGFEM discretisation of second–order PDE operators. To this end, we recall from Section 2 that the viscous component of the underlying PDE is given by


here, we have set the right-hand in (9) to zero, for simplicity, in order to concentrate on the discretisation of the viscous term. The semilinear DGFEM formulation of (9), shown in equation (2), is encapsulated by implementations of the abstract class



is a callable function, which defines the viscous flux function, u and v denote the trial and test functions, respectively, is the DGFEM penalisation coefficient, is the homogeneity tensor, and is the face normal.

class DGFemSIPG(DGFemViscousTerm):
  def interior_residual(self, dInt):
    G = self.G
    F_v, u, v, grad_v = self.F_v, self.U, self.V, self.grad_v_vec
    sig, n = self.sig, self.n
    residual = \
      - inner(tensor_jump(u, n), avg(hyper_tensor_T_product(G, grad_v)))*dInt \
      - inner(ufl_adhere_transpose(avg(F_v(U))), tensor_jump(v, n))*dInt \
      + inner(sig(’+’)*hyper_tensor_product(g_avg(G), tensor_jump(u, n)), tensor_jump(v, n))*dInt
    return residual
  def exterior_residual(self, u_gamma, dExt):
    G = self._make_boundary_G(self.G, u_gamma)
    F_v, u, v, grad_u, grad_v = self.F_v, self.U, self.V, grad(self.U), self.grad_v_vec
    n = self.n
    residual = \
      - inner(dg_outer(u - u_gamma, n), hyper_tensor_T_product(G, grad_v)) * dExt \
      - inner(hyper_tensor_product(G, grad_u), dg_outer(v, n)) * dExt \
      + inner(self.sig*hyper_tensor_product(G, dg_outer(u - u_gamma, n)), dg_outer(v, n)) * dExt
    return residual
Listing 4: The class DGFemViscousTerm provides the abstract interface for interior and exterior residual formulation using symbolic algebra. The class DGFemSIPG uses UFL to automatically formulate the interior and exterior terms of the SIPG formulation of the viscous component of (2).

Recalling the definition of the homogeneity tensor (4) and the homogeneity tensor product (5), is automatically computed using the function homogeneity_tensor(F_v, u) and the function hyper_tensor_product(G, tau) computes the homogeneity tensor product. These functions may then be employed for the generation of the corresponding DGFEM formulation. On the basis of these two functions, we introduce the abstract class DGFemViscousTerm; this offers the following three methods for handling each of the boundary components present in the DGFEM scheme:

  1. DGFemViscousTerm.interior_residual(dS) automatically generates terms associated with the interior boundaries present in (2);

  2. DGFemViscousTerm.exterior_residual(u_gamma, ds_i) generates the terms associated with exterior boundary component ds_i present in (2) with boundary condition ;

  3. Finally, DGFemViscousTerm.neumann_residual(gN, ds_i) generates any terms arising from Neumann boundary conditions with flux specification .

We demonstrate an example of the implementation of the DGFemViscousTerm with the SIPG method in Listing 4; for the sake of brevity, we have removed input and error checking. Further implementations such as the non–symmetric interior penalty (NIPG) and Baumann–Oden schemes are available in the classes DGFemNIPG and DGFemBO, respectively. The DGFEM formulation proposed by Bassi & Rebay [6] is challenging in the symbolic framework due to the requirement of solving elementwise problems for the local lifting operator. Such operations are not easily formulated in the UFL for use with DOLFIN, for example; the implementation of DGFEMs defined based on employing local lifting operators will be considered as part of our programme of future research. For a more detailed discussion of the formulation of lifting operators in DOLFIN, we refer to [34, Chapter 5].

An example of the UFL code required to generate the DGFEM semilinear residual discretisation of the quasi-linear second-order PDE problem


is shown in Listing 5.

def F_v(u, grad_u): return (u + 1)*grad(u)
G = homogeneity_tensor(F_v, u)
vt = DGFemSIPG(F_v, u, v, sig, G, n)
residual = dot((u + 1)*grad(u), grad(v))*dx \
         + vt.interior_residual(dS) \
         + vt.exterior_residual(g_D, ds_D) \
         + vt.neumann_residual(g_N, ds_N) \
         - f*v*dx
Listing 5: Example UFL code for the DGFEM discretisation of the quasi-linear PDE (11) using the DGFemSIPG utility class.

3.3 Automatic Generation of DGFEM Formulations

Even with the utility functions outlined in Sections 3.1 and 3.2, the specification of the DGFEM scheme can be very verbose; this is particularly the case when discretising systems comprising many PDE variables with multiple boundary conditions. In this section we propose a hierarchical scheme for the management of DGFEM formulations of PDE operators of increasing complexity.

3.3.1 Boundary Conditions

Firstly, we outline a simple framework for managing the boundary conditions enforced within a given PDE problem. We define the DGBC (discontinuous Galerkin boundary condition) abstract class from which the classes DGDirichletBC and DGNeumannBC inherit. These implementations simply serve to store the boundary condition and the boundary component over which the condition should be enforced. For example, applying a Dirichlet boundary condition as required by the quasi-linear PDE problem stated in (11) simply requires instantiation of DGDirichletBC(ds\_D, g\_D). Similarly, for the imposition of the Neumann boundary condition we construct DGNeumannBC(ds\_N, g\_N); note that Robin conditions can be implemented in an analogous manner. By generating a list of boundary conditions in this manner, we may easily formulate their imposition within the DGFEM scheme. As an example of a series of Dirichlet boundary conditions being automatically generated using DGFemViscousTerm, we refer to Listing 6.

bcs = [DGDirichletBC(bc_1, ds_1), DGDirichletBC(bc_2, ds_2), ...]
vt = DGFemViscousTerm(F_v, u, v, delta)
ext = sum(vt.exterior_residual(bc.get_function(), bc.get_boundary()) for bc in bcs)
Listing 6: Example of the automatic generation of the exterior residual terms of a given DGFemViscousTerm for a list of Dirichlet boundary conditions.

3.3.2 Abstract DGFEM Formulation

class HyperbolicOperator(DGFemFormulation):
    def __init__(self, mesh, V, bcs,
                 F_c=lambda u: u,
                 H=LocalLaxFriedrichs(lambda u, n: inner(u, n))):
        DGFemFormulation.__init__(self, mesh, V, bcs)
        self.F_c = F_c
        self.H = H
    def generate_fem_formulation(self, u, v, dx=None, dS=None):
        if dx is None:
            dx = Measure(’dx’, domain=self.mesh)
        if dS is None:
            dS = Measure(’dS’, domain=self.mesh)
        n = FacetNormal(self.mesh)
        residual = -inner(self.F_c(u), grad(v))*dx
        self.H.setup(self.F_c, u(’+’), u(’-’), n(’+’))
        residual += inner(self.H.interior(self.F_c, u(’+’), u(’-’), n(’+’)), (v(’+’) - v(’-’)))*dS
        for bc in self.dirichlet_bcs:
            gD = bc.get_function()
            dSD = bc.get_boundary()
            self.H.setup(self.F_c, u, gD, n)
            residual += inner(self.H.exterior(self.F_c, u, gD, n), v)*dSD
        for bc in self.neumann_bcs:
            dSN = bc.get_boundary()
            residual += inner(dot(self.F_c(u), n), v)*dSN
        return residual
Listing 7: The HyperbolicOperator class.

We encapsulate the abstraction of a DGFEM scheme in the class DGFemFormulation which prescribes one abstract method generate\_fem\_formulation. The DGFemFormulation constructor requires the mesh , the function space and the list of boundary conditions. This class will serve as the base class for the DGFEM formulation of all derived PDE operators. In this work we describe two direct children of DGFemFormulation: HyperbolicOperator and EllipticOperator. We stress that the flexibility of this design permits DGFEM formulations of other PDE operators, cf. Sections 4.6, 4.7, and 4.8 below.

The class HyperbolicOperator inherits DGFemFormulation and its purpose is to generate DGFEM formulations of the PDE operator . The class implementation is shown in Listing 7; here, generate\_fem\_formulation is overridden to construct the DGFEM formulation of the provided convective flux operator on the interior and exterior faces, as well as on the elements in the mesh. The numerical flux function provided in the constructor H, which implements ConvectiveFlux, is used for the DGFEM formulation. To highlight to modularity of this design, consider an extension of the HyperbolicOperator for the time-dependent Burgers’ equation


Setting , we may recast (12) in the following equivalent form


the implementation of the class SpacetimeBurgersOperator (using LocalLaxFriedrichs by default) is depicted in Listing 8. Here, we note that the derived class must simply specify the form of the convective flux and the flux function ; indeed, once these are defined, the automatic generation of the DGFEM formulation is then handled by the parent DGFemFormulation.

class SpacetimeBurgersOperator(HyperbolicOperator):
    def __init__(self, mesh, V, bcs, flux=None):
        def F_c(u):
            return as_vector((u**2/2, u))
        if flux is None:
            flux = LocalLaxFriedrichs(lambda u, n: u*n[0] + n[1])
        HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, flux)
Listing 8: The SpacetimeBurgersOperator class.

Secondly, we introduce the class EllipticOperator, which inherits DGFemFor\-mulation, and requires the specification of at instantiation. The overridden method generate\_fem\_formulation is then written to automatically generate the interior and exterior boundary, as well as the element, integration terms, implementing all of the concepts of the utility functions for elliptic operators in DGFemViscousTerm. The outline of the EllipticOperator class, which generates the SIPG formulation by default, is given in Listing 9. To highlight the modularity of EllipticOperator in Listing 10 we show the implementation of the class PoissonOperator which specifies , where is the diffusion coefficient. An example of the automatic generation of the DGFEM formulation of the quasi-linear elliptic PDE problem (11) is given in Listing 11.

class EllipticOperator(DGFemFormulation):
    def __init__(self, mesh, fspace, bcs, F_v, C_IP=10.0):
        DGFemFormulation.__init__(self, mesh, fspace, bcs)
        self.F_v = F_v
        self.C_IP = C_IP
    def generate_fem_formulation(self, u, v, dx=None, dS=None, vt=None):
        if dx is None:
            dx = Measure(’dx’, domain=self.mesh)
        if dS is None:
            dS = Measure(’dS’, domain=self.mesh)
        h = CellVolume(self.mesh)/FacetArea(self.mesh)
        n = FacetNormal(self.mesh)
        sigma = self.C_IP*Constant(max(self.fspace.ufl_element().degree()**2, 1))/h
        G = homogeneity_tensor(self.F_v, u)
        if vt is None:
            vt = DGFemSIPG(self.F_v, u, v, sigma, G, n)
        if inspect.isclass(vt):
            vt = vt(self.F_v, u, v, sigma, G, n)
        assert(isinstance(vt, DGFemViscousTerm))
        residual = inner(self.F_v(u, grad(u)), grad(v))*dx
        residual += vt.interior_residual(dS)
        for dbc in self.dirichlet_bcs:
            residual += vt.exterior_residual(dbc.get_function(), dbc.get_boundary())
        for dbc in self.neumann_bcs:
            residual += vt.neumann_residual(dbc.get_function(), dbc.get_boundary())
        return residual
Listing 9: The EllipticOperator class.

3.3.3 Hierarchy of PDE Operators

A natural extension of this framework is a hierarchy of automatically generated DGFEM operators. As shown in Listing 8, the SpacetimeBurgersOperator inherits from HyperbolicOperator. Consider now, for example, the compressible Navier Stokes equations. Here the inheritance chain may begin with the HyperbolicOperator, from which a CompressibleEulerOperator would inherit, and in turn a CompressibleNavierStokesOperator would inherit Compress\-ibleEulerOperator and additionally EllipticOperator for the viscosity terms. Further implementations of each member of this class hierarchy need not only be undertaken by inheritance. Operators deriving from models of large physical systems may more appropriately aggregate sub-operators as necessary. The derived classes must simply manage the function spaces and boundary conditions amongst the aggregated DGFemFormulation members. This framework significantly reduces the code required for subsequent development of DGFEM formulations of PDE operators of increasing complexity. By ensuring that each layer of the hierarchy is correctly verified and fully tested means that DGFEM formulations may be debugged in a very straightforward manner.

class PoissonOperator(EllipticOperator):
    def __init__(self, mesh, fspace, bcs, kappa=1):
        def F_v(u, grad_u):
            return kappa*grad_u
        EllipticOperator.__init__(self, mesh, fspace, bcs, F_v)
Listing 10: The PoissonOperator class need only inherit the EllipticOperator and define its own viscous flux, F\_v.
po = PoissonOperator(mesh, V, DGDirichletBC(ds, gD), kappa=u+1)
residual = po.generate_fem_formulation(u, v) - f*v*dx
Listing 11: Example of implementing the PoissonOperator utility class.

4 Examples

In this section we present a series of examples of increasing complexity to highlight the ease in which each PDE problem may be discretised using a DGFEM formulation, based on employing the class hierarchy proposed and implemented within this article. We stress that the verbosity and complexity of specifying each DGFEM discretisation is vastly reduced within this modular framework. To this end, we consider the discretisation of a simple scalar advection-diffusion equation, the compressible Euler equations, the compressible Navier-Stokes equations posed in both conserved and entropy variables, the indefinite Maxwell problem, and a hyperelasticity problem. In each case, for simplicity of presentation, we employ the SIPG DGFEM discretisation of the second-order PDE operator and a local Lax-Friedrichs flux for the numerical approximation of the convective terms. For brevity, in some of the examples given below only code snippets will be shown, however, full versions of the corresponding python scripts are available from

4.1 Example 1: Advection-diffusion problem

In this first example, we highlight the use of the classes HyperbolicOperator and EllipticOperator given in Listings 79, respectively. To this end, given , with boundary , consider the problem: find such that


where denotes the diffusion coefficient and . Thereby, setting and , the PDE problem (14), (15) can be written in the general form (1)–(3), where , , and . Moreover, it can easily be shown that the dissipation parameter arising in the local Lax-Friedrichs flux, is given by , while the homogeneity tensor , cf. (4), has entries , .

The specification of the numerical fluxes and , when , along with used with the local Lax-Friedrichs flux function within UFL syntax is provided in the code listing presented in Listing 12. Exploiting the software concepts outlined in Section 3 for the automatic computation of DGFEM formulations, the UFL code to construct and solve the resulting DGFEM approximation of  (14), (15) is presented in Listing 13, where we have specified that , , and ; thereby, the analytical solution to (14), (15) is given by . We note that on the final line of the code listing given in Listing 13, the call to DOLFIN’s solve function automatically computes the Gâteaux derivative of the DGFEM semilinear form , cf.  (2), employed for the numerical approximation of (14), (15) and utilises a Newton solver to evaluate the DGFEM solution; for further details, we refer to [2]. Finally, in Figure 1 we show the asymptotic behaviour of the underlying DGFEM on a sequence of uniformly refined triangular meshes with polynomial orders ; as we expect, we observe that the and norms of the error tend to zero at the respective rates of and as the mesh size tends to zero. The complete code employed for this numerical example is provided in the file

    def F_c(u):
        return b*u**2
    H = LocalLaxFriedrichs(lambda u, n: 2*u*dot(b, n))
    def F_v(u, grad_u):
        return (u + 1)*grad_u
Listing 12: Example 1: UFL representation of and of (14).
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, ’DG’, 1)
u, v = Function(V), TestFunction(V)
gD = Expression(’exp(x[0] - x[1])’, element=V.ufl_element())
f = Expression(’-4*exp(2*(x[0] - x[1])) - 2*exp(x[0] - x[1])’,
b = Constant((1, 1))
ho = HyperbolicOperator(mesh, V, DGDirichletBC(ds, gD), F_c, H)
eo = EllipticOperator(mesh, V, DGDirichletBC(ds, gD), F_v)
residual = ho.generate_fem_formulation(u, v) \
           + eo.generate_fem_formulation(u, v) \
           - f*v*dx
solve(residual == 0, u)
Listing 13: Example 1: Automatic DGFEM formulation for the numerical approximation of (14), (15), using the definitions of the fluxes in Listing 12.



Figure 1: Example 1: Convergence of the DGFEM with –refinement: (a) ; (b) .

4.2 Example 2: Compressible Euler Equations

In this second example, we consider the DGFEM discretisation of the compressible Euler equations: find such that


where the convective flux is given by


Here, , , , and denote the density, velocity vector, pressure, and specific total energy, respectively. The equation of state of an ideal gas is given by

where is the ratio of specific heat capacities at constant pressure, , and constant volume, ; for dry air, . Finally, denotes the identity matrix.

def F_c(U):
  rho, u1, u2, E = U[0], U[1]/U[0], U[2]/U[0], U[3]/U[0]
  p = (gamma - 1.0)*rho*(E - 0.5*(u1**2 + u2**2))
  H = E + p/rho
  return as_matrix([[rho*u1,        rho*u2       ],
                    [rho*u1**2 + p, rho*u1*u2    ],
                    [rho*u1*u2,     rho*u2**2 + p],
                    [rho*H*u1,      rho*H*u2     ]])
def alpha(U, n):
    rho, u1, u2, E = U[0], U[1]/U[0], U[2]/U[0], U[3]/U[0]
    p = (gamma - 1.0)*rho*(E - 0.5*(u1**2 + u2**2))
    u = as_vector([u1, u2])
    c = sqrt(gamma*p/rho)
    lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]
    return lambdas
Listing 14: Example 2: Specification of the convective flux and dissipation parameter for the two-dimensional compressible Euler equations in UFL.

The automatic construction of the DGFEM formulation for the numerical approximation of (16) is encapsulated within the class CompressibleEulerOperator. We note that this class simply inherits the HyperbolicOperator class, with the specification of the Euler flux and dissipation parameter required for the definition of the Lax Friedrichs flux, cf. Listing 14 for the case when . We note that in this setting , where is the speed of sound.

Figure 2: Example 2: Computational domain for Ringleb’s flow.

As an illustration of the exploitation of CompressibleEulerOperator, we consider the DGFEM approximation of Ringleb’s flow problem for which an analytical solution may be obtained using the hodograph method, cf. [9]. The DGFEM discretisation of this test case has also been studied in [18]. This problem represents a transonic flow in a curved channel domain, where the flow is mainly subsonic, with a small supersonic region near the right-hand-side wall; cf. Figure 2. Here, inflow/outflow boundary conditions are imposed on the top and bottom boundaries of , while a solid wall condition is specified on the left- and right-hand side boundaries. Following [21] the latter boundary conditions are enforced based on employing a symmetry technique through the specification of the boundary function . More precisely, we treat the walls as part of the Dirichlet boundary, whereby we set

thereby, here and .



Figure 3: Example 2: Convergence of the DGFEM with –refinement for Ringleb’s flow: (a) ; (b) .
gD = RinglebAnalyticalSoln(element=V.ufl_element(), domain=mesh)
u_vec = project(gD, V)
n = FacetNormal(mesh)
slip_proj = as_matrix(((1,            0,            0, 0),
                       (0,  1-2*n[0]**2, -2*n[0]*n[1], 0),
                       (0, -2*n[0]*n[1],  1-2*n[1]**2, 0),
                       (0,            0,            0, 1)))
slip_bc = slip_proj * u_vec
bcs = [DGDirichletBC(ds(0), gD), DGDirichletBC(ds(1), slip_bc)]
ceo = CompressibleEulerOperator(mesh, V, bcs)
residual = ceo.generate_fem_formulation(u_vec, v_vec)
solve(residual == 0, u_vec)
Listing 15: Example 2: Code snippet for the automatic generation of the DGFEM formulation for the numerical approximation of Ringleb’s flow.

Functions employed for the specification of the analytical solution of Ringleb’s flow, together with routines for the construction of the computational mesh are provided within the file The curved boundaries on the walls of the domain must be treated in a careful manner to ensure optimal convergence of the underlying DGFEM discretisation. Indeed, our computational experience suggests that a curved polynomial description of the boundary of order is necessary to ensure that the and norms of the error tend to zero at the optimal rates of and as the mesh size tends to zero, for fixed , cf. Figure 3. The complete python code for this numerical example is provided; a snippet of this code is depicted in Listing 15 in order to highlight the simplicity of the specification of the DGFEM for this example.

4.3 Example 3a: Compressible Navier-Stokes Equations

In this example we consider the DGFEM discretisation of the compressible Navier-Stokes equations: find such that


where is defined as in (17) and the viscous flux is given by


where denotes the temperature and is the thermal conductivity coefficient. The stress tensor is defined by

where is the dynamic viscosity coefficient. Finally, we note that , where is the Prandtl number.

def F_v(U, grad_U):
    rho, rhou, rhoE = conserved_variables(U)
    u = rhou/rho
    grad_rho = grad_U[0, :]
    grad_rhou = as_matrix([[grad_U[j,:] for j in range(1, dim + 1)]])[0]
    grad_rhoE = grad_U[-1,:]
    # Quotient rule to find grad(u) and grad(E)
    grad_u = as_matrix([[(grad_rhou[j,:]*rho - rhou[j]*grad_rho)/rho**2 for j in range(dim)]])[0]
    grad_E = (grad_rhoE*rho - rhoE*grad_rho)/rho**2
    tau = mu*(grad_u + grad_u.T - 2.0/3.0*(tr(grad_u))*Identity(dim))
    K_grad_T = mu*gamma/Pr*(grad_E - dot(u, grad_u))
    r = as_matrix([[0.0,                            0.0                             ],
                   [tau[0,0],                       tau[0,1]                        ],
                   [tau[1,0],                       tau[1,1]                        ],
                   [dot(tau[0,:], u) + K_grad_T[0], (dot(tau[1,:], u)) + K_grad_T[1]]])
    return r
Listing 16: Example 3a: Specification of the viscous flux for the two-dimensional compressible Navier-Stokes equations in UFL.

Given the UFL code in Listing 14, together with the specification of the viscous flux in Listing 16, in the case when , we have implemented the class CompressibleNavierStokesOperator class, which inherits both the CompressibleEuler\-Operator and EllipticOperator classes. Indeed, the CompressibleEulerOperator component is treated in an identical manner as in the previous section, while the EllipticOperator component is employed with the UFL specification of the viscous flux. On the basis of these classes, the homogeneity tensor and the resulting symbolic DGFEM formulation can then be automatically generated. Moreover, the Gâteaux derivative of the DGFEM formulation is automatically computed by UFL for use within the Newton solver managed in DOLFIN by invoking the call to solve. As a simple test, we consider the example outlined in [20]; namely, we set and select so that the analytical solution to (4.3) is given by


where . Furthermore, we set , and . The snippet of UFL code required to solve this problem is given in Listing 17; the complete code is provided in The orders of convergence of and are reported in Figure 4 as the mesh is uniformly refined for ; as in [20] we observe that and