Extendible and Efficient Python Framework for Solving Evolution Equations with Stabilized Discontinuous Galerkin Method

This paper discusses a Python interface for the recently published DUNE-FEM-DG module which provides highly efficient implementations of the Discontinuous Galerkin (DG) method for solving a wide range of non linear partial differential equations (PDE). Although the C++ interfaces of DUNE-FEM-DG are highly flexible and customizable, a solid knowledge of C++ is necessary to make use of this powerful tool. With this work easier user interfaces based on Python and the Unified Form Language are provided to open DUNE-FEM-DG for a broader audience. The Python interfaces are demonstrated for both parabolic and first order hyperbolic PDEs.

DUNE-FEM-DG has been used in several applications (see [13] for a detailed list), most notably a comparison with the production code of the German Weather Service COSMO has been carried out for test cases for atmospheric flow ( [7,45]). The focus of the implementation is on Runge-Kutta DG methods using mainly a matrix-free approach to handle implicit time discretizations which is a method especially used for convection dominated problems.
A strength of the DUNE-FEM-DG module is the general application area, i.e. convection dominated as well as diffusion dominated problems, for 1d, 2d, and 3d models, including parallelization and local grid adaptivity. The model includes the implementation of several discretizations for the second order terms such as Interior Penalty (and variants), Compact Discontinuous Galerkin 1 and 2, and Bassi-Rebay 1 and 2 as well as Local Discontinuous Galerkin. For the first order terms different numerical flux functions can be used and various limiter based stabilizations are available. For the time discretization a method of lines approach is adopted and a number of implicit, explicit, and IMEX schemes are available. Overall the module thus offers not only a strong base for building state of the art simulation tools but it also allows for the development of new methods and comparative studies of different method. Parallelization and adaptivity, including hp adaptivity, can be used seamlessly with all methods.
A shortcoming so far has been the template heavy and relatively complicated C++ user interfaces, leading to a steep learning curve for implementing new models and applications or coupling of such. Recent development (e.g. [18]) has therefore been focused on adding a Python layer on top of DUNE-FEM-DG allowing user friendly model description based on the Unified Form Language (UFL) [1] and code generation. While quite a few applications using UFL and also consider DG approximation, e.g., the FEniCS package [41], only a few examples exits that consider convection dominated evolution equations. For example, in [31] UFL is used to describe weak forms for the compressible Euler and Navier-Stokes equations but only in the stationary setting. Another example is [47] where shallow water applications are considered. A package that also provides Python bindings (not using UFL) but focuses on hyperbolic problems is [42], where higher order Finite-Volume schemes are the method of choice. To our best knowledge this is the first package which combines high level scripting support with efficient stabilized Discontinuous Galerkin methods for solving the full range of pure hyperbolic systems, through advection dominated problems, to diffusion equations.
In DUNE low level Python bindings were introduced for the DUNE grid interface in [22] and a detailed tutorial providing high level access to DUNE-FEM is also available [15]. In this work, we describe how these bindings can now be used. The new interface, together with the efficient and flexible C++ implementation of the DG methods available in DUNE-FEM-DG, make the rapid prototyping of new methods straightforward and provides a flexible framework to solve very complex coupled nonlinear evolution equations. Other works in this direction

Motivation and aim of this paper
This paper describes a collaborative effort to establish a test and research environment for DG methods based on DUNE and DUNE-FEM for advection-diffusion problems ranging from advection dominated to diffusion only problems. The aim here is to provide easy access to a comprehensive collection of existing methods and techniques to serve as a starting point for new developments and comparisons without having to reinvent the wheel. This is combined with the easy access to advanced features like the availability of different grid element types, not limited to dimensions less than or equal to three, hp adaptation, and parallel computation with both distributed (e.g. MPI) and shared memory (e.g. OpenMP) parallelization, with excellent scalability [36] even for large core counts and very good oncore efficiency in terms of floating point operations per second [13]. Python has become a widespread programming language in the scientific computing community including the use in industry we see great potential in improving the scientific development of DG methods by providing access to state of the art research tools which can be used through the Python scripting language, backed up by an efficient C++ backend.
The focus of this paper is to describe the different parts of the framework and how they can be modified by the user in order to tailor the code to their needs and research interests. In this paper we will mainly look at advection dominated evolution equations because most diffusion dominated problems can be often completely described within a variational framework for which the domain specific language UFL is very well suited and the standard code generation features available in DUNE-FEM are thus sufficient. For advection dominated problems additional ingredients are often required, e.g., for stabilization, that do not fit the variational framework. This aspect of the algorithm will be a central part of this paper. The development and improvement of the DG method in general or the high performance computing aspect of the underlying C++ framework, which has been investigated previously, are outside the scope of this work.
The paper is organized as follows. In Section 1 we briefly recall the main building blocks of the DG discretization of advection diffusion problems. In Section 2 we introduce the newly developed Python based model interface. In Section 3 we investigate the performance impact of using Python scripting and conclude with discussing the extensibility of the approach in Section 4. Additional code examples are available in Appendix A.

Governing Equations and Discretization
We consider a general class of time dependent nonlinear advection-diffusion-reaction problems for a vector valued function U : (0, T ) × Ω → R r with r ∈ N + components of the form in Ω ⊂ R d , d = 1, 2, 3. Suitable initial and boundary conditions have to be added. F c describes the convective flux, F v the viscous flux, S i a stiff source term and S e a non-stiff source term. Note that all the coefficients in the partial differential equation are allowed to depend explicitly on the spatial variable x and on time t but to simplify the presentation we suppress this dependency in our notation. Also note that any one of these terms is also allowed to be zero. For the discretization we use a method of lines approach based on first discretizing the differential operator in space using a DG approximation and then solving the resulting system of ODEs using a time stepping scheme.

Spatial Discretization
Given a tessellation T h of the domain Ω with ∪ K∈T h K = Ω we introduce a piecewise is a space containing all polynomials up to degree k.
Furthermore, we denote with Γ i the set of all intersections between two elements of the grid T h and accordingly with Γ the set of all intersections, also with the boundary of the domain Ω .
We seek U h ∈ V h be discretizing the spatial operator L (U) in (1) with either Dirichlet,  Neumann, or Robin type boundary conditions by defining for all test functions with the element integrals with S(U h ) = S i (U h ) + S e (U h ) and the surface integrals (by introducing appropriate numerical fluxes F c , F v for the convection and diffusion terms, respectively) with {{U}} e , [[U]] e denoting the classic average and jump of U over e, respectively. The convective numerical flux F c can be any appropriate numerical flux known for standard finite volume methods e.g. F c could be simply the local Lax-Friedrichs flux function (also known as Rusanov flux) where λ e is an estimate of the maximum wave speed on intersection e. One could also choose a more problem tailored flux (i.e. approximate Riemann solvers). Different options are implemented in DUNE-FEM-DG (cf. [13]). A wide range of diffusion fluxes F v can be found in the literature and many of these fluxes are available in DUNE-FEM-DG, for example, Interior Penalty and variants, Local DG, Compact DG 1 and 2 as well as Bassi-Rebay 1 and 2 (cf. [8,13]).

Temporal discretization
To solve the time dependent problem (1) we use a method of lines approach in which the DG method described above is first used to discretize the spatial operator and then a solver for ordinary differential equations is used for the time discretization, see [13]. After spatial discretization, the discrete solution U h (t) ∈ V h has the form U h (t, x) = ∑ i U i (t)ϕ i (x). We get a system of ODEs for the coefficients of U(t) which reads with f (U(t)) = M −1 L h (U h (t)), M being the mass matrix which is in our case is block diagonal or even the identity, depending on the choice of basis functions. U(0) is given by the projection of U 0 onto V h . A range of different ODE solvers are available most based around Strong Stability Preserving Runge-Kutta methods (SSP-RK) (for details see [13]). In Section 2 we show how to easily implement other Runge-Kutta methods into the code. The results and implementation techniques presented in this paper can be applied to explicit, implicit, or semi-implicit methods and mostly the non-linear system arising from any implicit treatment is solved using a matrix-free Newton-Krylov method.
When using semi-implicit time stepping, the operator L is split such that L (U) = L e (U) where L e is treated explicitly and is L i is treated implicitly. Of course other splittings are possible but currently not available through the Python bindings.

Stabilization
The RK-DG method is stable when applied to linear problems such as linear hyperbolic systems; however for nonlinear problems spurious oscillations occur near strong shocks or steep gradients. In this case the RK-DG method requires some extra stabilization. In fact, it is well known that only the first order scheme (k = 0) produces a monotonic structure in the shock region. Many approaches have been suggested to make this property available in higher order schemes, without introducing the amount of numerical viscosity, which is such a characteristic feature of first order schemes. Several approaches exist, among those are slope limiters [12,34,16,40,25] and for a comprehensive literature list we refer to [46]. Another popular stabilization technique is the artificial diffusion (viscosity) approach [27,35,44,29,23] and others. Further techniques exists, such as a posteriori techniques to stabilize the DG method [20] or order reduction [24]. The construction of the general stabilization approach implemented in DUNE-FEM-DG is described in detail in [16] and follows the slope limiting approach combined with a troubled cell indicator suggested in [40]. However, since multiple discretization techniques for the diffusion operator are available within DUNE-FEM-DG, the implementation of an artificial diffusion approach would be straight forward.
In the following we briefly recall the main steps since these are necessary to understand the code design decision later on.
A stabilized discrete operator is constructed by concatenation of the DG operator L h from (2) and a stabilization operator Π h , leading to a modified discrete spatial operator The stabilizations considered in this work can be computed element wise based only on data from neighboring elements thus not increasing the stencil of the operator. Given a DG function U h we call U * h = Π h (U h ) the stabilized DG function and we call U E = U h | E the restriction of a function U h on element E and denote withŪ E its average. Furthermore, we call e the intersection between two elements E, K for K ∈ N E , with N E being the set of neighbors of E and I E the set of intersections of E with it's neighbors.
The stabilized solution should fulfil the following requirements: 6. Maximum-minimum principle and monotonicity, i.e. in regions where the solution is not smooth the function U * E should only take values between min K∈N EŪ K and max K∈N EŪ K . These requirements are built into the stabilization operator in different ways. For example, requirement 2 and 3 form the so called troubled cell indicator that triggers whether a stabilized solution has to be computed. Requirement 5 simply limits the choice of available reconstruction methods, e.g. higher order finite volume reconstructions could not be used because of the potentially larger stencil that would be needed which is not desired here.
The stabilization is then defined in two steps 1. We define the set of troubled cells by with J E being a smoothness indicator and T OL a threshold that can be influenced by the user (see Figure 9 in Section 2.3). As default we use T OL = 1 and a smoothness indicator based on the superconvergence result from [40]: this indicator accumulates the integral of the jump of some scalar quantity derived from U h denoted with φ (U E ,U K ) over all inflow boundaries of E (defined as the intersection where some given velocity v and the normal n e have opposite signs), i.e.
α d (k) = 2 125 d5 k denotes a scaling factor, h E is the elements diameter and |e| the area of the intersection between the two elements E and K. Please note the slight derivation from the notation in [16] by denoting the smoothness indicator in (9) with J E instead of S E like in [16] In this case the DG solution on E needs to be altered until E is no longer a troubled cell. In this work we guarantee this by restricting ourselves to the reconstruction of limited linear functions based on the average values of U h on E and it's neighbors. Two reconstruction methods are implemented, one described in [16] and extended to general polyhedral cells in [37] and an optimization based strategy described in [9] based in ideas from [43]. Both strategies correspond to 2nd order MUSCL type finite volume schemes when used with piecewise constant basis functions.

Summary of building blocks and limitations
-Grid structure: with discontinuous Galerkin methods we can use any form of grid elements, from (axis aligned or general) cubes, simplices, to general polyhedrons). The best choice will depend on the application. In addition different local adaptation strategies can be considered as part of the underlying grid structure, e.g., conforming refinement using red-green or bisection strategies, for example, or simple nonconforming refinement. From a parallelization point of view DUNE-FEM-DG is restricted to DG methods that can be implemented using direct neighboring information only. On the one hand this is beneficial for many core architectures and on the other hand unstructured grids in DUNE only implement this ghost cell approach and arbitrary overlap is only available for Cartesian grids.
-Space: an advantage of discontinuous Galerkin methods is that there is little restriction on the combination of grids and discrete spaces one can use. For approximation purposes these spaces are generally chosen so that they contain the full polynomials of a given order on each element. This order is perhaps the most important property of the spaces. But even after having fixed the space, the actual set of basis functions to use to represent the space, can have a huge impact on the performance of the scheme. Central here is the structure of the local mass matrix on each element of the grid. This has to be inverted in each time step so efficiency here can be crucial. So a possible choice is to use a polynomial space with basis functions orthonormalized over the reference elements in the grid. If the geometric mapping between physical grid elements and these reference element is affine then the resulting local mass matrix is a very simple diagonal matrix. If the mapping is not affine then the mass matrix depends non trivially on the element under consideration and will often be dense. In this case another approach is to use Lagrange type basis functions where the interpolation points coincide with a suitable quadrature for the mass matrix. The case of non affine mapping is especially relevant for cube grids and here tensor product quadrature points can be used to construct the Lagrange type space. A possible approach is to use a tensor product Gaussian rule which results in a diagonal mass matrix for each element in the grid with a trivial dependency on the element geometry. Another typical approach found in the literature is often referred to a spectral discontinuous Galerkin method based on Lobatto-Gauss-Legendre quadrature [39,38]. Although this quadrature is not precise enough to exactly integrate products of the basis functions, it has been shown that diagonalizing (lumping) the mass matrix in this way leads to a well posed scheme. -Numerical fluxes: the discretization of both the diffusion and advective terms, involve boundary integrals where the discrete solution is not uniquely defined. Consequently, numerical fluxes have to be used to evaluate these terms. We consider fluxes here which depend on the traces of the discrete solution on both sides of the interface. A wide range of choices are available. In the package presented here a large number of fluxes for the diffusion are available, but since we are focusing on advection dominated problems we will not discuss this further. The simplest choice for the advective flux is the Local-Lax-Friedrichs or Rusanov flux from equation (5) which requires little additional input from the user and in combination with higher order schemes generally gives good results. -Time evolution: after the spatial discretization, the resulting ordinary differential equation is solved using a Runge-Kutta as already mentioned. We will only briefly discuss this aspect of the method in the following. -Stabilization: as was discussed in some detail in the previous section, our stabilization approach is based on a troubled cell indicator (combining a smoothness indicator with a check on physicallity at quadrature points) and a reconstruction process of the solution in troubled cells.
If not mentioned otherwise, the results presented in the following are obtained using a Cartesian grid with a forth order polynomial space spanned by orthonormal basis function. We use the local Lax Friedrichs flux (5) for the advection term and the CDG2 method [8] for the problems with diffusion. The troubled cell indicator is based on a smoothness indicator described above. The reconstruction is computed based on a linear programming problem as described in [9]. Finally, we use a standard third order SSP Runge-Kutta method for the time evolution [28].
In the following we will describe how the different building blocks making up the DG method can be provided by the user, starting with some very simple problems, requiring little input from the user and slowing building up the complexity. We are not able to provide the full code listings in each case but a tutorial with all the test cases presented here is available [17].
2 Customizing the DG method using the Python Interface First we need to choose the grid structure and the computational domain. For simplicity we concentrate here on problems defined on 2D intervals [x l , x r ]×[y b , y t ] and use an axis aligned cube grid capable of non-conforming refinement, with n x , n y elements in x and y direction, respectively. The space consists of piecewise polynomials of degree p = 4 spanned by an orthonormal basis over the reference element [0, 1] 2 . This can be setup with a few lines of Python code: from dune . alugrid import aluCubeGrid as grid from dune . fem . view import a d a p t i v e L e a f G r i d V i e w as view # needed for adaptive simulation from dune . fem . space import dgonb gridView = view ( grid ( Model . domain ) ) gridView . hierarchicalGrid . globalRefine ( 3 ) # refine a possibly coarse initial grid The Model class will be specified in the following. It contains information on the PDE but as seen above also provides a description of the domain, the number of components of the PDE to solve and the initial conditions U 0 defined as a UFL expression. Examples will be provided in the next sections. After setting up the grid, the space, and the discrete solution, we define the DG operator and the default Runge-Kutta time stepper: from dune . femdg import femDGOperator from dune . femdg . rk import femdgStepper operator = femDGOperator ( Model , space , limiter = None ) stepper = femdgStepper ( order =3 , operator = operator ) Again we use the Model class to provide all required information to the DG operator. The constructor of the DG operator takes a number of additional parameters which will described throughout this section as required (see also Appendix A). In our first example we will require no stabilization so the limiter argument is set to None. In the final line of code, we pass the constructed operator to the time stepper together with the desired order (as pointed out we use order 3 for all presented simulations). Finally, a simple loop is used to evolve the solution from the starting time (assumed to be 0) to the final time T (which is again a property of the Model class): The call method on the stepper evolves the solution U h from the current to the next time step and returns a suggested size ∆t for the next time step. In the following we will describe the Model class in detail and show how the above code snippets have to be modified to include additional features like stabilization for nonlinear advection problems or adaptivity. We will not describe each problem in detail, the provided code should give all the required information.

Linear advection problem
In the following we solve a linear advection problem and explain various options for stabilization and local grid adaptivity.

Model description
from dune . ufl import cell from dune . grid import cartesianDomain from ufl import * x = SpatialCoor dinate ( cell ( 2 )  Note that to use the local Lax-Friedrichs flux we need to define the maximum wave speed of the hyperbolic problem. This function is also used to provide an estimate for the time step based on the CFL condition, which is returned by the stepper. The stepper also provides a property deltaT which allows the used to fix a time step to use for the evolution. A result for the described three body problem is presented in Fig. 1 (left) and Fig. 2 (left).

Stabilization
It is well known that the DG method with a suitable numerical flux is stable when applied to linear advection problems like the one studied here. But it does produce localized over and undershoots around steep gradients. While this is not necessarily a problem for linear equations, it can be problematic in some instances, e.g., when negative values are not acceptable. More crucially this behavior leads to instabilities in the case of nonlinear problems. Therefore, stabilization mechanism have to be provided. As mentioned above we use a reconstruction approach combined with a troubled cell indicator. The standard troubled cell indicator is described in Section 1.3 and requires the user to provide some additional information in the Model class, i.e. a description of the jump of the solution φ (U,V ) between of two elements E and K and a suitable velocity v required to compute the smoothness indicator given in equations (9):  to the code. Another alternative is to use the approach detailed in [48,11]  The resulting plots using the scaling limiter and simply limiting in cells where the solution is below zero or above one, are shown in Fig. 3.

Adaptivity
Due to the regions of steep gradients and low regularity in the solution, advection dominated problems benefit a lot from the use of local grid refinement and coarsening. Although this is a more general feature of the DUNE-FEM framework on which the package is based, adaptivity is such an important tool to improve efficiency of schemes for the type of problems discussed here, that we will describe how to use it in some detail. The grid modification requires an indicator computed between time steps which provides the information for each cell whether to keep it as it is, refine it, or coarsen it. DUNE-FEM provides a very simple function dune . fem . markNeighbors ( indicator , refineTol , coarsenTol , minLevel , maxLevel ) Here indicator provides a number η E ≥ 0 on each element E which is used to determine if an element is to be refined or coarsened. An element E is refined, if it's level in the grid hierarchy is less then minLevel and η E >refineTol; it is coarsened if it's level is greater then maxLevel and η E <coarsenTol; otherwise it remains unchanged. In addition if a cell is to be refined all its neighboring cells are refined as well, so that important structures in the solution do not move out of refined regions during a time step. If the initial grid is suitably refined, then one can take maxLevel=0. The choice for minLevel will depend on the problem and has to be chosen carefully since increasing minLevel by one can potentially double the runtime simply due to the reduction in the time step due to the CFL condition. At the moment we do not provide any spatially varying time step control. We usually choose coarsenTol simply as a fraction of refineTol. In our simulation we use coarsenTol=0.2*refineTol. This leaves us needing to describe our choice for indicator and refineTol.
The simulations shown here use a residual type indicator based on a scalar quantity and flux provided in a subclass Indicator of the Model class. We consider the residual of an auxiliary PDE where for example η, F, S could simply be taken from one of the components of the PDE or could be based on an entropy, entropy flux pair. We now define the element residual to be where R vol is a discretized version of the interior residual indicating how accurate the discretized solution satisfies the auxiliary PDE at every interior point of the domain for two timesteps t n and t n+1 , and we have jump indicators This indicator is just one of many possible choices, e.g., simply taking the jump of U h over intersections will often lead to very good results as well. We have had very good experiences with this indicator as shown here see also [14] where we used a similar indicator and [20] where a similar indicator was derived.

Compressible Euler equations
We continue our presentation with a system of evolution equations. We focus on simulations of the compressible Euler (and later Navier-Stokes) equations. We show results for two standard test cases: the interaction of a shock with a low density region and the simulation of a Kelvin-Helmholtz instability with a density jump.

Model description
We start with the methods needed to describe the PDE, the local Lax-Friedrichs flux and time step control, the troubled cell indicator, and the residual indicator:    Now that we have set up the model class, the code presented for the advection problem for evolving the system and adapting the grid can remain unchanged.
Results for both test cases on locally adapted grids are shown in Fig. 6 and the left column of Fig. 8. It turns out that in the Kelvin-Helmholtz case the stabilization is almost completely determined by the physicality check since the smoothness indicator shown here is based on the pressure which in this case is continuous over the discontinuity. To increase the stabilization of the method it is either possible to reduce the tolerance in the troubled cell indicator (results shown on right of Fig. 9) or to use a different smoothness indicator all together, which is discussed in Section 2.3.

Adding source terms
So far we have in fact simulated the interaction of a shock wave with a column of low density gas and not in fact a bubble. To do the later would require to either extend the problem to 3D (discussed in the next section) or to simulate the problem using cylindrical coordinates. This requires adding a geometric source term to the right hand side of the Euler equations. As pointed out in the introduction our model can take two types of source term depending upon the desired treatment in the time stepping scheme. Here we want to treat the source explicitly so need to add an S_e method to the Model class:   Results are shown in Fig. 7.

Adding diffusion
Finally, we discuss the steps needed to add a diffusion term. This requires adding an additional method to the Model class. So to solve the compressible Navier-Stokes equations instead of the Euler equations the following method needs to be added (given a viscosity parameter µ):  As an example we repeat the simulation of the Kelvin-Helmholtz instability (see Fig. 8).
Note that with the default setting for the stepper, an IMEX scheme is used where the diffusion is treated implicitly and the advection explicitly with a time step given by the CFL condition.

User defined smoothness indicator
To exchange the smoothness indicator the construction of the operator has to be slightly changed. Assuming that the indicator is defined in a source file modalindicator.hh which defines a C++ class ModalIndicator which takes the C++ type of the discrete function U_h as template argument. This class needs to be derived from a pure virtual base class and override a single method. Then the construction of the operator needs to be changes to As an example we use here a smoothness indicator based on studying the decay properties of the modal expansion of the solution on each cell following the ideas presented in [35]. For the following implementation we assume that we are using a modal basis function The actual computation of the indicator is carried out in the method smoothnessIndicator but is slightly too long to include here directly but is shown in Appendix B. It is important to note that it requires very little knowledge of the DUNE programming environment or even C++ since it relies mainly on the local degrees of freedom vector provided by the argument uEn.

User defined time stepping schemes
The DUNE-FEM package provides a number of standard strong stability preserving Runge-Kutta (SSP RK) solvers including explicit and implicit methods and IMEX schemes method of degree one to four. In the literature there is a wide range of additional suitable RK methods (having low storage or better CFL constants using additional step for example). Furthermore, multistep methods can be used. Also there are a number of other packages providing implementations of timestepping methods. Since the computationally critical part of a DG method of the type described here lies in the computation of the spatial operator, the additional work needed for the timestepper can be carried out on the Python side with little impact. Furthermore, as pointed out in the introduction, it is often desirable to use Python for rapid prototyping and to then reimplement the finished algorithm in C++ after a first testing phase to avoid even the slightest impact on performance. The following code snippet shows how a multistage third order RK method taken from [33] can be easily implemented in Python and used to replace the stepper used so far:  Again the remaining code can remain unchanged. Results with this time stepping scheme are included in some of the comparisons shown in the next section (see Fig. 12 and Fig. 13).

Different grids and spaces
One of the strengths of the DUNE framework on which we are basing the software presented here, is that it can handle many different types of grid structures. Performing a 3D simulation can be as simple as changing the domain attribute in the Model class (Fig. 10).
It is also straightforward to perform the simulations on for example a simplicial grid instead of the cube grid used so far: from dune . alugrid import aluSimplexGrid as grid gridView = view ( grid ( Model . domain ) ) gridView . hierarchicalGrid . globalRefine ( 3 ) # refine a possibly coarse initial grid It is even possible to use a grid consisting of general polygonal elements, by simply importing the correct grid implementation: from dune . polygongrid import polygonGrid as grid A wide range of other grid types are available and a recent overview is given in [5] some examples are shown in Fig. 11. As pointed out in the previous section, the choice of the basis function set used to represent the discrete solution can strongly influence the efficiency of the simulation. So far we have used an orthonormal basis function set for the polynomial space over the reference cube [0, 1] d . If the grid elements are all affine mapping of [0, 1] d (i.e. parallelograms) then this is good choice, since it has the minimal number of degrees of freedom for a desired approximation accuracy. Also the mass matrix will be a very simple diagonal matrix on all elements in this case. These properties always hold for simplicial elements when an orthonormal polynomial basis over the reference simplex is used. As soon as the mapping between the reference element and a given element in the grid becomes non affine, both properties can be lost. In this case it might be necessary, to achieve the right approximation properties, to use a tensor product polynomial space, increasing the number of degrees of freedom per element considerably. Also an orthonormal set of basis function over the reference cube, will still lead to a dense mass matrix. This lead to a significant reduction in the efficiency of the method, if the local mass matrix on each element can not be stored due to memory restrictions. A possible solution to this problem is to use a Lagrange type set of basis functions with interpolation points coinciding with a quadrature rule over the reference cube. The obvious choice for this quadrature is a tensor product Gauss-Legendre rule. Since this quadrature is accurate up to order 2k + 1, it can be used to exactly compute the mass matrix and consequently, the mass matrix ϕ i ϕ j = ∑ q ω q ϕ i (x q )ϕ j (x q ) = ∑ q ω q δ iq δ jq = ω i δ i j is a diagonal matrix where the diagonal element only depend on the integration element of the reference mapping at the quadrature points. In addition, if all the other element integrals needed to evaluate the spatial operator use the same quadrature, the interpolation property of the basis functions can be used to further speedup evaluation. All of this is well known and for example investigated in [39]. In the software framework presented here, it is straightforward to switch to this representation of the discrete space by replacing the construction of the space object by from dune . fem . space import dglagrange space = dglagrange ( gridView , dimRange = Model . dimRange , order =4 , pointType = " gauss " ) If the operator is constructed using this space, the suitable Gaussian quadrature is chosen automatically. Note that this space is only well defined over a grid consisting of cubes.
It is also common practise to use the points of a tensor product Lobatto-Gauss-Legendre (LGL) quadrature (see e.g. [38]). In this case the evaluation of the intersection integrals in the spatial operator can also be implemented more efficiently. To still retain the simple structure of the mass matrix requires to use the same Lobatto-Gauss-Legendre rule to compute ϕ i ϕ j . Since this rule is only accurate up to order 2k − 1 this results in an underintegration of the mass matrix similar to mass lumping; this does not seem to influence the accuracy of the scheme. It is straightforward to use this space: space = dglagrange ( gridView , dimRange = Model . dimRange , order =4 , pointType = " lobatto " ) and again using this space in the construction of the spatial operator will result in the correct LGL quadrature being used. We compare L 2 errors on a sequence of non affine cube grids non affine GL non affine LGL non affine ONB simplex non affine LGL with 4stage RK Fig. 12: Two rarefaction wave problem simulated from t = 0.05 to 0.12 using different representation for the discrete space with polynomial order 4. The simulations were performed on a sequence of grids consisting of non affine cubes (split into 2 simplices for simplicial simulation). Shown are the errors measured in the L 2 norm versus the number of degrees of freedom (left) and the total runtime (right). The right plot also includes results from a simulation with the LGL method and the SSP3 time stepper implemented in Python as discussed previously. The results of this simulation are not included on the left since they are broadly in line with the LGL simulation using the 3 stage RK method. Due to the larger CFL constant this method is more efficient then the 3 stage method of the same order. The simulation on the simplicial grid produces a slightly better error on the same grid but requires twice as many degrees of freedom so that it is less efficient compared to the LGL or GL simulations. Note that the runtime with the ONB basis is significantly larger due to the additional cost of computing the inverse mass matrix on each element of the grid as discussed above. On an affine cube grid the runtime is comparable to the GL scheme on the same grid but this is not shown here.
using different sets of basis functions in Fig. 12, 13, and 14. We show both the error vs. the number of degrees of freedom (left) and vs. the runtime (right).

Reactive advection diffusion problem
We conclude with an example demonstrating the flexibility of the framework to combine different components of the DUNE-FEM package to construct a scheme for a more complex problem. As a simple example we use a chemical reaction type problem with linear advection and diffusion where the velocity field is given by discretizing the solution to an elliptic problem in a continuous Lagrange space. As mentioned this is still a simple problem but can be seen as a template for coupled problems, e.g., transport in porous media setting or where the flow is given by solving incompressible Navier Stokes equations.
Let us first compute the velocity given as the curl of the solution to a scalar elliptic problem: non affine GL non affine LGL non affine ONB simplex non affine LGL with 4stage RK Fig. 13: Sod's Riemann problem simulated from t = 0 to 0.2 using different representation for the discrete space with polynomial order 4. The simulations were performed on a sequence of grids consisting of non affine cubes (split into 2 simplices for the simplicial simulation). Shown are the errors measured in the L 2 norm versus the number of degrees of freedom (left) and the runtime (right). The right plot also includes results from a simulation with the LGL method and the SSP3 time stepper implemented in Python as discussed previously. The results of this simulation are not included on the left since they are broadly in line with the LGL simulation using the 3 stage RK method. Due to the larger CFL constant this method is more efficient then the 3 stage method of the same order. The simulation on the simplicial grid produces a sginificantly better error on the same grid but requires twice as many degrees of freedom. Overall it is more efficient then the LGL or GL simulations using the same three stage RK method. Note that the runtime with the ONB basis is significantly larger due to the additional cost of computing the inverse mass matrix on each element of the grid as discussed above. On an affine cube grid the runtime is comparable to the GL scheme on the same grid but this is not shown here. We use this velocity field to evolve three chemical components reacting through some nonlinear reaction term and include some small linear diffusion: from ufl import * from dune . ufl import DirichletBC from dune . fem . space import lagrange from dune . fem . scheme import galerkin Note that the source term includes both the chemical reaction and a source for the first two components. The third component is generated by the first two interacting.
Due to the diffusion we do not need any stabilization of the form used so far. However, in this case a reasonable assumption is that all components remain positive throughout the simulation, so the physicality check described above is still a useful feature. We can combine this with the scaling limiter already described for the advection problem. To this end we need to add bounds to the model and change to the construction call for the operator: Model . lowerBound = [0 ,0 , 0 ] operator = femDGOperator ( Model , space , limiter = " scaling " ) Note that by default the stepper switches to a IMEX Runge-Kutta scheme if the Model class contains both an advective and a diffusive flux. This behavior can be changed by using the rkType parameter in the constructor call for the stepper. A final remark concerning boundary conditions: here we use simple Dirichlet boundary conditions which are then used for both the diffusive and advective fluxes on the boundary as outside cell value. We saw in a previous example that we can also prescribe the advective flux on the boundary directly. In that example we used Model . boundary [ 3 ] = lambda t ,x ,u , n : noFlowFlux (u , n ) If this type of boundary conditions is used we also need to prescribe the flux for the diffusion term, so for an advection-diffusion problem we pass in a pair of fluxes at the boundary, e.g.,

Efficiency of Python based auto-generated models
While Python is easy to use, its flexibility can lead to some deficiencies when it comes to performance. In DUNE-Python [22] and DUNE-FemPy [15] a just-in-time compilation concept is used to create Python modules based on the static C++ type of every object used, i.e. the models described in the previous section are translated into C++ code based on the UFL descriptions in the various model methods which is then compiled and loaded as Python modules. This way we avoid virtualization of the DUNE interfaces and consequently one would expect very little performance impact as long as calls between Python and C++ are only done for long running methods. To verify this, we compare the performance of the approach shown here with the previously hand-coded pure C++ version described in [13]. As a test example we choose a standard Riemann problem for the Euler equations solved on a series of different grid resolutions using quadratic basis functions, a minmod limiter, and explicit RK3 time stepping. We use two different grid implementation, a dedicated Cartesian grid (SPGrid) and a fully unstructured grid (ALUGrid). The results are shown in Table 1.
We observe that for the Cartesian grid (SPGrid) using the Python frontend leads consistently to a small improvement of 2%. For the unstructured grid (ALUGrid) we observe a performance decrease of about 4-5%. This can be explained with the fact that for SPGrid all code can be inlined in the just-in-time compiled Python module. For ALUGrid, where a library exists, this is not so straight forward. In the future we will experiment with link time optimization and try to reduce implementation of small code snippets in the ALUGrid library.

Conclusion and outlook
In this paper we presented a comprehensive framework for the Discontinuous Galerkin (DG) Method with an easy to use Python interface for model description and solver configuration. The framework covers different variants of existing DG methods and can be easily extended to include improvements and new development of the methodology.
Although the DUNE-FEM-DG framework serves as a good starting point for highperformance DG solvers, the aim of this work is to demonstrate how the Python interface simplifies the simulation of a wide range of evolution equations as well as enabling the rapid prototyping of new methods, their testing and their comparison with other methods. For example, in almost all the tests presented here, computations were done using the same setting for the indicator and tolerance for adaptivity and the troubled cell detection, demonstrating the robustness of the default setup. Of course the default setup will only be a suitable starting point and so most of the building blocks of the discretization can be straightforwardly replaced from within the Python script, often without requiring any or not much C++ knowledge. Also we would like to point out the implemented DG solver can be directly used as higher (2nd) order Finite Volume (FV) scheme by replacing the DG space with a FV space consisting of piecewise constant polynomials.
The Python interfaces are not yet fully reflecting the possibilities on the C++ side, for example, the reconstruction in primary variables for the Euler equations is possible, but not yet available on the Python side. Another missing feature is that it is currently not possible on the Python side to easily exchange the average operators in the discretization of the diffusive terms as, for example, needed in some applications [14].
Another missing feature in Python is the assembly of Jacobian matrices during the nonlinear solve. This could be desirable for some applications or to test different existing solvers. Although the feature is available in the C++ code it needs a few code alterations in the underlying infrastructure package. This would also make it straightforward to then use other solver packages available in Python, e.g., scipy. For the compressible applications, which were the focus of this work, it is much more feasible to work with the Jacobian free nonlinear solvers and there robust and efficient preconditioning techniques are still a very active research topic. We are currently focusing on including ideas presented in [6] based on subcell finite volume multigrid preconditioning.
The extension of the DGSEM methods to simplicial grids (see [10] for an overview) could be implemented with a few minor modifications. Based on earlier work [21] the necessary basis function implementations for arbitrary quadrature points are available but need to be integrated into the DUNE-FEM and DUNE-FEM-DG framework.
While DUNE-FEM provides a number of Runge-Kutta methods, we have shown here that it is straightforward to add other time stepping algorithms on the Python side, a feature which should also allow to use other packages which provide bindings for Python, such as the Assimulo package [2]. For IMEX schemes the splitting is at the moment still a bit restricted focusing on the important case of implicitly treating the diffusion (and part of the source term) while using an explicit method for the advection. In some applications other types of splitting are of interest and will be made available in future releases.

A Interfaces
Signature of femdgOperator: def femDGOperator ( Model , space , limiter = " default " , advectionFlux = " default " , diffusionScheme = " cdg2 " , parameters = None ) The Model class contains the description of the mathematical model to solve and is described below. The space parameter is one of the available Discontinuous Galerkin spaces available in the DUNE-FEM framework. We have limiter equal to None,"minmod","scaling", taking advectionFlux equal to "default" to use Local-Lax-Friedrichs scheme, for Euler there are some other fluxes available, this parameter can also be used to pass in a user defined flux implementation. The parameter diffusionScheme can be used to change the diffusive flux. The parameters is a dictonary where additional information can be passed to the C++ code, e.g., the tolerance for the troubled cell indicator. Static methods on Model class passed to femdgOperator: The final ingredient is the time stepper: def femdgStepper ( * , order , operator , rkType = None , cfl = 0 . 45 , parameters = None ) If rkType is None or "default" then the type of Runge-Kutta method used will depend on the methods defined on the Model class used to construct the operator. If no diffusive flux F_v was defined an explicit RK method is used, if a diffusive flux but no advective flux F_c is defined an implicit method is used, while if both fluxes are present an IMEX scheme is employed. Again the parameter dictonary can be used to set further parameters read by the underlying C++ code.

B C++ code for modal troubled cell indicator
The full class for the modal indicator. The class needs to derive from a templated pure virtual base class Dune::Fem::TroubledCellIndicatorBase<DiscreteFunction> and override a single method. The main work is done in the private method smoothnessIndicator which sets up a least square problem to compute a smoothness indicator following [35]. We assume that the degrees of freedom uEn[i] are modal.