Acknowledgements

Adaptive discontinuous Galerkin methods

for nonlinear parabolic problems

Thesis submitted for the degree of

Doctor of Philosophy

at the University of Leicester

by

Stephen Arthur Metcalfe MMath

Department of Mathematics

University of Leicester

2014

Adaptive discontinuous Galerkin methods

for nonlinear parabolic problems

Stephen Arthur Metcalfe


This work is devoted to the study of a posteriori error estimation and adaptivity in parabolic problems with a particular focus on spatial discontinuous Galerkin (dG) discretisations.

We begin by deriving an a posteriori error estimator for a linear non-stationary convection-diffusion problem that is discretised with a backward Euler dG method. An adaptive algorithm is then proposed to utilise the error estimator. The effectiveness of both the error estimator and the proposed algorithm is shown through a series of numerical experiments.

Moving on to nonlinear problems, we investigate the numerical approximation of blow-up. To begin this study, we first look at the numerical approximation of blow-up in nonlinear ODEs through standard time stepping schemes. We then derive an a posteriori error estimator for an implicit-explicit (IMEX) dG discretisation of a semilinear parabolic PDE with quadratic nonlinearity. An adaptive algorithm is proposed that uses the error estimator to approach the blow-up time. The adaptive algorithm is then applied in a series of test cases to gauge the effectiveness of the error estimator.

Finally, we consider the adaptive numerical approximation of a nonlinear interface problem that is used to model the mass transfer of solutes through semi-permiable membranes. An a posteriori error estimator is proposed for the IMEX dG discretisation of the model and its effectiveness tested through a series of numerical experiments.

Acknowledgements

Foremost, I wish to thank my supervisors Dr. Andrea Cangiani and Dr. Emmanuil Georgoulis for their invaluable advice and guidance throughout the last four years. Additionally, I would like to thank Dr. Irene Kyza for inviting me to collaborate with her in Crete and for her insights into blow-up problems. I would also like to thank everyone in the mathematics department at the University of Leicester for the many fruitful discussions that helped contribute to this work. Finally, I acknowledge the funding of the EPSRC; without their support this work would not have been possible.

Contents

Chapter 1 Introduction

Partial differential equations are key in the modelling of various physical and biological phenomena. The solutions to PDEs are usually unavailable through analytical means, so numerical methods are employed in order to approximate the solution. Furthermore, a number of nonlinear PDE problems exhibit local multiscale behaviour such as boundary or interior layers, interfaces or even local space-time blow-up. Such local multiscale features require high local resolution of the numerical methods employed in their approximation. Hence, the use of numerical methods which can automatically detect and resolve such multiscale features is of interest.

Adaptive algorithms that are driven by a posteriori error estimators lie at the heart of finite element analysis. For elliptic problems, there are a wide variety of different error estimators available [3, 105]. Moreover, adaptive algorithms for elliptic problems are relatively well understood; at least for simple elliptic problems see [28, 37, 82] for the standard conforming finite element method and [18, 61] for the interior penalty dG method.

For linear parabolic problems, there are many error estimators available in the literature for popular discretisations (typically a standard time stepping scheme paired with a spatial finite element discretisation). These error estimators are usually composed of an initial condition estimator, a space estimator and a time estimator. However, generally speaking, it is unclear how to utilise each of these individual estimators to drive adaptivity. Furthermore, while mesh change is crucial for the efficient numerical approximation of mobile solutions, it is well known that careless mesh refinement and/or coarsening can lead to destabilisation of the finite element solution [12, 39]. While some progress has been made on the construction of adaptive algorithms for parabolic problems [30, 40, 41, 42, 43, 44, 68, 88, 90, 96], none of the algorithms in the literature have been shown to reduce the error estimator that they utilise at the correct rate of convergence with respect to the average number of degrees of freedom and the total number of time steps.

This work will investigate adaptive algorithms for spatial discontinuous Galerkin discretisations of parabolic problems with a focus on nonlinear problems. To achieve this, we shall look at adaptivity in the context of three different problems, to be detailed below.

In Chapter 2, we introduce notation and state some approximation results and general theorems that shall be used throughout this work. We also discuss the issue of robustness that arises in the a posteriori error estimation of discontinuous Galerkin discretisations of stationary convection-diffusion equations. We then introduce a robust error bound for the stationary problem, taken from [97], that shall be used extensively in this work.

Chapter 3 deals with linear non-stationary convection-diffusion equations. In particular, we derive an a posteriori error estimator for a backward Euler dG discretisation of the problem. The primary challenge in such a derivation is the robustness of the error estimator with respect to the diffusion parameter . We address this through the elliptic reconstruction framework of Makridakis and Nochetto [81] paired with a robust error estimator for the stationary problem [97] together with a robust treatment of the temporal residual. The proposed error estimator can be viewed as the analogue of that given in [104] by Verfürth but with a discontinuous Galerkin spatial discretisation instead of a streamline upwind Petrov-Galerkin spatial discretisation. As well as [104], there are other error estimators available in the literature for different discretisations of linear non-stationary convection-diffusion equations. In particular, we mention [36] wherein the authors produce a robust error estimator through a flux reconstruction approach and the work of Picasso and Prachittham [89] wherein they develop an error estimator for a Crank-Nicolson temporal discretisation; their discretisation and error estimator also provide for the use of anisotropic elements. There are also a variety of different a posteriori error estimators available for the norm, cf. [45, 63, 99]. In addition to the error estimator we also propose an adaptive algorithm, based on that given in [30], that utilises different parts of the error estimator to control space and time adaptivity. The effectiveness of both the error estimator and the proposed algorithm is then tested in four numerical experiments. It is also worth noting that adaptive algorithms designed specifically for non-stationary convection-diffusion problems are explored in [89, 99].

In Chapter 4, we investigate the numerical approximation of blow-up in ODEs. More specifically, we derive an a posteriori error estimator for an ODE with polynomial nonlinearity that is discretised using standard time stepping schemes. The biggest difficulty in the construction of such an error estimator is having to deal with a nonlinear error equation this can be handled through a local continuation argument. A continuation argument is a special type of proof by contradiction that is often used to prove existence results for nonlinear parabolic PDEs; such arguments have been instrumental in the derivation of a posteriori bounds for a variety of nonlinear parabolic problems [13, 57, 74]. A posteriori error estimators produced by a continuation argument are conditional in the sense that they only hold providing that the estimators involved are sufficiently small. In order to use the proposed error estimator to approximate the blow-up time, we investigate the design of suitable adaptive algorithms. To that end, two adaptive algorithms are proposed and then applied in two test cases under different time stepping schemes to compare their effectiveness. Although we choose to investigate the numerical approximation of blow-up through a posteriori error estimation, there are other ways of approaching this problem that have been published in the literature. In [66], the authors prove existence results for numerical approximations to a nonlinear ODE with a polynomial growth condition provided that the time step lengths are sufficiently small. For the particular case of a polynomial nonlinearity, they show that selecting the time step lengths in a certain way yields approach to the blow-up time. In [60], the authors transform an ODE with polynomial nonlinearity through an arc length transformation. They then use a forward Euler method to approximate the transformed equation and they show that their adaptive algorithm, which is based on their transformation plus a tolerance controlled ODE integrator, converges towards the blow-up time linearly with respect to the total number of time steps. Finally, in [98] the authors approximate a nonlinear ODE using a -method along with a temporal rescaling of the ODE and they show that their numerical solution has the same asymptotic behaviour as the exact solution.

In Chapter 5, we build on the results of the previous chapter by investigating blow-up in semilinear parabolic PDEs. In particular, we study blow-up in nonlinear non-stationary convection-diffusion equations that, based on the results of the previous chapter, are discretised using an IMEX dG method. An a posteriori error estimator is then proposed for this discretisation of the problem. In order to produce a viable error estimator for this problem, there are three major obstacles that need to be overcome: the lack of symmetry of the problem, the nonlinear error equation and the error due to non-conformity. A posteriori error estimation for blow-up in symmetric problems has been considered in [78, 79] by Kyza and Makridakis, however, such results are not easily generalised to non-symmetric problems; the lack of symmetry in the problem can, however, be dealt with through the Gagliardo-Nirenberg inequality. The nonlinear error equation is dealt with in a similar way to in Chapter 4 through a continuation argument and the error due to non-conformity can be dealt with via localised bounds for the non-conforming part of the error [34, 70, 71]. The proposed error estimator is utilised to approximate the blow-up time of the problem through an adaptive algorithm that is based on those given in Chapter 3 and Chapter 4. The adaptive algorithm is then applied to some test problems and the results are compared to those given in Chapter 4. It is worth noting that solution profiles close to the blow-up time can also be obtained through the rescaling algorithm of Berger and Kohn [16, 85] or the MMPDE method [20, 65]. There is also work looking at the numerical approximation of blow-up in the nonlinear Schrödinger equation and its generalisations [4, 31, 50, 75, 101]. Other numerical methods for approximating blow-up in a variety of different nonlinear PDEs can be found in [8, 33, 35, 49, 84].

In Chapter 6, we consider an IMEX dG discretisation of a nonlinear interface problem that was introduced in [25, 26], based on the works [51, 73, 91, 100, 114], to model the mass transfer of solutes through semi-permiable membranes. There are a variety of a posteriori error estimators of both residual and recovery type in the literature for different discretisations of interface problems. In particular, for the conforming finite element method [22, 23], the discontinuous Galerkin method [21] and the finite volume method [47, 83]. Based on these works, and the techniques used in previous chapters, we seek to derive a residual-based a posteriori error estimator for this discretisation of the model. The effectiveness of the error estimator is then tested through a series of numerical experiments utilising the adaptive algorithm that was developed in Chapter 3.

Finally, in Chapter 7 we summarise the results of this work and discuss ways in which this work could be extended.

Chapter 2 Preliminaries

2.1 Sobolev spaces

Let be a bounded Lipschitz domain with boundary . For , we define the norms by

and the respective spaces by

Note that is a Hilbert space with an inner product given by

When is the computational domain (to be defined later) then because both the norm and inner product over are used frequently in this thesis, the relevant subscripts are omitted. Given a multi-index , the weak derivative of order is given by

For and , the Sobolev space is given by

The spaces are equipped with the norms

The space together with the standard inner product is a Hilbert space which we shall denote by . Fractional Sobolev spaces ( in particular) are also of significant use when trying to make sense of boundary values, in the sense of traces, in Sobolev spaces; we refer to [1] for details. Whenever boundary values are used in this thesis, they are to be understood in the sense of traces. We define the space , which is the prototypical PDE solution space, by

where is some subset of with positive one-dimensional Hausdorff measure; if , we denote this space by . Finally, we let denote the space of all functions for which is continuous for all multi-indices with .

For , the spaces (where is a real Banach space with norm ) consist of all measurable functions for which

We also define . Finally, we denote by and , respectively, the spaces of continuous and Lipschitz continuous functions such that

2.2 Stationary convection-diffusion equation

Let the computational domain be a polygon with boundary , this assumption will be used throughout the rest of this thesis. We consider the model problem of finding such that

(2.1)

The variable functions are collectively referred to as the data of the problem. Problem (2.1) is the prototypical convection-diffusion equation. If the convection is constant then the scale of the solution to (2.1) can be characterised through the ratio of convection to diffusion as described by the Péclet number

When , (2.1) is advection dominated and can exhibit some or all of the following features (see [69, 93] for a detailed analysis):

  • The presence of ordinary layers, spatial areas containing steep gradients of the solution of width that usually occur near the outflow boundary as boundary layers.

  • The presence of parabolic layers, spatial areas containing moderate gradients of the solution of width that typically occur on the inflow boundary as boundary layers or as interior layers.

These complex spatial features can render the numerical approximation of the solution difficult as discussed in the next section.

In order to analyse (2.1) and its parabolic counterparts, we must make some assumptions on the data. We assume that: , , and . Furthermore, we require some additional assumptions in order to state standard coercivity and continuity results (see, e.g., [97, 110]). To that end, we assume that there are constants and such that

(2.2)

The weak form of (2.1) reads: find such that

(2.3)

where

(2.4)

2.3 Discontinuous Galerkin method

A finite element method is a numerical technique for finding approximate solutions to the weak formulation of PDEs characterised by the use of a subdivision of referred to as the mesh or triangulation. The mesh is a collection of elements with denoting a generic element. A finite element space is then constructed over the mesh and the weak form of the PDE is discretised; different choices of finite element space and different discretisations give rise to different finite element methods. The discretisation parameters are quantities related to convergence of the method, specifically, the diameters of elements in the mesh (and the lengths of time steps in parabolic problems).

When it comes to the finite element approximation of (2.3), the presence of layers introduces a certain amount of difficulty. In particular, the standard conforming finite element method performs poorly if an insufficient number of elements are placed in the vicinity of the layers resulting in unphysical oscillations. This issue can be solved with layer-adapted meshes such as Shishkin meshes (see [76] for an overview of this subject) but these special meshes require a priori knowledge of where the layer will occur. These problems led to the development of stabilised finite element methods for convection-diffusion equations [93] the most popular of which are the streamline upwind Petrov-Galerkin (SUPG) method [19] and the discontinuous Galerkin (dG) method. In this work, we focus upon a dG discretisation of (2.3) taken from [62] which is based upon a classical interior penalty discretisation of the diffusive term originally introduced in [7, 9, 87] and an upwind discretisation of the transport term first discussed in [80, 92].

In order to state the dG discretisation of (2.3), we need some additional notation. The mesh is assumed to be constructed via affine mappings with non-singular Jacobian where is the reference triangle or the reference square. The mesh is allowed to contain a uniformly fixed number of regular hanging nodes per edge. We define the finite element space

(2.5)

where is the space of polynomials of total degree if is the reference triangle, or the space of polynomials of degree in each variable if is the reference square. Let denote the set of all edges in the mesh and the set of all interior edges. We also denote the diameter of an element by and the length of an edge by . The outward unit normal to the boundary of an element is denoted by . We assume that the mesh is shape-regular, that is, there exists such that for all we have

where denotes the diameter of the largest ball that can be completely contained in .

In what follows, it will be useful to associate patches with each element . Specifically, we have the (elemental) patch which is the union of all elements that “neighbour” and the edge patch which is the union of all edges that intersect the boundary of . Formally, these are defined, respectively, by

We also associate an (elemental) patch with each edge which is the union of all elements whose boundary intersects , viz.,

Given an edge shared by two elements and , a vector field and a scalar field , we define jumps and averages of and across by

If , we set , , and , with denoting the outward unit normal to the boundary .

We define the inflow and outflow parts of the boundary , respectively, by

Similarly, the inflow and outflow parts of an element are defined as

With all the above notation at hand, the dG approximation to (2.3) reads as follows: find such that

(2.6)

where

(2.7)

The penalty parameter, , is set to in light of [62] so that the operator is coercive on (see below).

We note that the bilinear form is not well-defined for arguments in , but the bilinear form is and is equal to that appearing in (2.3). To analyse the dG discretisation, we introduce the quantities

These quantities define norms on . In the literature, is referred to as the energy norm while the quantity is referred to as a dual norm. It is easy to see that the bilinear form is coercive on , viz.,

(2.8)

for all , and is continuous in the following sense

(2.9)

for all and . Moreover, the discrete bilinear form is coercive for with respect to the energy norm, viz.,

(2.10)

The symbols and used above and throughout the rest of the thesis are used to describe inequalities that are true up to an unspecified positive constant that is independent of the data, the discretisation parameters, the exact solution and the dG solution.

2.4 Error bounds for the stationary problem

Let be the exact solution of a PDE and be some finite element approximation; an a posteriori error estimator, , is an approximation of the error, , in a certain norm such that . The error estimator must be computable and thus is allowed to depend upon the data, the discretisation parameters and the finite element solution but not the unknown solution . In order to discuss how good is at approximating , it is useful to introduce the notion of reliability and efficiency of an estimator. An a posteriori estimator, , is said to be reliable if there is , independent of the exact solution , such that

(2.11)

while is said to be efficient if there is , independent of the exact solution , such that

(2.12)

The constants appearing in (2.11)-(2.12) are often impossible to calculate explicitly which leads us to the useful notion of the effectivity index. If is known for particular data, may be calculated explicitly for different realisations of . Thus, we can compute the effectivity index - the ratio of to :

effectivity index

The effectivity index naturally leads to the notion of robustness. An error estimator, , is said to be robust (with respect to ) if the constants in (2.11)-(2.12) are always independent of the data, the discretisation parameters and the finite element solution . There is also the weaker notion of asymptotic robustness: is said to be asymptotically robust if it is robust once the discretisation parameters are sufficiently small. If is asymptotically robust then it successfully reproduces the convergence rate of with respect to the discretisation parameters.

Figure 2.1: Numerical approximation of a boundary layer in the pre-asymptotic regime (left) and the asymptotic regime (right).

Regarding the a posteriori error estimation of (2.3), many estimators exist for stabilised finite element schemes [5, 6, 17, 46, 54, 67, 77, 95, 97, 108, 110, 113] and the primary issue is robustness with respect to the small parameter . Before the layers have been covered by a sufficient number of elements, the pre-asymptotic regime, the error in the energy norm is relatively constant. When a sufficient number of elements have been placed within the layers, the asymptotic regime, the error in the energy norm starts to display the correct convergence rate with respect to the discretisation parameters. In order for an error estimator to be robust with respect to in the energy norm, it must accurately capture the behaviour of the energy norm in both the pre-asymptotic and asymptotic regimes. Standard, classical estimators in the literature significantly overestimate the error in the energy norm in the pre-asymptotic regime to allow the layers to be detected and refined but this means that they are not robust with respect to in the energy norm until the layers have been sufficiently resolved. One way out of this difficulty is to add an additional term, the dual norm, to the energy norm to account for the pre-asymptotic regime [94, 95]; robustness is then recovered in all regimes for the full norm [97, 110]. The question of whether a robust error estimator exists for the energy norm in all regimes is still open, however, recent results in this direction seem promising [53].

An a posteriori estimator for the stationary problem, inspired by [97], will be utilised in our analysis. More specifically, we have the following result whose proof is completely analogous to Theorem 3.2 in [97] and is therefore omitted for brevity.

Theorem 2.1 ().

For , let be such that

and consider such that

Then the following a posteriori bound holds for any :

2.5 Finite element approximation results

Error bounds for the approximation of functions in have been constructed in the finite element literature and will be utilised below.

Theorem 2.2 ().

Given , there is a finite element quasi-interpolant such that

Proof.

See [109]. ∎

A useful tool in a posteriori error estimation involving dG finite element spaces is the approximation of functions in the dG finite element space, , by functions in the conforming finite element space, . The next theorem gives bounds on such approximation errors.

Theorem 2.3 ().

Given , there exists a decomposition with and such that the following bounds hold for each element :

Proof.

The proof is based upon a weighted averaging of the dG solution around the nodes; see [70, 71] for the first two estimates and [34] for the final estimate. ∎

Given a function , the projection of onto , denoted by , is the unique solution of

(2.13)
Remark 2.1 ().

For the dG finite element space, can also be constructed by considering (2.13) elementwise instead of globally and then piecing together the resultant local functions.

The final theorem in this section gives bounds on the approximation error involved in the projection.

Theorem 2.4 ().

Given a function and its projection , the following bound holds for any element :

Proof.

We use the definition of the projection along with the Cauchy-Schwarz inequality to conclude that for any :

Therefore,

All that remains is to choose to give the bound presented in the theorem. In particular, the finite element interpolant in Theorem 2.2 suffices. ∎

2.6 Useful inequalities

In this section, we introduce general theorems that will be used throughout the rest of this thesis.

Theorem 2.5 (Young’s inequality).

Given , then for any we have

Proof.

Follows by expanding the inequality . ∎

Theorem 2.6 (Power mean inequality).

Given , then for any we have

Proof.

Follows from Jensen’s inequality. ∎

Theorem 2.7 (Multidimensional integration by parts).

Given and we have

Proof.

Follows from the divergence theorem. ∎

Theorem 2.8 (Poincaré-Friedrichs inequality).

For , we have the following bound:

Proof.

See [112]. ∎

Theorem 2.9 (Gagliardo-Nirenberg inequality).

For any and , we have the following bound:

Proof.

See [2] specifically or [86] for a larger class of inequalities. ∎

Theorem 2.10 (Trace inequality).

Given , the following bound holds for any :

Proof.

See Theorem 1.5.1.10. in [58]. ∎

Theorem 2.11 (Inverse estimate).

Given , the following bound holds for any :

Proof.

See [59]. ∎

The final two theorems in this section relate to a specific group of integral inequalities and provide an upper bound on the functions that satisfy such inequalities. The first theorem is (classical) Gronwall’s inequality while the second theorem is a variant of the first theorem that makes use of a lower order term.

Theorem 2.12 (Gronwall’s inequality).

Let and suppose that and . If for almost every we have

then

where . Additionally, if and are non-negative a.e. then

Proof.

Multiplying the inequality by yields

The product rule and the fundamental theorem of calculus imply that

Thus,

The primary result then follows by integrating over and noting that . ∎

Theorem 2.13 ().

Let and suppose that is a constant, and are non-negative functions and that is a non-negative function that satisfies

then

Proof.

See Theorem 21 in [38]. ∎

Chapter 3 A posteriori error estimation and adaptivity for non-stationary convection-diffusion problems

3.1 Non-stationary convection-diffusion equation

For , we consider the model problem of finding such that

(3.1)

The model problem (3.1) can display the same spatial features as (2.1), namely, boundary and interior layers; the addition of the time domain has the potential to make things more complicated, however. Indeed, for suitable data the solution may exhibit moving boundary or interior layers. Additionally, the temporal behaviour of (3.1) can also be strongly influenced by [93].

The standard weak formulation (see, e.g., [48]) of (3.1) reads: find such that for almost every we have

(3.2)

where denotes the bilinear form (2.7) with the data evaluated at the time . We also make some assumptions on the data: , , , and .

In order to ensure coercivity and continuity of for any , we must extend (2.2). To that end, we assume that there are constants and such that

(3.3)

3.2 Space-time discretisation

The semi-discrete discontinuous Galerkin approximation to (3.2) then reads as follows. For , set to be some projection of onto . Then, seek such that for almost every we have

(3.4)

We shall also consider a full discretisation of problem (3.2) by using a backward Euler method to approximate the time derivative.

To this end, consider a subdivision of into time intervals of lengths , …, such that for some then set and . Denote an initial triangulation by and then associate a triangulation to each time step which is assumed to have been obtained from by locally refining and coarsening . This restriction upon mesh change is made to avoid altering the mesh too much between time steps in an attempt to prevent degradation of the finite element solution, cf. [12, 39]. To each mesh , we assign the finite element space given by (2.5). We also set , and for brevity.

The fully-discrete dG method then reads as follows. Set to be a projection of onto . For , …, , find such that

(3.5)

for all . We shall take to be the orthogonal projection of onto although other projections onto can also be used.

3.3 Error bounds for the non-stationary problem

Here, we will present a posteriori error bounds for both the semi-discrete and fully-discrete schemes which can be found in [27]; these results are extended by the authors in [72] to convection-diffusion problems with nonlinear reaction term. Other a posteriori error estimators for different space-time discretisations of (3.2) can be found in [36, 45, 63, 89, 99, 104].

To devise our error bounds, we use the elliptic reconstruction approach originally introduced by Makridakis and Nochetto for the conforming finite element method [81] and extended to dG methods in [55, 56]; this effectively allows decomposition of the error into separate parabolic and elliptic parts and can be viewed as the a posteriori counterpart to Wheeler’s elliptic projection [111] from a priori analysis. Elliptic reconstruction allows us to bound the elliptic part of the error with any error estimator for the stationary problem (2.3) that currently exists in the literature; we shall use the bound in Theorem 2.1 taken from [97]. Given that we already have a reasonable spatial estimator to use, it is clear that the main challenge in the a posteriori estimation of (3.2) is thus related to the parabolic part of the error. In particular, the primary challenge is obtaining an error estimator that is either robust or at least asymptotically robust with respect to in the type norm