NORM PRECONDITIONERS FOR DISCONTINUOUS GALERKIN HP-FINITE ELEMENT

. We consider a norm-preconditioning approach for the solution of discontinuous Galerkin ﬁnite element discretizations of second order PDE with non-negative characteristic form. In particular, we perform an analysis for the general case of discontinuous hp -ﬁnite element discretizations. Our solution method is a norm-preconditioned three-term GMRES routine. We ﬁnd that for symmetric positive-deﬁnite diﬀusivity tensors the convergence of our solver is independent of discretization, while for the semideﬁnite case both theory and experiment indicate dependence on both h and p . Numerical results are included to illustrate performance on several test cases.

1. Introduction. Recent years have seen an increasing interest in a class of non-conforming finite element approximations of elliptic boundary-value problems with hyperbolic nature, usually referred to as discontinuous Galerkin finite element methods. Various families of discontinuous Galerkin finite element methods (DGFEMs) have been proposed, particularly for the numerical solution of convectiondiffusion problems, due to the many attractive properties they exhibit. In particular, DGFEMs admit good stability properties, they offer flexibility in the mesh design (irregular meshes are admissible) and in the imposition of boundary conditions (Dirichlet boundary conditions are weakly imposed), and they are increasingly popular in the context of hp-adaptive algorithms. The increase in popularity for DGFEMs has created a corresponding demand for developing linear solvers.
Existing approaches to solving systems arising in DGFEMs include domain decomposition, either non-overlapping [12], [13] or overlapping [23] and multigrid [19], [5]. Another favoured approach consists in reformulating the problem as a system of PDE which is then solved using a mixed finite element method for which block-preconditioners can be devised [29], [24], [22].
Remarkably, preconditioned Krylov methods feature rarely and have been mostly employed for the case of time-dependent problems, whether using time-discontinuous discretizations [10], or spacediscontinuous ones [3]. One of the reasons for this notable absence may be explained through the ill-conditioning the resulting linear systems suffer from [7] coupled with a need for preconditioner design.
In this work we extend results available for the case of standard finite element methods to the hpdiscontinuous case. More precisely, it is known that finite element matrices belong to a framework that allows the design of useful preconditioners [11], [26], as well as meaningful monitoring of convergence (e.g., see [1]). In particular, it is now known that the norms associated with the finite element analysis can be employed as preconditioners. The analysis of their performance was first carried out using the concept of norm-equivalence.
The concept of norm-equivalence was formally introduced in [11] in the context of preconditioning for standard finite element discretisation of elliptic problems. The authors of [11] realized that finite element matrices usually are provided with useful preconditioners in the form of the matrix-norm associated with the finite element spaces involved in the solution process. In particular, they were concerned with (and showed mesh-independence for) the distribution of eigenvalues and singular values of the norm-preconditioned system with a view to employing an expensive method, such as Conjugate Gradients on the Normal Equations (CGNE). However, modern methods such as GMRES, cannot be analysed in terms of the distribution of the spectrum alone. Thus, norm-equivalence had to be replaced by a stronger property, field-of-value equivalence, which can be shown to be relevant in the convergence analysis of GMRES. This is the concept we employ in order to identify a useful preconditioner as well as to derive useful analytical results regarding convergence properties of our iterative scheme.
The paper is structured as follows. We first introduce the formulation of the problem, together with notation. In section 4 we derive some (continuous) norm-equivalences which we need in order to perform the analysis of our preconditioners. Section 5 describes our iterative approach and includes theoretical convergence results which are then verified in section 6.

Model Problem.
Let Ω be a bounded open (curvilinear) polygonal domain in R d , and let Γ ∂ signify the union of its one-dimensional open edges, which are assumed to be sufficiently smooth (in a sense defined rigorously later). We consider the convection-diffusion-reaction equation where f ∈ L 2 (Ω), c ∈ L ∞ (Ω), b is a vector function whose entries are Lipschitz continuous real-valued functions onΩ, andā is the symmetric diffusion tensor whose entries are bounded, piecewise continuous real-valued functions defined onΩ, with Under this hypothesis, (2.1) is termed a partial differential equation with nonnegative characteristic form. By n we denote the unit outward normal vector to Γ ∂ . We define The sets Γ − and Γ + are referred to as inflow and outflow boundary, respectively. We can also see that If Γ 0 has positive one-dimensional Hausdorff measure, we also decompose Γ 0 into two parts Γ D , Γ N and we impose Dirichlet and Neumann boundary conditions, respectively, via where we adopt the (physically reasonable) hypothesis that b · n ≥ 0 on Γ N , whenever the latter is nonempty. Existence and uniqueness of solutions (in the weak sense) for the corresponding homogeneous problem were considered by Houston & Süli [21] (cf. also [27] and the references therein, for existence and uniqueness results for classical solutions), under the assumption that there exists a positive constant γ 0 such that 3. Discontinuous Galerkin Finite Element Method. We denote by L p (ω), 1 ≤ p ≤ ∞, the standard Lebesgue spaces, ω ⊂ R d with corresponding norms · L p (ω) ; the norm of L 2 (ω) will be denoted by · ω for brevity. We also denote by H s (ω), the standard Hilbertian Sobolev space of index s ≥ 0 of real-valued functions defined on ω ⊂ R d .
Let T be a subdivision of Ω into disjoint open elements κ ∈ T such that each side of κ has at most one regular hanging node. We assume that the subdivision T is constructed via mappings F κ , where F κ :κ := (−1, 1) 2 → κ is a C 1 -diffeomorphism, with non-singular Jacobian. The above mappings are assumed to be constructed so as to ensure that the union of the closures of the elements κ ∈ T forms a covering of the closure of Ω, i.e.,Ω = ∪ κ∈Tκ .
We assign to the subdivision T the broken Sobolev space of order s, H s (Ω, T ) := u ∈ L 2 (Ω) : u| κ ∈ H s (κ) for all κ ∈ T , 2 equipped with the standard broken Sobolev norm. For a nonnegative integer p, we denote by Q p (κ), the set of all tensor-product polynomials on κ of degree p in each coordinate direction. To each κ ∈ T we assign nonnegative integers p κ (the local polynomial degree), we collect the p κ and F κ into piecewise continuous functions p : Ω → R and F : Ω → R, with p| κ = p κ and F| κ = F κ , κ ∈ T , respectively, and consider the finite element space Similarly, we introduce the mesh quantities h ⊥ : ∪ κ∈T ∂κ → R and r : Ω → R where (h ⊥ )| e⊂∂κ = |κ|/|e|, for all elemental faces e ⊂ Γ, and (r ⊥ )| κ = ρ κ , κ ∈ T , respectively, where ρ κ denotes the radius of the largest inscribed circle in the element κ.
We assume that the discretisation parameters satisfy the following conditions: on each elemental face e = ∂κ∩∂κ ′ , with e having positive (d−1)-dimensional measure, and κ, κ ′ two generic neighbouring elements of the subdivision T , we have for all e ⊂ Γ, with c 1 , c 2 , c 3 positive constants. By Γ we denote the union of all one-dimensional element faces associated with the subdivision T (including the boundary). Further we decompose Γ into three disjoint subsets Γ Next, we introduce some trace operators. Let κ, κ ′ be two (generic) elements sharing a face e :=κ ∩κ ′ ⊂ Γ int . Define the outward normal unit vectors n + and n − on e corresponding to ∂κ and ∂κ ′ , respectively. For functions q : Ω → R and φ : Ω → R 2 that may be discontinuous across Γ, we define the following quantities. For q + := q| ∂κ , q − := q| ∂κ ′ and φ + := φ| ∂κ , φ − := φ| ∂κ ′ , we set Further, we divide the union of all open edges ∂κ of every element κ into two subsets where n(·) denotes the unit outward normal vector function associated with the element κ; we call these the inflow and outflow parts of ∂κ respectively. Then, for every element κ ∈ T , we denote by u + κ the trace of a function u on ∂κ taken from within the element κ (interior trace). We also define the exterior trace u − κ of u ∈ H 1 (Ω, T ) for almost all x ∈ ∂ − κ\Γ to be the interior trace u + κ ′ of u on the element(s) κ ′ that share the edges contained in ∂ − κ\Γ of the boundary of element κ. Then, the upwind jump of u across ∂ − κ\Γ is defined by We note that this definition of jump is not the same as the one above; here the sign of the jump depends on the direction of the flow, whereas [ [·] ] depends only on the element-numbering. We note that the subscript in this definition will be suppressed when no confusion is likely to occur.
The broken weak formulation of the problem (2.1), (2.3), from which the interior penalty DGFEM will emerge, reads for θ ∈ {−1, 1}, with the function σ defined by where p and h ⊥ as above, and a : Ω → R with a| κ = | √ā | 2 L ∞ (κ) , κ ∈ T where | · | 2 denotes the matrix-2-norm, and C σ positive constant. Here and in the sequel, we assume that there exists a constant α > 0, such that for all pairs of neighbouring elements κ and κ ′ having a common face e, whenever a is non-zero.
Then, the interior penalty DGFEM for the problem (2.1), (2.3) is defined as follows: We refer to the DGFEM with θ = −1 as the symmetric interior penalty DGFEM (SIPG), whereas for θ = 1 the DGFEM will be referred to as the non-symmetric interior penalty DGFEM (NIPG). This terminology stems from the fact that when b ≡ 0, the bilinear form B(·, ·) is symmetric if and only if θ = −1. Various types of error analysis for the variants of interior penalty DGFEMs can be found in [2,28,20,14,18,16]. We make some assumptions on the regularity of the solution and on the functions in the finite element space S. We assume that p κ i ≥ 1, i = 1, 2, κ ∈ T , whenever diffusion is present, in order to ensure that the matrix of the system of linear algebraic equations that arises from (3.4) is nonsingular. When the analytical solution u ∈ A, the following Galerkin orthogonality property holds: B(u − u DG , v) = 0 for every v ∈ S. If the continuity assumptions involved in the definition of A are violated, as is the case, for example, in an elliptic transmission problem, the DGFEM has to be modified accordingly.
We conclude this section by introducing the notion of inverse property for the diffusion tensorā. Definition 3.1 (inverse property). We shall say that the tensorā has the inverse property if an inverse inequality of the form holds, for all v ∈ Q p (κ), where e is a face of the element κ, and C inv is a positive constant independent ofā and the discretisation parameters. We refer to [15] for a detailed discussion on diffusion tensors satisfying the inverse property. For a discussion on how to circumvent the necessity of the inverse property assumption, by modifying the DGFEM, we refer to [16].
We are interested in studying the symmetric and the skew-symmetric parts of the bilinear form B(·, ·). For this reason, we rewrite the numerical fluxes for the convection part of the boundary as described in the following lemma.
Lemma 4.1. Using the notation above, the following identity holds: Recalling the definitions of [ [·] ] and { {·} } on the boundary Γ ∂ , along with −(b·n) = |b·n| and (b·n) = |b·n| on the inflow and outflow parts of the boundary, respectively, it is immediate that Summing (4.2) and (4.3), the result follows. 2 Splitting ideas for the numerical fluxes of DGFEM for first order hyperbolic problems are presented in [6] when the convection is in divergence form; here we modified the argument presented in that work for the case where the convection is of the form b · ∇u (cf. also Section 2.7 in [14]).
Motivated by identity (4.1), we decompose B(·, ·) into symmetric and skew-symmetric components. Lemma 4.2. The bilinear form can be decomposed into symmetric and skew-symmetric parts:

4)
and Proof. Straightforward calculation yields that with B symm (u, v) as defined in (4.4). Integration by parts of the second term in the first integral on the right-hand side of (4.5) yields The result now follows by making use of the (standard) identity (see, e.g., [6]) 2 Note that for the skew-symmetric part of the bilinear form we have B skew (u, u) = 0, for all u ∈ H 3/2+ε (Ω, T ). Lemma 4.3. Let the domain Ω, its subdivision T and the positive function σ be defined as above, with the constant C σ sufficiently large when θ = −1 (see [2,14] for details). Let alsoā satisfy the inverse property. Then, the following equivalence holds

5)
for all u ∈ S, where κ 1 and κ 2 are positive constants independent of the mesh parameters. Proof. We have Thus, when θ = 1, (4.5) holds with κ 1 = 1 = κ 2 . When θ = −1, we work as follows. Using the Cauchy-Schwarz inequality together with the inverse inequality (3.5), we obtain using Young's inequality and choosing C σ sufficiently large. Thus, the left inequality in (4.5) holds with κ 1 = 1 2 and the right inequality holds with κ 2 = 3 2 . 2 Along the same lines, we can show the continuity for the full bilinear form B(·, ·). Theorem 4.4. Let the domain Ω, its subdivision T and the positive function σ be defined as above, with the constant C σ sufficiently large when θ = −1 (see [2,14] for details). Let alsoā satisfy the inverse property. Further assume that c 0 (x) ≥ γ κ 0 > 0 for all x ∈ κ, for each κ ∈ T . Then, for u, v ∈ S and C > 0 constant, independent of the problem and the mesh parameters, we have on the elements where √ā is not positive definite. Proof. We have Using Cauchy-Schwarz inequality together with the inverse inequality (3.5), and working as in (4.6), we conclude that and similarly for the last term on the right-hand side of (4.7). Next, we bound each of the terms of B skew (·, ·) separately, for every κ ∈ T . To this end, we consider 2 cases: Case 1.ā is positive definite. We have for e ⊂ Γ int , where e =κ + ∩κ − , for two (generic) elements κ + , κ − ∈ T , we have and similarly for the skew-symmetric counterparts. If θ = 1, the remaining terms of the skew-symmetric part of the bilinear form are bounded using (4.8). Case 2.ā is positive semi-definite, but not positive definite. We have Also, for e ⊂ Γ int , where e =κ + ∩κ − , for two (generic) elements κ + , κ − ∈ T , we have where the last inequality follows by observing that r| κ = r| e⊂∂κ ≤ h ⊥ | e⊂∂κ . Similarly we can bound the rest of the terms. If θ = 1, the remaining terms of the skew-symmetric part of the bilinear form are bounded using (4.8). Summing up the resulting bounds, and using the discrete version of the Cauchy-Schwarz inequality, the result follows. 2 Next we consider the case when c 0 ≡ 0 on Ω. This includes the case of pure convection-diffusion where the wind b is an incompressible vector field.
Theorem 4.5. Let the domain Ω, its subdivision T and the positive function σ be defined as above, with the constant C σ sufficiently large when θ = −1 (see [2,14] for details). We assume that the diffusion tensorā is positive definite, and that c 0 ≡ 0 in Ω. Let also the mesh-regularity assumption hold for all elemental faces e ⊂ Γ int , for some constantC > 0 independent of the problem and the mesh parameters. Then, for u, v ∈ S, we have for C > 0 constant independent of the problem and the mesh parameters Proof. The symmetric part of the bilinear form along with the face integrals that contain the diffusion tensor can be bounded as in the proof of Theorem 4.4. We bound each of the terms of B skew (·, ·) separately, for every κ ∈ T . We have Also, for e ⊂ Γ int , where e =κ + ∩κ − , for two (generic) elements κ + , κ − ∈ T , we have and similarly for the skew-symmetric counterparts. Summing up and using the discrete version of the Cauchy-Schwarz inequality, we deduce Using the Poincaré-Friedrichs inequality for broken H 1 -functions presented in [4] (inequality (1.8) therein; following the notation of that work we have chosen Γ to be the Dirichlet boundary Γ D ), along with the mesh-regularity assumption, we deduce L ∞ (Ω) | w| 2 , which yields the result. 2 5. Norm preconditioners. We now turn to the issue of preconditioner-design by first considering an analysis tool for GMRES. In the following, we use the standard definition x 2 H = x, x H = x T Hx where H ∈ R n×n is symmetric and positive-definite and x ∈ R n . We also introduce the following definition.
Definition 5.1. Field-of-values (F OV ) equivalence Non-singular matrices A, B ∈ R n×n are said to be F OV -equivalent if there exist constants ξ 1 , ξ 2 independent of n such that for all x ∈ R n \ {0} We write We remark here that F OV -equivalence implies that ξ 1 < λ(AB −1 ) < ξ 2 , where λ(AB −1 ) denote the eigenvalues of AB −1 . This is usually a desirable property in the context of preconditioning. The above definition was introduced with the following convergence result in mind.
Lemma 5.2. If A ≈ H B the GMRES algorithm converges with respect to ·, · H in a number of iterations independent of n. Moreover, the residuals satisfy [9], [30, Thm 6.7] where ξ 1 , ξ 2 are the constants in Definition 5.1. The requirements of Definition 5.1 can be shown to be satisfied for the matrices resulting from our DG formulation. We first note that the coercivity result of Lemma 4.3 and the continuity result of Thms 4.4 have the following discrete counterparts where η 1 = κ 1 and η 2 is the continuity constant either in Thm 4.4 or in Thm 4.5. We denoted here the matrix representation of the discrete bilinear form B(·, ·) by K and by H the discrete representation of norm | ·| . Furthermore, the equivalence (4.5) has the following discrete version where K s = (K + K T )/2. The following two results taken from [26] are sufficient to ensure that H is F OV -equivalent to K and therefore represents a good preconditioning candidate. Given equivalence (5.3), one deduces that (5.2) holds also with H replaced by K s so that we get the following Proposition 5.4. Let (5.2), (5.3) hold. Then K ≈ K −1 s K s . The last two results suggest that both H and K s are candidates for optimal preconditioning inside a GMRES algorithm (cf. Lemma 5.2). However, only the latter can be employed in an efficient fashion, as we describe below.

Three-term GMRES.
It is well-known that while GMRES is one of the most robust methods available, it is not also the most efficient. In particular, the construction and storage of orthonormal Arnoldi vectors is the main hindrance in a practical context. A short (m-)term recurrence for GMRES (which entails the storage of only m vectors) is not guaranteed to exist for any given matrix. However, there is a certain class of matrices which affords a 3-term recurrence, as the following result shows. For a proof see [1]. See also [8] for related work on preconditioning with the symmetric part of a matrix.
Lemma 5.5. Let A = I + S, where S = −S T . Then the GMRES algorithm applied to matrix A is a 3-term recurrence. Given the above result, it is straightforward to choose a preconditioner for K since where S is a skew-symmetric matrix. From an implementation point of view, employing GMRES with system matrix K −1/2 s KK −1/2 s is equivalent to running GMRES in the K s -inner product and using K s as a left-preconditioner (see [26] for more details). This is our approach below. We end with the theoretical bound on GMRES convergence for the DG convection-diffusion problem.
Theorem 5.6. Let (5.2) hold. Then the residuals of the 3-term GMRES algorithm in the K s -inner product satisfy The above bound is discretization-independent only if the diffusion tensorā is positive definite. In this case, the constant η 2 is independent of both h and p, though it remains dependent on the PDE coefficients. For the case when the diffusion tensor is positive semi-definite we cannot expect optimal performance. While the constant η 1 is still innocuous, the continuity constant η 2 depends effectively on both discretization and coefficients. In particular, for uniform meshes of size h with constant degree p of the approximating polynomials, we essentially have at points where the diffusion vanishes This is a pessimistic estimate, as we shall see below.
6. Numerical Examples. We now demonstrate the validity of our theoretical bounds on two classes of problems. The first corresponds to a diffusion tensor which is positive-definite. For this class we chose to present standard convection-diffusion problems, given their ubiquitous nature in fluid modelling. We note that for non-diagonal diffusion tensors the results are expected to be similar. The second class of problems corresponds to a singular diffusion tensor. In this case, performance should deteriorate with the discretization.
6.1. The positive definite case. Our first example is that of a standard convection-diffusion operator. In this caseā = ǫI 2 is positive definite. We chose the work with two values for b: constant and circular. We solved −ǫ∆u + b · ∇u = f for (x, y) ∈ (0, 1) 2 , with b = (1, 1), subject to a Dirichlet boundary condition, which, along with the forcing function f , is chosen so that the analytical solution is This problem was considered in [20] (Example 3). The solution exhibits boundary layer behaviour along x = 1 and y = 1, and the layers become steeper as ǫ → 0. Note that the theory developed above includes this case, as here ∇ · b = 0 on Ω. We solved the problem for a range of ǫ. Discretisations for a range of uniform h (meshsize) and p (degree of polynomial approximation) were employed. The results are presented  As predicted by theory, the number of iterations is independent of discretisation parameters. The dependence of ǫ of bound (5.4) gives also a description of convergence behaviour.  12  13  7  1  10,000  36  40  29  40,000  124  117  69  5,625  18  17  12  2  22,500  61  59  60  90,000  235  231  137  10,000  39  29  23  3  40,000  112  114  100  160,000 > 300 > 300 > 300 Table 6.2 GMRES iterations for DGFEM-discretisation of the convection diffusion problem with ILU(10 −2 ) preconditioning.
For comparison purposes, we included the corresponding GMRES runs for the choice of a blackbox preconditioner such as ILU. The results are presented in Table 6.2. We can see that while the number of iterations is low for some values of the parameters, the overall convergence behaviour is quite undesirable, with iteration counts growing with both discretisation parameters. Thus, while the number of iterations appears to be decreasing with ǫ, it is exactly for this range that the discretisation parameters have to be increased in order to resolve layers. The resulting convergence behaviour becomes rapidly too costly to implement in practice. We note here that the ILU preconditioner is implemented with a standard full GMRES routine, which means that the storage increases with every iteration.
The solution exhibits boundary layer behaviour along x = 1 and the layers become steeper as ǫ → 0. We note that c 0 = −1/2∇ · b ≡ 0 on Ω. The iteration count for this case is presented in Table 6.3. We remark here that the same behaviour is expected and observed.  Table 6.3 GMRES iterations for DGFEM-discretisation of the convection diffusion problem with circular wind and with preconditioner Ks.
6.2. The positive semi-definite case. The second class of examples we consider is that for symmetric positive semidefinite diffusion tensors. We first consider the Grušin-type boundary-value problem −u xx − 16x 6 u yy + (1 − y)u y = f in Ω ≡ (−1, 1) 2 , u = 0 on ∂Ω, (6.1) with f is chosen so that the analytical solution is The diffusion tensor is positive-definite everywhere except at x = 0. Note that the analytical solution does not belong globally to H 1 (Ω) due to a singularity of the gradient at the origin; nevertheless, it is analytic inΩ\{(0, 0)}. The iteration count for this case is presented in Table 6.4. As expected, the results depend on both discretization parameters. However, while the dependence of the number of iterations is roughly like h −1 , the p−dependence is better than the predicted p 2 . p n its 2,304  22  1  10,000  39  40,000  67   5,184  42  2  22,500  74  90,000  121   9,216  62  3 40,000 107 160,000 183 Table 6.4 GMRES iterations for DGFEM-discretisation of the degenerate convection diffusion problem with preconditioner Ks.
Finally, we consider the following equation on Ω = (−1, 1) 2 −x 2 u yy + u x + u = 0, for − 1 ≤ x ≤ 1, y > 0, , for − 1 ≤ x ≤ 1, y > 0; sin 1 2 π(1 + y) exp(−x), for − 1 ≤ x ≤ 1, y ≤ 0, along with an appropriate Dirichlet boundary condition. This problem is of changing-type, as there exists a second order term for y > 0, which is no longer present for y ≤ 0. Moreover, we can easily verify that its analytical solution u exhibits a discontinuity along y = 0, although the derivative of u, in the direction normal to this line of discontinuity in u, is continuous across y = 0. We test the performance of DGFEM with θ = 1 by employing various meshes. As mentioned at the end of Section 3, we have to modify the method by setting σ| e = 0 for all element edges e ⊂ (−1, 1) × {0}, where σ e denotes the discontinuity-penalisation parameter; this is done in order to avoid penalising physical discontinuities. Note that the diffusive flux ( √ā ∇u)·n is still continuous across y = 0, and thus the method still applies. The results are shown in Table 6.5. While the dependence on parameters h and p has not changed radically, the actual iteration count is much greater, highlighting the fact that coercivity is weaker in that part of the domain for which y ≤ 0.