An hp-Adaptive Iterative Linearization Discontinuous-Galerkin FEM for Quasilinear Elliptic Boundary Value Problems

In this article we consider the a posteriori error analysis of hp-version discontinuous Galerkin finite element methods for the numerical solution of a second-order quasilinear elliptic boundary value problem of strongly monotone type. In particular, we employ and analyze a practical solution scheme based on exploiting a discrete Kačanov iterative linearization. The resulting a posteriori error bound explicitly takes into account the three sources of error: discretization, linearization, and linear solver errors. Numerical experiments are presented to demonstrate the practical performance of the proposed hp-adaptive refinement strategy.


Introduction
In this article, we consider the a posteriori error analysis, in a natural meshdependent energy norm, for a class of interior-penalty hp-version discontinuous Galerkin finite element methods (DGFEMs) for the numerical solution of the following quasilinear elliptic boundary value problem: Here, Ω ⊂ R 2 is a bounded polygon with a Lipschitz continuous boundary Γ , and f ∈ L 2 (Ω ), where for an open set D ⊆ Ω , we signify by L 2 (D) the space of all square integrable functions on D. Additionally, we assume that the nonlinearity µ satisfies the following assumptions: (A1) µ ∈ C 0 (Ω × [0, ∞)); (A2) there exist positive constants m µ and M µ such that m µ (t − s) ≤ µ(x x x,t)t − µ(x x x, s)s ≤ M µ (t − s), t ≥ s ≥ 0, x x x ∈Ω . We remark that, if µ satisfies (A2), there exist constants β ≥ α > 0, such that for all vectors v v v, w w w ∈ R 2 , and all x x x ∈ Ω , see [14,Lemma 2.1]. For ease of notation, in the sequel, we will simply write µ(s) instead of µ(x x x, s), thereby suppressing the explicit dependence of µ on x x x ∈ Ω . The weak formulation of (1) is to find u ∈ H 1 0 (Ω ) such that where, given w ∈ H 1 0 (Ω ), we define the bilinear form A(w; u, v) = Ω µ(|∇w|)∇u · ∇v dx x x, u, v ∈ H 1 0 (Ω ), as well as the L 2 (Ω )-inner product (v, w) L 2 (Ω ) = Ω vw dx x x, v, w ∈ L 2 (Ω ). Here, H 1 0 (Ω ) is the standard Sobolev space of first order, with zero trace along Γ , equipped with the norm v H 1 0 (Ω ) = ∇v L 2 (Ω ) , v ∈ H 1 0 (Ω ). Under the assumptions (A1)-(A2) above, it is elementary to show that the form A is strongly monotone and Lipschitz continuous in the sense that and respectively. From these properties, classical monotone operator theory implies existence and uniqueness of a solution of (3); see, e.g., [17,Theorem 3.3.23]. The exploitation of automatic adaptive hp-refinement algorithms has the potential to compute numerical solutions to partial differential equations (PDEs) in a highly efficient manner, often leading to exponential rates of convergence as the underlying finite element space is enriched; see, e.g., [9,16]. The key tool required to design such strategies is the derivation of a posteriori estimates for the Galerkin discretization errors; in recent years such bounds have been extended to the context of linearization and/or linear solver errors, cf. [1,2,4,5,7,12]. In the present article we consider the derivation of an hp-version a posteriori error bound for the DGFEM approximation of the second-order quasilinear elliptic PDE problem stated in (1). To this end, we employ the interior penalty DGFEM proposed in [8], cf. also [11], together with a discrete Kačanov iterative linearization scheme, cf. [6]. Based on the analysis undertaken in [11], together with the use of a suitable reconstruction operator, cf. [13,15], we derive a fully computable bound for the error, measured in terms of a suitable DGFEM energy norm, which separately accounts for the three main sources of error: discretization, linearization, and linear solver errors. On the basis of this a posteriori bound, we design and implement an hp-adaptive refinement algorithm which automatically controls each of these error contributions as the underlying finite element space is enriched. Numerical experiments highlighting the practical performance of the proposed adaptive strategy are presented.

Discrete hp-discontinuous Galerkin spaces
Let T h be a partition of Ω into disjoint open and shape-regular elements κ such that Ω = κ∈T h κ. We assume that each κ ∈ T h is an affine image of a given master element κ, which is either the open triangle {(x, y) : −1 < x < 1, −1 < y < −x}) or the open square (−1, 1) 2 in R 2 . By h κ we denote the element diameter of κ ∈ T h , and n n n κ signifies the unit outward normal vector to κ. We allow T h to be 1-irregular, i.e., each edge of any one element κ ∈ T h contains at most one hanging node (which, for simplicity, we assume to be the midpoint of the corresponding edge). In this context, we suppose that T h is regularly reducible (cf. [18,Section 7.1] and [11]), i.e., there exists a shape-regular conforming (regular) mesh T h (consisting of triangles and parallelograms) such that the closure of each element in T h is a union of closures of elements of T h , and that there exists a constant C > 0, independent of the element sizes, such that for any two elements κ ∈ T h and κ ∈ T h with κ ⊆ κ we have h κ/h κ ≤ C. Note that these assumptions imply that T h is of bounded local variation, i.e., there exists a constant ρ 1 ≥ 1, independent of the element sizes, such that ρ −1 1 ≤ h κ ♯ /h κ ♭ ≤ ρ 1 , for any pair of elements κ ♯ , κ ♭ ∈ T h which share a common edge e = ∂ κ ♯ ∩ ∂ κ ♭ . Moreover, let us consider the set E of all one-dimensional open edges of all elements κ ∈ T h . Further, we denote by E I the set of all edges e ∈ E that are contained in the open domain Ω (interior edges). Additionally, we introduce E B to be the set of boundary edges consisting of all e ∈ E that are contained in Γ .
For any integer p ∈ N 0 , we denote by P p (κ) the set of polynomials of total degree p on κ. Similarly, when κ is a quadrilateral, we also consider Q p (κ), the set of all tensor-product polynomials on κ of degree p in each coordinate direction. To each κ ∈ T h we assign a polynomial degree p κ (local approximation order). We collect the local polynomial degrees in a vector p p p = {p κ : κ ∈ T h }, and then introduce the hp-DGFEM space with S being either P or Q. We shall suppose that the polynomial degree vector p, with p κ ≥ 1 for each κ ∈ T , has bounded local variation, i.e., there exists a constant ρ 2 ≥ 1, independent of the local element sizes and p p p, such that, for any pair of neighbouring elements κ ♯ , κ ♭ ∈ T h , we have ρ −1 Evidently, since functions in V DG (T h , p p p) do not need to be continuous, we have that Π κ,p κ = Π T h ,p p p | κ , where, for κ ∈ T h , we let Π κ,p κ be the L 2 -projection onto S p κ (κ).

Nonlinear hp-DGFEM formulation
Let κ ♯ and κ ♭ be two adjacent elements of T h , and x x x an arbitrary point on the interior edge e ∈ E I given by e = (∂ κ ♯ ∩ ∂ κ ♭ ) • . Furthermore, let v and q q q be scalar-and vector-valued functions, respectively, that are sufficiently smooth inside each element κ ♯ , κ ♭ . Then, the averages of v and q q q at x x x ∈ e are given by v Similarly, the jumps of v and q q q at x x x ∈ e are given by = q q q| κ ♯ · n n n κ ♯ + q q q| κ ♭ · n n n κ ♭ , respectively. On a boundary edge e ∈ E B , we set v = v, q q q = q q q and [[v]] = vn n n, with n n n denoting the unit outward normal vector on the boundary Γ .
Furthermore, we introduce the edge functions h, p ∈ L ∞ (E ), which, for an edge e ∈ E , are given by h| e := h e and p| e = p | e , with h e denoting the length of e. In addition, we define the discontinuity penalisation function σ ∈ L ∞ (E ) given by σ = With this notation, following [8], we introduce the interior penalty DGFEM discretization of (3) by: where, for given w ∈ V DG (T h , p p p), we define the DGFEM bilinear form where θ ∈ [−1, 1] is a method parameter. Referring to [8, Theorem 2.5], provided that γ ≥ 1 is chosen sufficiently large (independent of the local element sizes and of the polynomial degree distribution), the existence and uniqueness of the DGFEM solution u DG ∈ V DG (T h , p p p) satisfying (5) is guaranteed.

Assumption 1
In the sequel, we suppose that there exists a computable a posteriori error estimate of the form u − u DG DG ≤ η(u DG , f ), where u ∈ H 1 0 (Ω ) is the solution of (1), and u DG is its hp-DGFEM approximation defined in (5).

Remark 1.
In the article [11,Theorem 3.2] it has been proved that such a bound does indeed exist. More precisely, we have that where, the local error indicators η κ , κ ∈ T h , are defined by and e is a data oscillation term. For κ ∈ T h and e ∈ E I , we have O 0,e , where I denotes a generic identity operator. Here, we write p p p − 1 := {p κ − 1} κ∈T h . Additionally, we denote by Π E ,p p p−1 | e the L 2 -projector onto Pp e −1 (e), where we let p e = max{p κ ♯ , p κ ♭ }, with κ ♯ , κ ♭ ∈ T h , e = ∂ κ ♯ ∩ ∂ κ ♭ . Moreover, C > 0 in (6) is a constant that is independent of the local element sizes, the polynomial degree vector p p p, and the parameters γ and θ .

Iterative DGFEM
In order to provide a practical solution scheme for the nonlinear hp-DGFEM system (5) we propose a linearization approach based on a discrete Kačanov fixed point iteration, see, e.g., [6]. To this end, we begin by selecting an initial guess u 0 DG ∈ V DG (T h , p p p). Then, for n ≥ 1, given u n−1 DG ∈ V DG (T h , p p p), we solve the linear hp-DGFEM formulation, defined by for u n DG ∈ V DG (T h , p p p). We emphasize that, in actual computations, the linear system (8) may be solved by an iterative algorithm, thereby generating an approximate numerical solution u n DG ∈ V DG (T h , p p p), with u n DG ≈ u n DG . This means that, in practice, instead of computing the sequence {u n DG } n≥0 obtained from the iteration (8), an inexact sequence { u n DG } n≥0 is generated such that From a mathematical view point, this (inexact) iterative linearization DGFEM approach gives rise to three different sources of error: 1. Discretization error, which is expressed by the residual 2. Linearization error, which is given in terms of the residual ψ n DG ∈ V DG (T h , p p p): p p p). (11) We observe that, if (1) is linear, then we immediately obtain ψ n DG = 0. 3. Linear solver error, which is described by a residual λ n DG ∈ V DG (T h , p p p): Note that, if (8) is solved exactly, then we have u n−1 DG = u n−1 DG and u n DG = u n DG , and it follows that λ n DG = 0.

Remark 2.
Since V DG (T h , p p p) may not need to be continuous along element interfaces, the linearization and linear solver residuals ψ n DG and λ n DG , respectively, can be computed elementwise, i.e., in parallel, and, hence, at a low computational cost.
The aim of the analysis in the following section is to investigate the above residuals, and then to provide a computable a posteriori error estimate for the error u − u n DG DG between the solution u of (1) and u n DG ∈ V DG (T h , p p p).

A posteriori error estimation
In order to bound the residual ρ n DG in (10), we apply an elliptic reconstruction technique along the lines of the works [13,15], see also [7]. Specifically, we define an auxiliary function u n ∈ H 1 0 (Ω ) to be the unique solution of the weak formulation where ψ n DG and λ n DG are the linearization and linear solver residuals from (11) and (12), respectively. Upon adding (11) and (12), we notice that In particular, we observe that u n DG is the DGFEM approximation of u n based on employing the (nonlinear) DGFEM scheme defined in (5). In particular, we may exploit the a posteriori error estimate in Assumption 1 to infer the computable bound We now turn to bounding the elliptic reconstruction error u − u n ∈ H 1 0 (Ω ); to this end, we first observe that u − u n DG = u − u n H 1 0 (Ω ) . Then, employing the strong monotonicity property (4), and recalling the weak formulation (3), we obtain Employing the Cauchy-Schwarz inequality, together with the Poincaré-Friedrichs inequality, v L 2 (Ω ) ≤ C PF ∇v L 2 (Ω ) for all v ∈ H 1 0 (Ω ), where C PF > 0 is a constant, we deduce that where the linearization and linear solver residuals are given, respectively, by Ψ n Summarizing the above analysis leads to the following result. (9), for n ≥ 1, the following a posteriori error bound holds:

Theorem 1. Suppose that Assumption 1 is satisfied. Then, given a sequence of (possibly inexact) DGFEM approximations
Here, u is the analytical solution of (1), ψ n DG and λ n DG are the residuals defined in (11) and (12), respectively, and α > 0 is the constant occurring in (2) and (4).
Proof. The result follows immediately upon application of the triangle inequality, i.e., u − u n DG DG ≤ u − u n DG + u n − u n DG DG , and inserting the bounds (13) and (14). Remark 3. We note that the above analysis naturally applies to other finite element schemes, provided that Assumption 1 is satisfied.

Adaptive iterative hp-DGFEM procedure
In this section we introduce an automatic hp-refinement algorithm which ensures that each of the three components of the error, namely discretization, linearization, and linear solver, are controlled in a suitable fashion. To this end, we propose the following strategy, cf. [12].

Algorithm 2 Given a (coarse) starting mesh T h , with an associated (low-order)
polynomial degree distribution p p p, and an initial guess u 0 ∈ V DG (T h , p p p). Set n ← 1.  (15) is fulfilled then the space V DG (T h , p p p) is adaptively hp-refined based on first marking elements for refinement according to the size of the local element indicators η κ , cf. (7). To this end, we exploit the maximal strategy whereby elements are marked for refinement which satisfy the condition η κ > 1 /3 max κ∈T h η κ . Secondly, once an element κ ∈ T h has been marked for refinement, we undertake either local mesh subdivision or local polynomial enrichment based on employing the hp-refinement criterion developed within the article [10]. Finally, when (15) is not fulfilled, rather than determining which source of error, i.e., the (computable) quantities Ψ n DG or Λ n DG from (11) and (12), respectively, is dominant, we choose to always undertake a further linearization step, and hence a further linear solver step is also computed, since this ensures that the most up to date approximation u n DG is employed at all times.

Application to Quasilinear Elliptic PDEs
In this section we present numerical experiments to highlight the performance of the proposed iterative hp-refinement procedure outlined in Algorithm 2. To this end, we set the interior penalty parameter constant γ to 10 and the steering parameter ϒ to 1 /4. The solution of the resulting set of linear equations is computed using an ILU(0) preconditioned GMRES algorithm. For the first numerical experiment, we let Ω = (0, 1) 2 and define the nonlinear coefficient as µ(|∇u|) = 2 + (1 + |∇u|) −1 . The right-hand forcing function f is selected so that the analytical solution to (1) is given by u(x, y) = x(1 − x)y(1 − y)(1 − 2y)e −20(2x−1) 2 . In Figure 1 we present a comparison of the actual error measured in terms of the energy norm versus the square root of the number of degrees of freedom in V DG (T h , p p p). From Figure 1(a) we clearly observe that exponential convergence of the proposed hp-refinement strategy as V DG (T h , p p p) in enriched. Furthermore, in Figure 1(b) we plot the individual residual error indicators; for this smooth problem, we observe that the discretization indicator (denoted as η n in the figure) is always dominant, while the linearization and linear solver residuals (denoted as Ψ n and λ n , respectively) are roughly of the same magnitude.
Secondly, we let Ω denote the L-shaped domain (−1, 1) 2 \ [0, 1) × (−1, 0] ⊂ R 2 and select µ(|∇u|) = 1 + exp(−|∇u| 2 ). By writing (r, ϕ) to denote the system of polar coordinates, we choose the forcing function f and an inhomogeneous boundary condition such that the analytical solution to (1) is u = r 2/3 sin ( 2 /3ϕ), cf. [3]. In  ergy norm versus the third root of the number of degrees of freedom in V DG (T h , p p p); as before we again attain exponential convergence of the proposed hp-refinement strategy as V DG (T h , p p p) is adaptively refined, though convergence of the a posteriori error estimator is no longer monotonic. Indeed, from Figure 2(b), we observe that once an hp-mesh refinement has been undertaken, then several linearization/solver steps may be required to ensure that the numerical solution has been computed to a sufficient accuracy before future refinements may be undertaken.

Conclusions
In this article we have derived a computable hp-version a posteriori error bound for the DGFEM approximation of a second-order quasilinear elliptic PDE problem, whereby a discrete Kačanov iterative linearization scheme is employed. The resulting computable upper bound directly takes into account discretization error, as well as the errors stemming from linearization and the underlying linear solver. Numer-ical experiments highlighting the performance of this bound within an automatic hp-refinement algorithm are presented.