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++, DUNE , deal.II  and OpenFOAM , 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. .
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 , and Manycore Form Compiler , 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  we originally developed easy-to-use DGFEM utility functions as part of our inhouse software package AptoPy , for application with our own FEM package AptoFEM . AptoPy is written in Python and exploits the open source symbolic algebra package SymPy . We stress that these choices are entirely user-dependent; indeed, in the past we have employed the symbolic algebra packages REDUCE , 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 https://bitbucket.org/nate-sime/dolfin_dg. 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. .
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,  and , 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 ; 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 .
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. . 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 , 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 which specifies we construct the abstraction of the numerical flux function as shown in Listing 1. The methods and will be called to construct the flux function on interior and exterior faces, respectively. The method is provided for the inheriting class to initialise any members prior to calls to and . This design allows for the flux function to be different on the two types of faces present in the mesh.
For simplicity, here we consider two implementations of the of an element , we define the local-Lax Friedrichs flux and HLLE flux byclass, 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
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 , 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 and 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, which specifies
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 defining the HLLE numerical flux follows in an analogous manner.
In many cases an analytical expression for the eigenvalues of the Jacobi matrix , , may be difficult to evaluate; thereby, packages such as SymPy  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 . In Listing 3 we show an example of this method applied to the convective component of the incompressible Navier-Stokes equations, i.e., .
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, is the DGFEM penalisation coefficient, is the homogeneity tensor, and is the face normal.and denote the trial and test functions, respectively,
Recalling the definition of the homogeneity tensor (4) and the homogeneity tensor product (5), is automatically computed using the function and the function 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 ; this offers the following three methods for handling each of the boundary components present in the DGFEM scheme:
present in (2);automatically generates terms associated with the interior boundaries
2) with boundary condition ;generates the terms associated with exterior boundary component present in (
Finally, .generates any terms arising from Neumann boundary conditions with flux specification
We demonstrate an example of the implementation of the 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 and , respectively. The DGFEM formulation proposed by Bassi & Rebay  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].with the SIPG method in Listing
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.
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 11) simply requires instantiation of . Similarly, for the imposition of the Neumann boundary condition we construct ; 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 , we refer to Listing 6.(discontinuous Galerkin boundary condition) abstract class from which the classes and 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 (
3.3.2 Abstract DGFEM Formulation
We encapsulate the abstraction of a DGFEM scheme in the class , 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 : and . 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.which prescribes one abstract method The constructor requires the mesh
The class . The class implementation is shown in Listing 7; here, 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 , which implements , is used for the DGFEM formulation. To highlight to modularity of this design, consider an extension of the for the time-dependent Burgers’ equationinherits and its purpose is to generate DGFEM formulations of the PDE operator
Setting , we may recast (12) in the following equivalent form
the implementation of the class 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 .(using by default) is depicted in Listing
Secondly, we introduce the class at instantiation. The overridden method 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 . The outline of the class, which generates the SIPG formulation by default, is given in Listing 9. To highlight the modularity of in Listing 10 we show the implementation of the class 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., which inherits , and requires the specification of
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 inherits from . Consider now, for example, the compressible Navier Stokes equations. Here the inheritance chain may begin with the , from which a would inherit, and in turn a would inherit and additionally 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 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.
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 https://bitbucket.org/nate-sime/dolfin_dg.
4.1 Example 1: Advection-diffusion problem
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 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 . 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 advection_diffusion.py.
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.
The automatic construction of the DGFEM formulation for the numerical approximation of (16) is encapsulated within the class . We note that this class simply inherits the 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.
As an illustration of the exploitation of 9]. The DGFEM discretisation of this test case has also been studied in . 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  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, we consider the DGFEM approximation of Ringleb’s flow problem for which an analytical solution may be obtained using the hodograph method, cf. [
thereby, here and .
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 ringleb.py. 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 ringleb_example.py; 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.
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 class, which inherits both the and classes. Indeed, the component is treated in an identical manner as in the previous section, while the 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 . As a simple test, we consider the example outlined in ; 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 compressible_navierstokes_square.py. The orders of convergence of and are reported in Figure 4 as the mesh is uniformly refined for ; as in  we observe that and