Robust interior penalty discontinuous Galerkin methods

by   Zhaonan Dong, et al.
University of Leicester

Classical interior penalty discontinuous Galerkin (IPDG) methods for diffusion problems require a number of assumptions on the local variation of mesh-size, polynomial degree, and of the diffusion coefficient to determine the values of the, so-called, discontinuity-penalization parameter and/or to perform error analysis. Variants of IPDG methods involving weighted averages of the gradient of the approximate solution have been proposed in the context of high-contrast diffusion coefficients to mitigate the dependence of the contrast in the stability and in the error analysis. Here, we present a new IPDG method, involving carefully constructed weighted averages of the gradient of the approximate solution, which is shown to be robust even for the most extreme simultaneous local mesh, polynomial degree and diffusion coefficient variation scenarios, without resulting in unreasonably large penalization. The new method, henceforth termed as robust IPDG (RIPDG), offers typically significantly better conditioning than the standard IPDG method when applied to scenarios with strong mesh/polynomial degree/diffusion local variation. On the other hand, when using uniform meshes, constant polynomial degree and for problems with constant diffusion coefficients, the RIPDG method is identical to the classical IPDG. Numerical experiments indicate the favourable performance of the new RIPDG method over the classical version in terms of conditioning and error.



There are no comments yet.


page 16


A p-robust polygonal discontinuous Galerkin method with minus one stabilization

We introduce a new stabilization for discontinuous Galerkin methods for ...

Residual-based a posteriori error estimates for 𝐡𝐩-discontinuous Galerkin discretisations of the biharmonic problem

We introduce a residual-based a posteriori error estimator for a novel h...

A posteriori error analysis of a local adaptive discontinuous Galerkin method for convection-diffusion-reaction equations

A local adaptive discontinuous Galerkin method for convection-diffusion-...

A discontinuous Galerkin method by patch reconstruction for convection-diffusion-reaction problems over polytopic meshes

In this article, using the weighted discrete least-squares, we propose a...

Implementation of high-order, discontinuous Galerkin time stepping for fractional diffusion problems

The discontinuous Galerkin dG method provides a robust and flexible tech...

Moving Mesh with Streamline Upwind Petrov-Galerkin (MM-SUPG) Method for Convection-Diffusion Problems

We investigate the effect of the streamline upwind Petrov-Galerkin metho...
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

Discontinuous Galerkin (DG) methods are a staple in the cascade of Galerkin approaches for the numerical solution of PDEs in divergence form. For the discretization of elliptic operators, the class of interior penalty DG (IPDG) methods is arguably the most popular choice due to its versatility and simplicity in implementation. The consistency of the IPDG method is achieved by the inclusion of specific terms involving integrals of solution fluxes on element interfaces. The basic, yet revolutionary, idea in IPDG methods is the further addition of penalization terms involving the jump of the approximate solution (or its derivatives for higher order elliptic PDEs) to control the departure from the PDE solution space and, therefore, to achieve stability in a consistent fashion. The IPDG method design concept was presented in full by Baker [3] and, briefly later, by Wheeler [22] and Arnold [1], thereby providing a consistent alternative to the classical inconsistent finite element method with penalty of Babuška [2]. These classical IPDG methods and their -version variants [19, 16] are widely used to this day. An important generalization of this classical setting, due to Dryja [10], Burman & Zunino [5] and Ern, Stephansen & Zunino [11], is the use of weighted averages on the normal flux functions in the IPDG definition to counteract high-contrast diffusion coefficients.

To fix ideas, we shall focus on the basic second order elliptic boundary value problem on an open domain , : find such that


with data ,

and a positive definite diffusion tensor

. Although considerable extensions to the problem can be incorporated immediately to the developments presented below, we refrain from exhausting the generalization possibilities in the interest of clarity of exposition.

Classical IPDG methods for the problem (1.1) involve a user-defined quantity, the so-called penalty or discontinuity-penalization parameter. The judicious choice of this parameter is crucial for the stability of the IPDG method, but also for its practical behaviour regarding approximation quality and conditioning. In particular, it is well known that excessive penalty parameter typically leads to increased ill-conditioning of the IPDG stiffness matrix. At the same time, the penalty parameter has to be chosen large enough to ensure stability of the method. In this work, we focus on the in the behaviour of interior penalty methods in the presence of highly locally variable meshes, local polynomial degrees, and/or diffusion coefficients . As we will see below, in such extreme local variation scenarios of discretization parameters and/or coefficients, the classical IPDG method typically requires excessively large penalization to retain stability. In fact,

it is possible that the classical IPDG penalty parameter degenerates to infinity while the dimension of the approximation spaces (i.e., the cardinality of the numerical degrees of freedom) remains finite

. This, in turn may be detrimental to the conditioning of the problem and, in certain cases, to the quality of the approximation.

To resolve this state of affairs, we present a new, to the best of our knowledge, weighted-type interior penalty discontinuous Galerkin method which is provably stable with a new choice of the discontinuity-penalization parameter and is robust with respect to extreme local variations in the discretization parameters and in the diffusion coefficient. The method, termed henceforth as robust IPDG (RIPDG) replaces the classical average-of-the-normal-flux term(s) in the IPDG by suitable, explicitly constructed weighted averages. This modification results into a robust and concrete selection of the discontinuity-penalization parameter, without resulting into excessive penalization in extreme combined local mesh, polynomial degree and coefficient variation scenarios. We note that a balanced choice for the discontinuity-penalization parameter is particularly important when no conforming subspace is available. Moreover, the new RIPDG method allows for the first time in the interior penalty literature for the proof of a priori error bounds without any local-quasi-uniformity/local bounded variation assumptions for the discretization parameters. As a result, the corresponding RIPDG penalty parameter remains finite for fixed finite dimension of the underlying approximation space.

As we will see below, the modifications required to implement RIPDG starting from an available IPDG code are minimal and trivial. We envisage that the RIPDG method will be highly effective in the numerical approximation of complex, multiscale problems characterized by extreme local physical features necessitating highly non-uniform Galerkin approximation spaces. At the other end of the spectrum, when using uniform meshes, constant polynomial degree to solve problems with constant diffusion coefficients RIPDG is identical to the classical IPDG. Numerical experiments indicate the favourable performance of the new RIPDG method over the classical version in terms of conditioning and error. In particular, despite extensive testing, we have not been able to identify an example whereby the RIPDG method is inferior in terms of conditioning and/or error in various norms when compared to the classical IPDG method. On the contrary, in all cases considered the RIPDG method performs better than the classical IPDG in one or more of the above indicators of performance. Although we expect that there exist cases whereby IPDG outperforms RIPDG exist, we have not yet been able to identify such example despite extensive testing.

The remainder of this work is structured as follows. In Section 2 we review the classical IPDG method and the weighted variant from [10, 5, 12], along with the basic argument for proving coercivity and continuity of the respective bilinear form. In Section 3, we present the new RIPDG method and discuss its construction and properties, as well as the motivating differences in the choice of the penalty parameter. In Sections 4 and 5 we provide extensions of RIPDG to polytopic meshes and to degenerate elliptic problems, respectively. In Section 6, we present the basic Strang’s-type Lemma satisfied by IPDG and RIPDG, showcasing that in the latter case the resulting a priori error bound is independent of any local variation in meshsize, polynomial degree and diffusion contrast. We conclude with a series of numerical experiments showcasing the benefits in using weighted averages in scenarios with high such local variations.

2. Interior penalty discontinuous Galerkin methods

We consider meshes consisting of mutually disjoint open simplicial and/or box-type elements so that . Let also the diameter of , and define the mesh-function by , . Let also denote the number of -dimensional faces of the element , i.e., for a -dimensional simplex and for a -dimensional box-type element, and define by , . Further, we let denote the mesh skeleton and set . The mesh skeleton is decomposed into –dimensional simplicial denoting the mesh faces, shared by at most two elements. These are distinct from elemental interfaces , herein defined as the simply-connected components of the intersection between the boundaries of two neighbouring elements. As such, hanging nodes/edges are naturally permitted: an interface between two elements may consist of more than one faces. Let also denote the broken gradient defined as , .

The discontinuous Galerkin space with respect to is defined as

with denoting the space of -variate polynomials of total degree up to on , and with , . The local elemental polynomial spaces employed within are defined in the physical coordinate system, i.e., without mapping from a given reference or canonical frame (cf.  [6, 7, 8]). All developments discussed below are also valid for the more ‘classical’ case of the discontinuous Galerkin spaces constructed through the usual mapping from a suitable reference element , viz.,

Let and be two adjacent elements of sharing a face . For and

element-wise continuous scalar- and vector-valued functions, respectively, we define the

weighted average across by

for with , respectively, with denoting the trace of from the element on , and correspondingly for . Also, we define the jump across by

with denoting the unit outward normal of on . We collect all weight pairs in a weight function , with , for each face . On a boundary face , we set and , respectively, with the unit outward normal to . In view of the latter, we extend to by setting with there.

The (weighted) interior penalty discontinuous Galerkin method reads: find , (or in ,) such that




with and the, so-called, discontinuity-penalization or penalty function, whose precise definition, along with that of the weights is a central concern in this work. Note that, selecting , for all , we retrieve the classical IPDG method of Wheeler [22] and of Arnold [1], while the case or for all gives the original method of Baker [3]. Moreover, in [10, 5], the weighted average choice with

with , was proposed for the case of element-wise constant scalar diffusion coefficients, while in [11] the selection for element-wise constant diffusion tensors was provided, by combining the ideas from [10, 5] and [13]. These choices result to robust a priori error analysis with respect to the relative sizes of diffusion locally. The main contribution of this work is the proposition of a different recipe for , which will yield robust dependence not only with respect to the local variation of the diffusion, but also with respect to the local variations in meshsize and polynomial degree. This will, in turn, yield a robust and rigorous choice for the discontinuity-penalization parameter .

The choice in (2.1) yields the symmetric version, while the choice the non-symmetric version of the weighted IPDG method. Although, we shall focus on the symmetric version in this work, as the most popular choice, the developments presented below are also valid for the whole range , with trivial modifications of the arguments and of the formulas.

To highlight the importance in the careful selection of , we consider a -dimensional simplicial triangulation , we set , for all faces , and we study the coercivity of the bilinear form on . To that end, for , have



the usual arithmetic average used in standard IPDG methods. To further estimate the last term on the right-hand side of (

2.2) from below, we employ a trace inverse estimate, such as the one below.

Lemma 2.1.

Let a simplicial -dimensional simplex and one of its faces . For any , we have


This follows immediately by combining a careful scaling argument with the main result from [21]. ∎

We refer to [8] for a generalization of the above estimate to general curved polygonal/polyhedral elements with arbitrary number of faces. Note that for shape-regular elements, (2.3) yields the familiar version with the constant depending on the aspect ratio of the element and the dimension . If, instead, we have , then will also depend on the elemental maps if the latter are non-affine.

Returning to (2.2), elementary calculations, along with (2.3), for faces give

with ,

and the number of faces of the element , defined above; note that , i.e., the discontinuous Galerkin space defined by the same mesh having polynomial basis of one degree less than .

These developments, in conjunction with (2.2), imply

upon setting . Therefore, the selection , for all faces gives the coercivity estimate


The above choice of the discontinuity-penalization parameter is, to the best of our knowledge, the sharpest general estimate for simplicial meshes. An analogous choice of was proposed in [6, 8] for significantly more complex element shapes. We note that hanging nodes are permissible by viewing a simplex as a polytopic element with many co-planar faces and by modifying accordingly.

Let us define now the extension of the bilinear form for the special case at hand, to accept arguments from the larger space , viz., with

for all , with denoting the orthogonal -projection operator onto . Completely analogous arguments to the ones presented for the proof of (2.4), along with the stability of in the -norm, result in the continuity bound

Remark 2.2.

If a different Galerkin space than is used, one has to modify the coercivity and continuity analysis and, correspondingly, the definition of accordingly. In particular, it may not be the case anymore that , but rather that resides on the Galerkin space itself.

3. A robust interior penalty discontinuous Galerkin method

We now propose a robust interior penalty discontinuous Galerkin method by taking advantage of appropriately selected weights in (2.1). In particular, we will select diffusion tensor, mesh and local polynomial degree dependent weights. To that end, upon considering (2.1), we have


Setting , using (2.3), and working as before, we get, for ,

Setting now

for , and selecting


we deduce

A reasonable selection for the discontinuity-penalization parameter is, therefore,


with in the case of a boundary face.

Definition 3.1.

The robust interior penalty discontinuous Galerkin method (RIPDG) is the Galerkin procedure defined by (2.1) with for given by (3.2), and the discontinuity-penalization parameter selected as in (3.3).

Remark 3.2.

The above argument still applies for different choices of element-wise polynomial Galerkin spaces by replacing suitably by an available upper bound of the respective inverse estimate constant for the space at hand.

Remark 3.3 (On inverse estimate constants).

For different choices of element-wise Galerkin spaces, we may end up with unspecified/hard to estimate universal constants in the inverse estimate. We note carefully that any such constants cancel out in the definition of the weights and, thus, the weights are independent from such constants, ensuring the practical selection of weights in various approximation space scenarios. Of course, as is the case for the standard IPDG method, the RIPDG penalty parameter will depend proportionally on such universal constants.

The proof of (2.5) is also then immediate upon considering the inconsistent formulation

for all , resulting in the continuity bound


We now compare the effects on the size of of the classical choice and the of the one proposed herein for the averaging operator. Any significantly different behaviour is bound to occur in extreme mesh and/or polynomial degree local variation scenarios. To that end, we consider the problem (1.1) for with

, the identity matrix. The latter’s solution is approximated via (

2.1) for and also via (2.1) with the choice (3.2), over a mesh consisting of two quadrilateral elements

for and with uniform local polynomial degree . The discontinuity penalization parameter will differ for the two variants only on the sole interior face . For the classical IPDG method ((2.1) with ), we calculate:

whereas for the choice (3.2), we get

So, as , we have , while This extreme scenario highlights vastly different penalization pattern between the two variants. In the new choice of the weighted version of IPDG ((2.1) with (3.2)), we observe robust behaviour with respect to the local mesh variation.

Completely analogously, we now set , so that and are identical in terms of shape, and we set and . We then compute

as , whereas

Having the penalty parameter converging to infinity may become an issue in terms of the conditioning of the resulting linear systems. Therefore, when a conforming subspace of the DG space with the same approximation properties exists, e.g., for the case of a simplicial mesh, one expects only conditioning differences between the two variants. However, in the more practically interesting case whereby there is no conforming subspace of full order, e.g., for the case of on a mesh consisting of box-type (quadrliateral/hexahedral) elements, may be compromising in terms of approximation also.

4. Extension to general polytopic elements

The classical IPDG has been put forward as a method of choice for the development of Galerkin methods on meshes consisting of extremely general polygonal/polyhedral elements; see, e.g., [6, 7, 8, 4] and the references therein. A key development for practical applications has been their ability to admit provably stable approximations on elements with arbitrary number of faces , possibly with . The weight selection in the robust IPDG method presented above depends on the number of elemental faces. Applying the method in the form (2.1) with the selection (3.2) is may lead to spurious, excessive over-penalization as .

Fortunately, it is often possible to use non-overlapping simplicial subdivisions of the element to effectively arrived at a bounded . More specifically, a number of extremely mild geometric conditions are presented in [8] (cf. also [7] for earlier ideas in this vein) so that a polygonal/polyhedral element , possibly containing curved faces, admits an inverse estimate of the form


for any and for a mutually disjoint subdivision of , with a concretely determined expression; see [8, Lemma 4.21] for the exact formula. A key idea in [8] is that each is star-shaped with respect to one interior point of , , and may contain an arbitrary number of (possibly very small in size) -dimensional faces. Following the same line of argument as in Section 3, we can define a robust IPDG method on meshes consisting of general, possibly curved, polygonal/polyhedral elements by selecting

and retaining the definitions of (3.2) and (3.3). We note carefully that in the case of general polytopic meshes (4.1) may contain an unknown constant which, nevertheless, cancels out in the definition of the weights as discussed in Remark 3.3, yielding explicit choice for the weights .

5. Robust IPDG for degenerate elliptic problems

It is possible to have well-posed elliptic problems of the form (1.1) with diffusion tensors for which , or even , on lower-dimensional manifolds of [18]. In such cases, the choice of weights from (3.2) has to be revisited. To that end, we consider the inconsistent robust interior penalty discontinuous Galerkin method reading: find , (or in ,) such that




for , with (or, ,) denoting the orthogonal -projection operator onto the dG space. The method was first presented in [13] for the case and with replaced by . For simplicity of the presentation, we shall again focus on the symmetric case ; the cases follow via trivial modifications of the arguments given below.

From (2.3), and working as before, we get, for ,

using the stability of the orthogonal -projection operator. Setting now

for , and selecting, again,


we deduce


for every face . If the limit (from within the element interior to the face) is valid a.e. on for exactly one of , we have , with the other weight in the pair tending to zero; also, in this case we have .

Further, if we have from within the element interior to the face for both traces a.e. on , we have and on .

6. A priori error analysis

We now show that standard a priori error bounds hold for all the robust IPDG methods described in the previous section, assuming only sufficient regularity of the exact solution. In particular, no ‘local bounded variation/local quasi-uniformity’ mesh and polynomial degree assumptions are required, as is standard in the error analysis of the classical IPDG method with .

More specifically, let , . Then, the coercivity, the continuity along with standard, elementary calculations imply


with for the robust IPDG from Section 3, or for the variant from Section 5.

To estimate the inconsistency term, assuming further that , we make use of the best approximation estimate from [9] (see also [16, 14] for earlier results on box-type elements): there exists a constant , independent of and of such that


for any face of a reference simplicial element . Employing a standard scaling argument, (6.2) implies for an element the following best approximation estimate:


Upon observing the identity (or for the method of Section 5), we have, respectively, for any ,

using (6.3) in the first step and the definition of in the penultimate step. The last estimate concerns the robust IPDG method from Section 3 and holds under the regularity assumption .

Correspondingly, for the method of Section 5, assuming that the is such that , (instead of as before,) we arrive at

Combining the above developments with (6.1), we deduce


for the robust IPDG of Section 3, or