Adaptive discontinuous Galerkin approximationsto fourth order parabolic problems

# Adaptive discontinuous Galerkin approximations to fourth order parabolic problems

## Abstract

An adaptive algorithm, based on residual type a posteriori indicators of errors measured in and norms, for a numerical scheme consisting of implicit Euler method in time and discontinuous Galerkin method in space for linear parabolic fourth order problems is presented. The a posteriori analysis is performed for convex domains in two and three space dimensions for local spatial polynomial degrees . The a posteriori estimates are then used within an adaptive algorithm, highlighting their relevance in practical computations, which results into substantial reduction of computational effort.

## 1 Introduction

Fourth order parabolic equations and corresponding initial-boundary value problems appear in the modelling in areas as diverse as biology, phase-field modelling and image processing to name a few. In most cases of practical interest one has to resort to numerical methods for their solution, due to complex geometry and/or the presence of non-linearities.

During the last five decades, finite element methods (FEMs) have been widely used to numerically solve fourth order elliptic or parabolic problems; see, e.g., [4, 15, 12, 16, 18, 10, 33] and the references therein for earlier works. There are, generally speaking, three families of FEMs developed for fourth order problems: conforming, mixed and non-conforming. The classical conforming methods (see, e.g., [15] and the references therein) require the construction of complicated elements with a number of degrees of freedom devoted to ensuring -continuity across the element interfaces. This results into limitations in the applicability of conforming methods on general, possibly irregular, meshes [36] and their non-trivial extensions to dimensions three (or higher). Mixed methods (see, e.g., [12, 16] and the references therein), whereby the fourth order operator is first transformed into a system of second order operators are widely used in practice, but they require very careful treatment in the imposition of essential and natural boundary conditions. Non-conforming methods for fourth order problems were first presented by [4] and then further developed in [18, 10, 33, 20] and other works. The key idea in non-conforming methods is the use of penalties to ensure convergence into the natural energy space of the variational problem, despite finite element basis functions being either just continuous (-interior penalty procedures; see, e.g., [18, 10]) or completely discontinuous (discontinuous Galerkin interior penalty procedures; see, e.g., [4, 33, 20, 21]).

Adaptive FEMs based on a posteriori error estimates has been an active field of research in recent years, especially for second order elliptic and parabolic problems. For the case of fourth order elliptic problems a posteriori error estimators and indicators have been developed, e.g., in [14, 38, 1, 34, 6, 13, 8, 27, 24]. A posteriori bounds and adaptive algorithms for parabolic fourth order problems are far less developed in the literature. For instance, the development of adaptive algorithms based on various types of a posteriori indicators for the Cahn-Hilliard fourth order parabolic problem can be found in [31, 22, 5]. Error control for variational methods for fourth order parabolic equations has been predominantly focused to space-discrete mixed or conforming formulations. The recent work [31] deals with goal-oriented error estimation for the fully discrete Cahn-Hilliard problem. Therefore, the development of adaptive algorithms based on a posteriori estimators for fully discrete methods for fourth order parabolic problems is still largely an unexplored area.

Advances in a posteriori error analysis of fully discrete schemes with non-conforming spatial discretizations of second order parabolic problems have been recently presented [19, 25]. In [25], an adaptive algorithm based on the derived a posteriori estimates is also considered. Local residual a posteriori error bounds for semi-discrete conforming and mixed spatial discretizations ffor the Cahn-Hilliard problem and the Hele-Shaw flow are presented in [22]. Finally, a posteriori error estimates in an -type norm and adaptive algorithms for fully discrete schemes with discontinuous Galerkin methods for fourth order problems are proposed in [39]. The derivation of reliability bounds in [39] is based on the elliptic reconstruction framework of Makridakis and Nochetto [32]; we also refer to [29, 25] for some relevant extensions.

This work is concerned with the derivation of a posteriori error estimates in weaker than -norms and their use within an adaptive algorithm for a class of discontinuous Galerkin interior penalty methods for a fully discrete approximation of the problem:

 ut+Δ2u = f in Ω×(0,T], (1) u=∇u⋅n = 0 in ∂Ω×(0,T] and (2) u = u0 in Ω×{0} (3)

with , a convex polygonal domain with boundary . More specifically, we derive a posteriori error estimators for the error measured in and norms for a numerical scheme consisting of discontinuous Galerkin method in space and simple backward-Euler time-stepping for the problem (1) - (3). The a posteriori analysis is performed for convex domains (as is usual for these norms) in two and three space dimensions for local spatial polynomial degrees . To enable the optimality of the a posteriori estimators in the and norms, the elliptic reconstruction framework is employed. Moreover, the -norm analysis employs a special test function construction inspired from the a priori analysis of FEMs for wave problems in [2]. Somewhat surprisingly, the use of this special testing, in conjunction with the elliptic reconstruction, results into the derivation of -norm a posteriori estimators via a standard energy argument. The efficiency of the a posteriori estimators is assessed numerically. The reliability bounds are used within two variants of a space-time adaptive algorithm. The adaptive algorithm is able to achieve the same error reduction with far fewer degrees of freedom compared to uniform meshes, thereby highlighting the relevance of the derived a posteriori estimates in practical computations. The simple model problem (1) - (3) appears to be sufficient in highlighting some of the challenges in the error estimation and adaptivity of finite element methods for more complex fourth order parabolic problems. It appears that the derived a posteriori bounds and the respective adaptive algorithms can be modified in a straightforward fashion to include the original dG method of Baker [4] and -interior penalty methods [18, 10].

The remaining of this work is organised as follows. In Section 2, notation is introduced and some standard results needed in the subsequent analysis are recalled. The discontinuous Galerkin (dG) method for the biharmonic problem, along with the derivation of posteriori error bounds for the dG approximation of the biharmonic problem in -norm are derived in Section 3. The respective fully discrete scheme for the parabolic model problem (1) - (3) is given in Section 4, while Section 5 contains the derivation of residual type a posteriori estimates of errors in and norms for the fully discrete scheme. The efficiency and reliability of the a posteriori estimators is tested on a range of uniform meshes in Section 6. The adaptive algorithm utilizing the a posteriori estimates in a series of numerical experiments are also presented in Section 6. Some concluding remarks regarding the results and possible extensions are given in Section 7.

## 2 Notation and preliminaries

The standard Hilbertian Lebesgue space is denoted by , for a domain , , with corresponding inner product and norm ; when , we shall drop the subscript writing and , respectively. We also denote by , the standard Hilbertian Sobolev space of index of real-valued functions defined on , along with the corresponding norm and seminorm and , respectively. For , we also define the spaces , consisting of all measurable functions , for which

 ∥v∥Lp(0,T;Hs(ω)):=(∫T0∥v(t)∥ps,ωdt)1/p<+∞,for1≤p<+∞,∥v∥L∞(0,T;X):=esssup0≤t≤T∥v(t)∥s,ω<+∞,forp=+∞. (4)

Let be a subdivision of into disjoint elements . The subdivision is assumed to be shape-regular (see, e.g., p.124 in [15]) and is constructed via smooth mappings with uniformly bounded Jacobian throughout the mesh family considered, where is the reference element. The above mappings are assumed to be constructed so as to ensure and that the elemental edges are straight segments (i.e., lines or planes). Note that we also use the expression edge to mean side when .

The broken Laplacian, , is defined element-wise by for all .

For a nonnegative integer , we denote by , the set of all polynomials of total degree at most , if is the reference simplex, or of degree at most in each variable, if is the reference hypercube. We consider the finite element space

 Sr:={v∈L2(Ω):v|κ∘Fκ∈Pr(^κ),κ∈T}. (5)

By we denote the union of all -dimensional element edges associated with the subdivision , including the boundary. Further, we decompose into two disjoint subsets , where .

For two (generic) elements sharing an edge , we define the outward normal unit vectors and on corresponding to and , respectively. For functions and , that may be discontinuous across , we define the following quantities. For , , , and , we set

 {v}:=12(v++v−), {\boldmathq}:=12(\boldmathq++\boldmathq−),[[v]]:=q+\boldmathn++q−\boldmathn−, [% \boldmathq]:=\boldmathq+⋅\boldmathn++\boldmathq−⋅\boldmathn−;

if , these definitions are modified to , , , . With the above definitions, it is easy to verify the identity

 Missing dimension or its units for \hskip (6)

with denoting the outward normal unit vector on , corresponding to .

We define the element size , where is the -dimensional Lebesgue measure; we collect the element sizes into the into the element-wise constant function , with , and on . Also, for two (generic) elements , sharing an edge , we define .

As we shall be dealing with mesh adaptive algorithms below, we assume that all sequences of meshes considered in this work are locally quasi-uniform, i.e., there exists constant , independent of , such that, for any pair of elements and in which share an edge,

 c−1≤hκ+/hκ−≤c. (7)

Finally we recall a series of some (standard) results used throughout this work; their proofs can be found, e.g., in [15, 17, 7, 9, 11].

###### Lemma 2.1 (approximation property)

Let and be a subdivision of , . Then there exists a constant , independent of , such that for any and , there exists , with and

 |u−p|j,κ≤Capp hm−jκ |u|m,κ ,0≤j≤m. (8)
###### Lemma 2.2 (inverse estimate)

There exists a constant , independent of , such that

 |p|j,κ≤Cinv hi−jκ |p|i,κ ,0≤i≤j≤2, (9)

for all , with .

###### Lemma 2.3 (trace inequality)

For every , with , there exists a constant independent of such that

 ||u||20,∂κ≤Ctr(h−1κ||u||20,κ+hκ|u|21,κ). (10)
###### Lemma 2.4 (Poincaré-Friedrichs inequality [11])

There exists a constant , independent of , such that for any , with for all , we have

 ||u||20,Ω+|u|21,Ω≤Cpf(|u|22,Ω+||\boldmathh−3/2[[u]]||20,Γ+||% \boldmathh−1/2[∇u]||20,Γ). (11)

## 3 Discontinuous Galerkin method for the biharmonic problem

We consider the biharmonic equation

 Δ2~u=ϕin Ω, (12)

with homogeneous essential boundary conditions

 ~u=0,∇~u⋅\boldmathn=0on ∂Ω, (13)

where denotes the unit outward normal vector to and . Then the regularity of the problem implies that [26].

Upon defining the lifting operator by

 ∫ΩL(ν)ψdx=∫Γ([[ν]]⋅{∇ψ}−{ψ}[∇ν])ds∀ψ∈Sr, (14)

the (symmetric) interior penalty discontinuous Galerkin (dG) method for (12), (13) is given by:

 find ~uh∈Sr such thatB(~uh,vh)=l(vh)∀vh∈Sr, (15)

where the bilinear form and the linear form are given by

 B(w,v):= ∫Ω(ΔhwΔhv+L(w)Δhv+ΔhwL(v))dx+Bp(w,v) (16)

with

 Bp(w,v):=∫Γ(σ[[w]]⋅[[v]]+ξ[∇w][∇v])ds,

and

 l(v):=∫Ωϕvdx, (17)

respectively, for . The piecewise constant discontinuity penalization parameters are given by

 σ|e=σ0(\boldmathh|e)−3,ξ|e=ξ0(\boldmathh|e)−1, (18)

respectively, where and . To guarantee the stability of the IPDG method defined in (15), and must be selected sufficiently large.

Note that this formulation is inconsistent for trial and test functions belonging to the solution space . However, when , in view of (14), (16) gives

 B(w,v)= ∫ΩΔhwΔhvdx+∫Γ({∇Δw}⋅[[v]]+{∇Δv}⋅[[w]] (19) −{Δw}[∇v]−{Δv}[∇w]+σ[[w]]⋅[[v]]+ξ[∇w][∇v])ds;

therefore, (15) coincides with the symmetric version interior penalty method presented in [37]. For the bilinear form in (16) we have the continuity and coercivity with respect to the energy norm on defined by

 |||w|||=(||Δhw||2Ω+||√σ[[w]]||2Γ+||√ξ[∇w]||2Γ)12. (20)
###### Lemma 3.1 ([23])

For sufficiently large and there exist positive constants and , depending only on the mesh parameters such that

 |B(u,v)|≤C\it{cont}|||u||||||v|||∀u,v∈Sand (21)
 B(u,u)≥C\it{coer}|||u|||2∀u∈S. (22)

An a posteriori bound for the energy norm error of the dG method (15) for (12) – (13) has been considered in [24]. Now, we shall present an a posteriori bound for the -norm error (cf. [35] for a corresponding result for the second order problem).

###### Theorem 3.2 (L2-a posteriori bounds for the elliptic problem)

Let be the solution of (12)–(13), be dG approximation (15) associated with the mesh . Then, there exists a positive constant , independent of , , and , such that

 ∥~u−~uh∥≤E(T,~uh,ϕ), (23)

where

 E(T,~uh,ϕ) :=C(???)(∥% \boldmathh4−λ/2(ϕ−Δ2h~uh)∥2+∥\boldmathh(7−λ)/2[∇Δ~uh]∥2Γ\rm int+∥\boldmathh(5−λ)/2[[Δ~uh]]∥2Γ\rm int (24) :=+∑e∈Γ(h3−λe(1+ξ20)∥[∇~uh]∥2+h1−λe(1+σ20)∥[[~uh]]∥2)

and .

###### Proof.

The dual problem

 Δ2z=~u−~uh=:~ein Ω, (25)

with homogeneous essential boundary conditions clearly satisfies and, therefore, the following regularity estimate holds

 ∥z∥4,Ω≤Creg∥~e∥. (26)

Using (25), integrating by parts twice, applying (6) and (14) as well as the regularity of the dual solution, , we have

 ∥~e∥2= ∑κ∈T∫κΔ2z~edx=∫ΩΔhzΔh~edx−∫Γ[∇~e]{Δz}ds+∫Γ[[~e]]⋅{∇Δz}ds. (27)

By using the fact that is a weak solution and integrating the term involving by parts, we arrive at,

 ∥~e∥2= B(~u,z)−∫ΩΔhzΔh~uhdx+∫Γ[∇~uh]{Δz}ds−∫Γ[[~uh]]⋅{∇Δz}ds (28) = l(z)−∫ΩzΔ2~uhdx+∫Γ\rm int{z}[∇Δh~uh]ds−∫Γ\rm int{∇z}⋅[[Δh~uh]]ds +∫Γ[∇~uh]{Δz}ds−∫Γ[[~uh]]⋅{∇Δz}ds.

Using the standard orthogonal -projection, , of , we can derive the following identity by integrating by parts and using (6) and (14) as follows,

 0=l(−Πz)−B(~uh,−Πz) (29) = ∑κ∈T∫κ((ϕ−Δ2~uh)(−Πz)−L(~uh)Δh(−Πz))dx+∫Γ\rm int% {−Πz}[∇Δh~uh]ds −∫Γ\rm int{∇(−Πz)}⋅[[Δh~uh]]ds−∫Γ(σ[[~uh]]⋅[[−Πz]]+ξ[∇~uh][∇(−Πz)])ds.

Using (14) in (29) and combining (28) and (29), we get

 ∥~e∥2= ∥~e∥2+l(−Πz)−B(~uh,−Πz) (30) = ∫Ω(ϕ−Δ2h~uh)(z−Πz)dx+∫Γ\rm int({z−Πz}[∇Δ~uh]ds−[[Δ~uh]]⋅{∇(z−Πz)})ds −∫Γ[[~uh]]⋅({∇Δ(z−Πz)}+σ0\boldmathh−3[[z−Πz]])ds +∫Γ[∇(~uh)]({Δ(z−Πz)}+ξ0\boldmathh−1[∇(z−Πz)])ds.

The assertion then follows by applying Young’s inequality, the trace inequality (10) where appropriate, the approximation property (8) and the regularity of the dual problem on each of the terms on the right hand side of (30). ∎

###### Remark 3.3

If a smooth subspace of the finite element space exists, such as Argyris elements in two dimensions, or corresponding constructions in three dimensions, it is possible to establish an a posteriori bound without dependence on penalty parameters; indeed, these terms would vanish from (30) if the projection, , could be defined onto the smooth subspace of .

###### Remark 3.4

It is interesting to note that the a posteriori error bound of (23) reflects the suboptimal -norm error convergence of the dG method when quadratic polynomials are applied. Similar behaviour is observed theoretically and numerically in [23] and in [37] in the context of the a priori error analysis of the same method.

## 4 DG method for the parabolic problem

Throughout the remaining of this work, we shall denote by the weak solution of the problem (1)–(3) in variational form: find such that

 ⟨ut,ϕ⟩+B(u,ϕ) =⟨f,ϕ⟩∀ϕ∈H20(Ω), (31) u =u0∈L2(Ω) in Ω×{0}.

We consider a subdivision of the time interval to be the family of intervals ; , with , and , with local time-step . Associated with this time-subdivision, let , , be a sequence of meshes which are assumed to be compatible, in the sense that for any two consecutive meshes and , can be obtained from by locally coarsening some of its elements and then locally refining some (possibly other) elements. The finite element space corresponding to will be denoted by and the respective dG bilinear form by . The backward Euler-dG method for approximating (31) is then given by: for each , find

 Un∈Sr,n such that ⟨Un−Un−1λn,V⟩+Bn(Un,V)=⟨~fn,V⟩∀V∈Sr,n, (32)

where and for is a piecewise polynomial of degree in time -projection in time of the source function . In practice, it suffices to take to achieve a first-order-in-time convergent method. We also set , with is the orthogonal -projection operator onto the finite element space .

## 5 A posteriori bounds for the parabolic problem

We shall derive a posteriori error bounds for the backward Euler-dG method (32) measured in - and -norms. To this end, we shall employ an energy argument (with carefully defined test functions) in conjunction with the elliptic reconstruction technique [32, 29, 25].

We begin by extending the sequence of numerical solutions into a continuous piecewise linear function of time

 U(0)=Π0u0 and U(t):=t−tn−1λnUn+tn−tλnUn−1 (33)

for and . Further, the discrete elliptic operator is defined by

 for ϕ∈Sr,n,⟨Anϕ,χ⟩=Bn(ϕ,χ)∀χ∈Sr,n. (34)

We now give definitions of the estimators involved in the estimation of the parabolic part of the error. Estimators at time step are denoted by subscript and subscript will be used for the cases of - and -bounds presented below, respectively.

###### Definition 5.1 (estimators for the parabolic error)

We define the coarsening or mesh-change estimators by

 γ∞,n:=1λn∥(I−Πn)Un−1∥2,γ2,n:=∥(I−Πn)Un−1∥2+n−1∑i=1∥(Πi−Πi−1)Ui−1∥2, (35)

the time-error evolution estimators by

 η∞,n:=∥gn−gn−1∥2λn,η2,n:=∥gn−gn−1∥2λ2n+n−1∑i=1λ2i∥gi−gi−1∥2, (36)

; the data approximation error in time estimators by

 β∞,n:=∫tntn−1∥~fn−f∥2dt,β2,n:=λn∫tntn−1∥~fn−f∥2dt, (37)

and an additional space estimator given by

 ~η∞,n:=E(^Tn,Un−Un−1,gn−gn−1)2, (38)

where finest common coarsening of and for each .

Using the notation above, we are ready to state the main result.

###### Theorem 5.2 (a posteriori bound)

Let be the solution of (31), be the approximation obtained by the dG method (32) and defined by (33). Then there exist positive constants and , independent of , , and , for any such that

 ∥e∥L∞(0,T;L2(Ω)) ≤ C∞(∥e(0)∥+(N∑n=1(γ∞,n+η∞,n+β∞,n)λn)12 (39) C∞(+(N∑n=1~η∞,n)12+max0≤n≤N{E(Tn,Un,gn)}) ∥e∥L2(0,T;L2(Ω)) ≤ C2(∥e(0)∥+(N∑n=1(γ2,n+η2,n+β2,n)λn)12 (40) C2(∥e(0)∥+N∑n=1(η1,2,n+(N∑n=1E(Tn,Un,gn)2λn)12).

The proof of this theorem will be the content of the remaining of this section, split into a number of intermediate results.

We begin by defining the elliptic reconstruction , of to be the solution of the elliptic problem

 Bn(ωn,v)=⟨gn,v⟩∀v∈H20(Ω) (41)

where, as above, . We note that under the assumptions on the domain , we also have . We also extend the elliptic reconstruction into a continuous piecewise linear-in-time function

 ω(t):=t−tn−1λnωn+tn−tλnωn−1 (42)

for and . Finally, we introduce the error decomposition

 e:=U−u=ρ−ϵ, where ϵ:=ω−U, and ρ:=ω−u, (43)

where and are understood as the parabolic and elliptic error, respectively, and we set .

###### Lemma 5.3 (Error identity)

For all , , we have

 ⟨ρt,v⟩+Bn(ρ,v)=⟨ϵt,v⟩+⟨(I−Πn)Ut,v⟩+t−tnλn⟨gn−gn−1,v⟩+⟨~fn−f,v⟩, (44)

for any , with denoting the identity mapping.

###### Proof.

Firstly, from (41) and (34) we have

 Bn(ωn,v)−⟨~fn,v⟩=Bn(Un,Πnv)−⟨Πn~fn,Πnv⟩. (45)

Also, using the method (32) and the definition of the -projection we deduce

 ⟨Ut,v⟩=⟨(I−Πn)Ut,v⟩+⟨Ut,Πnv⟩=⟨(I−Πn)Ut,v⟩−(Bn(Un,Πnv)−⟨Πn~fn,Πnv⟩). (46)

For the elliptic reconstruction error we also have,

 Bn(ω−ωn,v)=t−tnλn⟨gn−gn−1,v⟩. (47)

Lastly, for the terms on the left hand side of (44), we compute

 ⟨et,v⟩+Bn(ρ,v)=⟨Ut,v⟩+Bn(ω,v)−(⟨ut,v⟩+Bn(u,v))=⟨Ut,v⟩+Bn(ω,v)−⟨f,v⟩. (48)

Using (45), (45), and (47) in (48), along with the identity completes the proof. ∎

The a posteriori bounds (39) and (40) will be derived by selecting special test functions in the energy identity (44) above, along with estimation of the terms on the right-hand side of (44). More specifically, we consider the following two test functions, for the case, and

 ¯v(t,⋅):=∫Ttρ(s,⋅)ds,t∈[0,T], (49)

for the case; this choice is motivated by Baker [3], who used a similar construction for the proof of a priori bounds for the second order wave problem. The latter test function has most notably the following properties:

 ¯v∈H4(Ω)∩H20(Ω) as ρ∈H4(Ω)∩H20(Ω)a.e. in [0,T], (50) ¯v(T,⋅)=0=Δ¯v(T,⋅),∇¯v(T,⋅)=0,and ¯vt(t,⋅)=−ρ(t,⋅),a.e. in [0,T].

Next, we consider two auxiliary functions which are needed in the consequent proofs. More specifically, on each interval , for , we define

 ~G(t):=(I−Πn)U+ψn, with ψn:=−(I−Πn)Un−1+ψn−1,ψ0:=0, (51)

and

 Missing or unrecognized delimiter for \left (52)

We note that, for each , we then have , ,

 ~Gt(t):=(I−Πn)Ut, and Gt(t):=t−tnλn(gn−gn−1). (53)

The following estimates will be used in the proof of Theorem 5.2.

###### Lemma 5.4

Let . Then, we have

 ∫τ0⟨~Gt,ρ⟩dt ≤ N∑n=1∥(Πn−I)Un−1∥max0≤t≤T∥ρ∥ (54) ∫τ0⟨Gt,ρ⟩dt ≤ N∑n=1λn∥g