Ani Miraçi

Jan Papež

Martin Vohralík

A multilevel algebraic error estimator and the corresponding iterative solver with p-robust behavior *

Keywords: finite element method, stable decomposition, multilevel method, Schwarz method, a posteriori estimate, p-robustness

In this work, we consider conforming finite element discretizations of arbitrary polynomial degree p ≥ 1 of the Poisson problem. We propose a multilevel a posteriori estimator of the algebraic error. We prove that this estimator is reliable and efficient (represents a two-sided bound of the error), with a constant independent of the degree p. We next design a multilevel iterative algebraic solver from our estimator and show that this solver contracts the algebraic error on each iteration by a factor bounded independently of p. Actually, we show that these two results are equivalent. The p-robustness results rely on the work of Schöberl et al. [IMA J. Numer. Anal., 28 (2008), pp. 1-24] for one given mesh. We combine this with the design of an algebraic residual lifting constructed over a hierarchy of nested unstructured, possibly highly graded, simplicial meshes. The lifting includes a global coarse-level solve with the lowest polynomial degree one together with local contributions from the subsequent mesh levels. These contributions, of the highest polynomial degree p on the finest mesh, are given as solutions of mutually independent local Dirichlet problems posed over overlapping patches of elements around vertices. The construction of this lifting can be seen as one geometric V-cycle multigrid step with zero pre-and one postsmoothing by (damped) additive Schwarz (block Jacobi). One particular feature of our approach is the optimal choice of the step-size generated from the algebraic residual lifting. Numerical tests are presented to illustrate the theoretical findings.

Introduction

The finite element method (FEM) is a widespread approach for discretizing problems given in the form of partial differential equations, and has been used in engineering for more than fifty years. For a thorough overview on the topic, we refer the reader to, e.g., Ciarlet [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF], Ern and Guermond [START_REF] Ern | Theory and practice of finite elements[END_REF], and Brenner and Scott [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]. Many iterative methods have been suggested to treat the linear systems arising from finite element discretizations; see e.g., Bramble et al. [START_REF] Bramble | The construction of preconditioners for elliptic problems by substructuring[END_REF] and [START_REF] Bramble | Parallel multilevel preconditioners[END_REF], Hackbusch [START_REF] Hackbusch | Multi-grid methods and applications[END_REF], Bank, Dupont, and Yserentant [START_REF] Bank | The hierarchical basis multigrid method[END_REF], Brandt, McCormick, and Ruge [START_REF] Brandt | Algebraic multigrid (AMG) for sparse matrix equations[END_REF], Oswald [START_REF] Oswald | Multilevel finite element approximation[END_REF], or Zhang [START_REF] Zhang | Multilevel Schwarz methods[END_REF], and the references therein. A systematic description of iterative solvers is given by Xu in [START_REF] Xu | Iterative methods by space decomposition and subspace correction[END_REF]. For convergence results on unstructured and graded meshes, we refer the reader to, e.g., Wu and Chen [START_REF] Wu | Uniform convergence of multigrid V-cycle on adaptively refined finite element meshes for second order elliptic problems[END_REF], Hiptmair, Wu, and Zheng [START_REF] Hiptmair | Uniform convergence of adaptive multigrid methods for elliptic problems and Maxwell's equations[END_REF], Chen, Nochetto, and Xu [START_REF] Chen | Optimal multilevel methods for graded bisection grids[END_REF], and Xu, Chen, and Nochetto [START_REF] Xu | Optimal multilevel methods for H(grad), H(curl), and H(div) systems on graded and unstructured grids[END_REF]. The convergence of these methods is typically robust with respect to the size of the mesh (h-robustness). In fact, this is one of the key advantages of multigrid methods. For the conjugate gradient method on the other hand, h-robustness is not intrinsic; this problem can be bypassed with the development of appropriate preconditioners.

If we are to consider methods of arbitrary approximation polynomial degree, an additional question arises: how does the polynomial degree p affect the performance of the method? In this regard, results for p-version FEM include Foresti et al. [START_REF] Foresti | Multilevel solution method for the p-version of finite elements[END_REF] for two-dimensional domains, Mandel [START_REF] Mandel | Two-level domain decomposition preconditioning for the p-version finite element method in three dimensions[END_REF] for three-dimensional domains, and Babuška et al. [START_REF] Babuška | Efficient preconditioning for the p-version finite element method in two dimensions[END_REF] for two-dimensional domains. For the latter, the condition number of the preconditioned system grows at most by 1 + log 2 (p), and a generalization of this work for hp-FEM is given by Ainsworth [START_REF] Ainsworth | A preconditioner based on domain decomposition for h-p finite-element approximation on quasi-uniform meshes[END_REF], where the p-dependence is still present. An early version of a polynomial-degree robust (p-robust) solver was introduced by Quarteroni and Sacchi Landriani [START_REF] Quarteroni | Domain decomposition preconditioners for the spectral collocation method[END_REF] for a specific domain configuration (decomposable into rectangles without internal cross points). A notable development on p-robustness was later made by Pavarino [START_REF] Pavarino | Additive Schwarz methods for the p-version finite element method[END_REF] for quadrilateral/hexahedral meshes, where the author introduced a p-robust additive Schwarz method. The generalization of this result for triangular/tetrahedral meshes was achieved by Schöberl et al. [START_REF] Schöberl | Additive Schwarz preconditioning for p-version triangular and tetrahedral finite elements[END_REF], once more by introducing an additive Schwarz preconditioner. More recent works were carried out based on these approaches. In Antonietti et al. [START_REF] Antonietti | A uniform additive Schwarz preconditioner for high-order discontinuous Galerkin approximations of elliptic problems[END_REF] (see also the references therein), the p-robust approach for rectangular/hexahedral meshes was used for high-order discontinuous Galerkin (DG) methods; moreover the spectral bounds of the preconditioner are also robust with respect to the method's penalization coefficient. We also mention the introduction of multilevel preconditioners yielded by block Gauss-Seidel smoothers in Kanschat [START_REF] Kanschat | Robust smoothers for high-order discontinuous Galerkin discretizations of advectiondiffusion problems[END_REF] for rectangular/hexahedral meshes and DG discretizations. Further multilevel approaches for rectangular/hexahedral meshes using overlapping or nonoverlapping Schwarz smoothers can be found in, e.g., Janssen and Kanschat [START_REF] Janssen | Adaptive multilevel methods with local smoothing for H 1 -and H curlconforming high order finite element methods[END_REF] and Lucero Lorca and Kanschat [START_REF] Lucero Lorca | Multilevel Schwarz preconditioners for singularly perturbed symmetric reaction-diffusion systems[END_REF]. For a study on more general meshes, see, e.g., Antonietti and Pennesi [START_REF] Antonietti | V -cycle multigrid algorithms for discontinuous Galerkin methods on non-nested polytopic meshes[END_REF], where a multigrid approach behaves p-robustly under the condition that the number of smoothing steps (depending on p) is chosen sufficiently large. Another notable contribution is the design of algebraic multigrid methods (AMG) via aggregation techniques; see, e.g., Notay and Napov [START_REF] Notay | A massively parallel solver for discrete Poisson-like problems[END_REF], Bastian, Blatt, and Scheichl [START_REF] Bastian | Algebraic multigrid for discontinuous Galerkin discretizations of heterogeneous elliptic problems[END_REF], and the references therein. The numerical results of the latter give a satisfactory indication of p-robustness.

An associated topic is the development of estimates of the algebraic error. In this regard, a posteriori tools have primarily been used to estimate the algebraic error for existing solvers. One particular goal is the development of reliable stopping criteria, allowing one to avoid unnecessary iterations. This is achieved with a combination of a posteriori error estimators for the discretization error. Some contributions on this matter (see also references therein) include Becker, Johnson, and Rannacher [START_REF] Becker | Adaptive error control for multigrid finite element methods[END_REF] where adaptive error control is achieved for a multigrid solver, Bornemann and Deuflhard [START_REF] Bornemann | The cascadic multigrid method for elliptic problems[END_REF], where a one-way multigrid method is presented by integrating an adaptive stopping criterion based on a posteriori tools. Further developments were made in Meidner, Rannacher, and Vihharev [START_REF] Meidner | Goal-oriented error control of the iterative solution of finite element equations[END_REF], where goal-oriented error estimates are used in the framework of the dual weighted residual (DWR) method. In Jiránek, Strakoš, and Vohralík [START_REF] Jiránek | A posteriori error estimates including algebraic error and stopping criteria for iterative solvers[END_REF] and later in Papež, Strakoš, and Vohralík [START_REF] Papež | Estimating and localizing the algebraic and total numerical errors using flux reconstructions[END_REF], upper and lower bounds for both the algebraic and total errors are computed, which allow one to derive guaranteed upper and lower bounds on the discretization error, and consecutively construct safe stopping criteria for iterative algebraic solvers. Arioli, Georgoulis, and Loghin [START_REF] Arioli | Stopping criteria for adaptive finite element solvers[END_REF] propose practical stopping criteria which guarantee that the considered inexact adaptive FEM algorithm converges for inexact solvers of Krylov subspace type. To the best of the authors' knowledge, though, dedicated proofs of efficiency of a posteriori estimators of the algebraic error have so far not been presented.

In this work, we present an a posteriori algebraic error estimator and a multilevel iterative solver associated to it. The cornerstone of their definitions lies in the multilevel construction of a residual algebraic lifting, motivated partly by the approach of Papež et al. [START_REF] Papež | Sharp algebraic and total a posteriori error bounds for h and p finite elements via a multilevel approach[END_REF]. The lifting can be seen as an approximation of the algebraic error by continuous piecewise polynomials of degree p, obtained by a V-cycle multigrid method with no presmoothing step and a single postsmoothing step. The coarse correction is given by a lowest polynomial degree (piecewise affine) function. Our smoothing is chosen in the family of damped additive Schwarz (block Jacobi) methods applied to overlapping subdomains composed of patches of elements (two options for defining the patches will be given in due time) and corresponds to local Dirichlet problems with the highest p-degree on the finest mesh. Note that additive Schwarz-type smoothing allows for a parallelizable implementation at each level of the mesh hierarchy. Once this lifting is built, the a posteriori estimator is easily derived as a natural guaranteed lower bound on the algebraic error, following [START_REF] Papež | Sharp algebraic and total a posteriori error bounds for h and p finite elements via a multilevel approach[END_REF] and the references therein. As our first main result, we show that up to a p-robust constant, the estimator is also an upper bound on the error.

Our solver is then defined as a linear iterative method. Because we have at hand the residual lifting, which approximates the algebraic error, we use it as a descent direction (the asymmetric, since no presmoothing is used, approach in defining the lifting will not be a problem for the analysis). The step-size is then chosen by a line search in the direction of the lifting. Our choice presents a resemblance to the conjugate gradient method, in that we choose the step-size that ensures the best error contraction in the energy norm at the next iteration. Other precedents of the use of optimal step-sizes include, e.g., Canuto and Quarteroni [START_REF] Canuto | Preconditioned minimal residual methods for Chebyshev spectral calculations[END_REF], and in the multigrid setting Heinrichs [START_REF] Heinrichs | Line relaxation for spectral multigrid methods[END_REF]. As our second main result, we prove that this solver contracts the error at each iteration by a p-robust constant. Actually, we also show that the p-robust efficiency of the estimator is equivalent to the p-robust convergence of the solver. All these results are defined for a general hierarchy of nested, unstructured, possibly highly refined (graded) matching simplicial meshes, and no assumption beyond u ∈ H 1 0 (Ω) is imposed on the weak solution. The work is structured as follows. In section 2, we introduce the setting in which we will be working as well as the notation employed throughout the paper. Then, we introduce our multilevel residual lifting construction in section 3, following Papež et al. [START_REF] Papež | Sharp algebraic and total a posteriori error bounds for h and p finite elements via a multilevel approach[END_REF]. In section 4, we present the a posteriori estimator on the algebraic error and the corresponding multilevel solver based on the residual lifting. Our main results are presented in the form of two theorems in section 5, together with a corollary establishing their equivalence. Another important corollary is the equivalence of the algebraic error with a computable estimator which is localized levelwise as well as patchwise. We provide numerical experiments in section 6, focusing mainly on showcasing p-robustness, in agreement with our theoretical results, and on a comparison with several existing approaches. We also introduce a weighted restricted additive Schwarz variant of our solver. The proofs of our main results are given in section 7. In particular, for the stable decomposition estimate, the p-robust result on one level introduced by Schöberl et al. [START_REF] Schöberl | Additive Schwarz preconditioning for p-version triangular and tetrahedral finite elements[END_REF] is crucial. We also rely on the multilevel stable splitting of Xu, Chen, and Nochetto [START_REF] Xu | Optimal multilevel methods for H(grad), H(curl), and H(div) systems on graded and unstructured grids[END_REF] for p = 1 to obtain acceptable estimates with respect to the number of levels. Finally, section 8 brings forth our conclusions and outlook for future work.

Setting

We will consider in this work the Poisson problem defined over Ω ⊂ R d , d ∈ {1, 2, 3}, an open bounded polytope with a Lipschitz-continuous boundary.

Model problem

Let f ∈ L 2 (Ω) be the source term. We consider the following problem: find u : Ω → R such that

-∆u = f in Ω, u = 0 on ∂Ω. (2.1)
In the weak formulation, we search for u ∈ H 1 0 (Ω) such that

(∇u, ∇v) = (f, v) ∀v ∈ H 1 0 (Ω), (2.2)
where (•, •) is the L 2 (Ω) or [L 2 (Ω)] d scalar product. The existence and uniqueness of the solution of (2.2) follows from the Riesz representation theorem.

Finite element discretization

Let T J be a given simplicial mesh of Ω. Fixing an integer p ≥ 1, we introduce the finite element space of continuous piecewise p-degree polynomials

V p J := P p (T J ) ∩ H 1 0 (Ω), (2.3)
where

P p (T J ) := {v J ∈ L 2 (Ω), v J | K ∈ P p (K) ∀K ∈ T J }. We set N J := dim(V p J ). The discrete problem consists in finding u J ∈ V p J such that (∇u J , ∇v J ) = (f, v J ) ∀v J ∈ V p J .
(2.4)

Algebraic system, approximate solution, and algebraic residual

If one introduces ψ l J , 1 ≤ l ≤ N J , a basis of V p J , then problem (2.4) is equivalent to solving a system of linear algebraic equations. Denoting by (A J ) lm := (∇ψ m J , ∇ψ l J ) the symmetric, positive definite (stiffness) matrix and by (F J ) l := (f, ψ l J ) the right-hand side (load) vector, one obtains

u J = N J m=1 (U J ) m ψ m J , where U J ∈ R N J is the solution of A J U J = F J .
For any approximation U i J ∈ R N J of U J given by an arbitrary algebraic solver at iteration step i ≥ 0, the associated continuous piecewise polynomial of degree p is u

i J = N J m=1 (U i J ) m ψ m J ∈ V p J .
The associated algebraic residual vector is given by R i J := F J -A J U i J . Note, however, that R i J depends on the choice of the basis functions ψ l J , 1 ≤ l ≤ N J . To avoid this dependence, we work instead with a residual functional on V p J given by

v J → (f, v J ) -(∇u i J , ∇v J ) ∈ R, v J ∈ V p J . (2.5)
We emphasize that the forthcoming results are independent of the choice of the basis.

A hierarchy of meshes

We consider a hierarchy of nested matching simplicial meshes {T j } 0≤j≤J , J ≥ 1, where T J was introduced in section 2.2, and where T j is a refinement of T j-1 , 1 ≤ j ≤ J. For any element K on a given mesh, we denote h K := diam(K) and by V K the set of its vertices. We also denote h j := max K∈Tj h K for 0 ≤ j ≤ J.

Hereafter, we shall always assume that our meshes are shape regular.

Assumption 2.1 (Shape regularity).

There exists κ T > 0 such that

max K∈Tj h K ρ K ≤ κ T for all 0 ≤ j ≤ J, (2.6)
where ρ K denotes the diameter of the largest ball inscribed in K.

Additionally to the above assumption, we will treat below two specific cases. In the first one, we suppose quasi-uniformity of the meshes in the hierarchy and that the strength of refinement is bounded. In the second case, we suppose that the meshes are generated by a series of bisections, e.g., the newest vertex bisection; cf. Sewell [START_REF] Sewell | Automatic generation of triangulations for piecewise polynomial approximation[END_REF].

A hierarchy of quasi-uniform meshes

We assume quasi-uniformity and that the hierarchy of meshes is such that the size of each parent element is comparable to the size of each of its children. Assumption 2.2 (Maximum refinement strength and mesh quasi-uniformity). There exists 0 < C ref ≤ 1, a fixed positive real number such that for any j ∈ {1, . . . , J}, ∀K ∈ T j-1 , and for any K * ∈ T j such that K * ⊂ K, there holds

C ref h K ≤ h K * ≤ h K .
(

2.7)

There further exists C qu , a fixed positive real number such that for any j ∈ {0, . . . , J} and ∀K ∈ T j , there holds

C qu h j ≤ h K ≤ h j .
(2.8)

A hierarchy of graded bisection meshes

In the case of graded mesh hierarchies obtained by bisection, one refinement of an edge of T j-1 , for j ∈ {1, . . . , J}, gives a new finer mesh T j . We denote by B j ⊂ V j the set consisting of the new vertex obtained after the bisection together with its two neighbors on the refinement edge, cf. Figure 1 for an illustration when d = 2. We denote by h Bj the maximal diameter of elements having a vertex in B j . This setting is described by the following.

Assumption 2.3 (Local quasi-uniformity of bisection-generated meshes). T 0 is a conforming quasi-uniform mesh with parameter C 0 qu . The graded conforming mesh T J is generated from T 0 by a series of bisections. There exists a fixed positive real number C loc qu such that for any j ∈ {1, . . . , J}, there holds

C loc qu h Bj ≤ h K ≤ h Bj ∀K ∈ T j such that a vertex of K belongs to B j .
(2.9)

Figure 1: Illustration of the set B j . The mesh T j-1 and its refinement T j are given by full and dotted lines, respectively.

A hierarchy of spaces

In the following, we will need to introduce a hierarchy of finite element spaces associated to the mesh hierarchy. For this purpose, let p ∈ {1, . . . , p} be a polynomial degree between 1 and p that we employ for the intermediate levels. In particular, let the following hold:

for 1 ≤ j ≤ J -1, V p j := P p (T j ) ∩ H 1 0 (Ω) (p th order spaces), (2.10a)

for j = 0, V 1 0 = P 1 (T 0 ) ∩ H 1 0 (Ω) (lowest-order space), (2.10b)
where

P p (T j ) := {v j ∈ L 2 (Ω), v j | K ∈ P p (K) ∀K ∈ T j } for 1 ≤ j ≤ J -1. Note that V 1 0 ⊂ V p 1 ⊂ • • • ⊂ V p J-1 ⊂ V p J
. Let V j be the set of vertices of the mesh T j . We denote by ψ j,a the standard hat function associated to the vertex a ∈ V j , 0 ≤ j ≤ J. This is the piecewise affine function with respect to the mesh T j that takes value 1 in the vertex a and vanishes in all other jth level vertices of V j .

Two types of patches

For the following, we define two types of patches of elements. In order to facilitate the work with both, we introduce a switching parameter s ∈ {0, 1}. First, given a vertex a ∈ V j-s , j ∈ {1, . . . , J}, we denote by T a j,s the patch formed by all elements of the mesh T j-s sharing the vertex a, i.e., T a j,s := {K ∈ T j-s , a ∈ V K }.

(2.11)

We also denote by ω a j,s the open patch subdomain corresponding to T a j,s . An illustration is given in Figure 2 (left) for "small" patches s = 0 and (right) for "large" patches s = 1. Then the associated local space V a j,s is given by V a j,s := P p (T j ) ∩ H 1 0 (ω a j,s ), j ∈ {1, . . . , J -1} and V a J,s := P p (T J ) ∩ H 1 0 (ω a J,s ).

(2.12)

Note that V a j,s are continuous piecewise polynomial spaces with respect to the mesh T j for both s = 0 and s = 1, the support being bigger in the latter case. An illustration is also given in Figure 2. j,0 associated to the "small" patch T b j,0 (left) and for the space V a j,1 associated to the "large" patch T a j,1 (right). The mesh T j-1 and its refinement T j are defined in bold and dotted lines, respectively.

Multilevel lifting of the algebraic residual

In the spirit of Papež et al. [START_REF] Papež | Sharp algebraic and total a posteriori error bounds for h and p finite elements via a multilevel approach[END_REF], we design a multilevel lifting of the algebraic residual given by (2.5). This lifting will lead to the construction of an a posteriori error estimator; it will also serve as a descent direction for the solver we introduce in the next section.

Exact algebraic residual lifting

For illustration and theoretical analysis later, we introduce the following definition. Definition 3.1 (Exact residual lifting). Let u i J ∈ V p J be arbitrary. We introduce ρi J,alg ∈ V p J as the solution of the residual problem

(∇ρ i J,alg , ∇v J ) = (f, v J ) -(∇u i J , ∇v J ) ∀v J ∈ V p J , (3.1)
so that ρi J,alg = u J -u i J .

(3.2)

Coarse solve

The first step of our construction is to solve a global lowest-order problem on the coarsest mesh. Let

u i J ∈ V p J be given. Recalling that V 1 0 = P 1 (T 0 ) ∩ H 1 0 (Ω), we define ρ i 0 ∈ V 1 0 by (∇ρ i 0 , ∇v 0 ) = (f, v 0 ) -(∇u i J , ∇v 0 ) ∀v 0 ∈ V 1 0 . (3.3)
Note that due to (3.3) and (3.1), we have

(∇ρ i 0 , ∇v 0 ) = (∇ρ i J,alg , ∇v 0 ) ∀v 0 ∈ V 1 0 , (3.4)
so that ρ i 0 is the orthogonal projection of ρi J,alg onto the coarsest space V 1 0 .

Multilevel algebraic residual lifting

Let us now introduce our hierarchical construction of the algebraic residual lifting ρ i J,alg ∈ V p J that is hopefully close to ρi J,alg . The construction relies on the use of a coarse solution of (3.3) and on local contributions arising from all the finer mesh levels. These local contributions are defined on patch subdomains ω a j,s . We denote by (•, •) ω a j,s the L 2 (ω a j,s ) or [L 2 (ω a j,s )] d scalar product. Since we consider two definitions of patches with switching parameter s ∈ {0, 1} (see section 2.6), two constructions of ρ i J,alg are implied. Definition 3.2 (Construction of the algebraic residual lifting). Let w 1 , w 2 ∈ R ∪ {∞} be damping weights satisfying the conditions

1 ≤ w 1 < 6J(d + 1) and w 2 ≥ max 1, 5J 2 (d + 1) 2 w 1 (6J(d + 1) -w 1 ) . (3.5) Let u i J ∈ V p J be arbitrary. We introduce ρ i J,alg ∈ V p J by ρ i J,alg := ρ i 0 + J j=1 ρ i j , (3.6)
where

ρ i 0 ∈ V 1 0 solves (3.
3) and ρ i j ∈ V p j , for j ∈ {1, . . . , J -1}, and ρ i J ∈ V p J are given by

ρ i j := 1 w 1 a∈Vj-s ρ i j,a , 1 ≤ j ≤ J, (3.7)
with the local contributions ρ i j,a ∈ V a j,s given by patch problems, ∀v j,a ∈ V a j,s

(∇ρ i j,a , ∇v j,a ) ω a j,s = (f, v j,a ) ω a j,s -(∇u i J , ∇v j,a ) ω a j,s - 1 w 2 j-1 k=0 (∇ρ i k , ∇v j,a ) ω a j,s . (3.8) Remark 3.3 (Construction of ρ i J,alg
). The construction (3.6)-(3.8) of ρ i J,alg can be seen as an approximation of ρi J,alg from (3.1) by one iteration of a V-cycle multigrid, with no presmoothing and a single postsmoothing step, corresponding to a "damped" additive Schwarz iteration, with the damping factor determined by the weights w 1 and w 2 . The subdomains for the Schwarz method correspond to the patch domains where the local problems in (3.8) are defined. Two patch options as in Figure 2 are considered. In particular, for p = 1 and "small" patches, s = 0 (Figure 2, left), this corresponds to one step of the Jacobi (diagonal) smoother, whereas when p = p > 1, the smoother is block Jacobi. A weighted variant of Definition 3.2 is tested in section 6.

Remark 3.4 (Value of the damping parameter). Condition (3.5) is based on the proofs in section 7 below, where the use of appropriate damping seems crucial. This is what is also indicated numerically to be needed in our approach. Possible combinations of the damping weights satisfying (3.5) include, for example,

w 1 = J(d + 1) and w 2 = 1, (3.9a
)

w 1 = d + 1 and w 2 = J, (3.9b
)

w 1 = w 2 = J(d + 1), (3.9c
)

w 1 = 1 and w 2 = ∞, (3.9d
)

w 1 = 4 √ J and w 2 = ∞. (3.9e)
Examples (3.9a)-(3.9c) above result in a procedure that is additive patchwise and multiplicative levelwise. Examples (3.9d)-(3.9e), in turn, result in a completely additive patchwise and levelwise procedure, which is fully parallelizable. We also note that when the intermediate polynomial degree is p = 1 and for any choice with w 2 = 1, the smoothing resulting from Definition 3.2 is local with respect to mesh T 0 for graded meshes; it is actually only performed there where the meshes T j , j ≥ 1, are different from T 0 .

An a posteriori estimator on the algebraic error and a multilevel solver

We now present how the residual lifting ρ i J,alg of Definition 3.2 can be used to define an a posteriori estimator as well as a multilevel solver.

A posteriori estimate on the algebraic error

We begin by introducing η i alg , an a posteriori estimator defined using the residual lifting ρ i J,alg .

Definition 4.1 (Lower bound algebraic error estimator). Let u i J ∈ V p J be arbitrary and let ρ i J,alg be the algebraic residual lifting given by Definition 3.2. If ρ i J,alg = 0, we define the lower bound algebraic error estimator η i alg := 0. Otherwise, set

η i alg := (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) ∇ρ i J,alg . (4.1)
The estimator η i alg is immediately a guaranteed lower bound on the algebraic error; cf., e.g., [38, Theorem 5.3].

Lemma 4.2 (Guaranteed lower bound on the algebraic error). There holds

∇(u J -u i J ) ≥ η i alg . (4.2) Proof. Note that if ρ i J,alg = 0, then ∇(u J -u i J ) ≥ 0 = η i alg . Otherwise ∇(u J -u i J ) = max v J ∈V p J , ∇v J =0 (∇(u J -u i J ), ∇v J ) ∇v J ≥ (∇(u J -u i J ), ∇ρ i J,alg ) ∇ρ i J,alg (2.4) (4.1) = η i alg .

Multilevel solver

We will now reuse the construction of ρ i J,alg given in Definition 3.2 to obtain an approximation of u J on a next step in view of constructing a multilevel solver. Note that for any u i J ∈ V p J , the lifting ρ i J,alg is built to approximate the algebraic error ρi J,alg given in (3.1), where we recall that u J = u i J + ρi J,alg . Thus, it seems reasonable to consider a linear iterative solver of the form

u i+1 J := u i J + λρ i J,alg , (4.3)
where λ ∈ R is a real parameter. The optimal choice of λ is given below.

:= [(f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg )]/ ∇ρ i J,alg
2 leads to minimal algebraic error with respect to the energy norm.

Proof. We write the algebraic error of the next iteration as a function of λ,

∇(u J -u i+1 J ) 2 = ∇(u J -u i J ) 2 -2λ(∇(u J -u i J ), ∇ρ i J,alg ) + λ 2 ∇ρ i J,alg 2 , (4.4)
and realize that this function has a minimum at

λ min = (∇(u J -u i J ), ∇ρ i J,alg ) ∇ρ i J,alg 2 (2.4) = (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) ∇ρ i J,alg 2 (4.1) = η i alg ∇ρ i J,alg . (4.5)
We are now ready to define our multilevel solver.

Definition 4.4 (Multilevel solver).

1. Initialize u 0 J ∈ V 1 0 as the solution of (∇u 0 J , ∇v 0 ) = (f, v 0 ) ∀v 0 ∈ V 1 0 .
2. Let i ≥ 0 be the iteration counter and let ρ i J,alg be constructed from u i J following Definition 3.2. When ρ i J,alg = 0, set u i+1 J := u i J and stop; then actually u i+1 J = u i J = u J . Otherwise, let

u i+1 J := u i J + (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) ∇ρ i J,alg 2 ρ i J,alg . (4.6)
Remark 4.5 (Multilevel solver). Note that the solver of Definition 4.4 is not initialized randomly but via a coarse solve. The descent direction is the residual lifting ρ i J,alg , constructed via a single V-cycle iteration with no presmoothing and one postsmoothing step, and the step-size is optimized via the line search (4.5). This minimalist and asymmetrical procedure will not be an issue in the forthcoming analysis. Remark 4.6 (Cost of one iteration). On each iteration of the developed solver, there are costs which correspond to those of standard multigrid methods: coarse solve (here with the lowest polynomial degree) and interlevel transfer operations. The crucial difference is in the smoothing cost. While we prove below that our solver is p-robust and only mildly depends on h (since J ∼ |logh|), meaning the number of iterations will not degrade when p increases, the sizes of the local matrices used to solve the local problems (3.8) increase (in 2D approximately as p 2 ). This induces a significant computational, but perfectly parallelizable, cost for higher p. Other cheaper options may be developed to bypass the local problems, for example, in the spirit of Papež and Vohralík [START_REF] Papež | Inexpensive guaranteed and efficient upper bounds on the algebraic error in finite element discretizations[END_REF]. Recall, however, that there is only one smoothing per iteration in our approach.

Main results

In this section, we present the main results concerning our a posteriori estimator η i alg of Definition 4.1 and our multilevel solver of Definition 4.4. We shall also see how these two main results are related.

For the estimator the following holds.

Theorem 5.1 (p-robust reliable and efficient bound on the algebraic error). Let Assumption 2.1 hold, together with either Assumption 2.2 or 2.3. Let u J ∈ V p J be the (unknown) solution of (2.4) and let u i J ∈ V p J be arbitrary, i ≥ 0. Let η i alg be given by Definition 4.1. Then, in addition to ∇(u J -u i J ) ≥ η i alg of (4.2), there holds

η i alg ≥ β ∇(u J -u i J ) , (5.1)
where 0 < β < 1 only depends on the space dimension d, the mesh shape regularity parameter κ T , and the number of mesh levels J, as well as on the mesh refinement parameter C ref and the quasi-uniformity parameter C qu if Assumption 2.2 holds, or on the coarse mesh and the local quasi-uniformity parameters C 0 qu and C loc qu if Assumption 2.3 holds. For all weights satisfying (3.5), there holds β ≥ J -5/2 β * with β * independent of the number of levels J. Better bounds hold for the weights of Remark 3.4; see Example 7.8 below for details.

The theorem allows us to write η i alg as a two-sided bound of the algebraic error (up to the generic constant β for the upper bound), meaning that the estimator is robustly efficient with respect to the polynomial degree p. We can also reinterpret this result as follows.

Remark 5.2 (Angle between the error and the descent direction). Note that if we rewrite (5.1) by plugging in the Definition 4.1 of η i alg , when u J -u i J = 0 and ρ i J,alg = 0, we have

η i alg ∇(u J -u i J ) (4.1) = (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) ∇(u J -u i J ) ∇ρ i J,alg (2.4)
= (∇(u J -u i J ), ∇ρ i J,alg ) ∇(u J -u i J ) ∇ρ i J,alg (5.1)
≥ β > 0.

This can be compared to classical results in line search methods (see, e.g., Nocedal and Wright [35, Chapter 3.2]) of boundedness away from zero of the cosine of the angle between the vector to be minimized (here u J -u i J ) and the descent direction (here the lifting ρ i J,alg ).

For the solver, in turn, we have the following.

Theorem 5.3 (p-robust error contraction of the multilevel solver). Let Assumption 2.1 hold, together with either Assumption 2.2 or 2.3. Let u J ∈ V p J be the (unknown) solution of (2.4) and let u i J ∈ V p J be arbitrary, i ≥ 0. Take u i+1 J to be constructed from u i J using one step of the multilevel solver of Definition 4.4 by (4.6). Then there holds

∇(u J -u i+1 J ) ≤ α ∇(u J -u i J ) , (5.2)
where 0 < α < 1 is given by α = 1 -β 2 with β the constant from (5.1).

In the above theorem, α is a bound on the algebraic error contraction factor at each step i. Looking at the dependencies of α, we see that the solver of Definition 4.4 contracts the algebraic error at each iteration step in a robust way with respect to the polynomial degree p.

Theorems 5.1 and 5.3 are connected as follows.

Corollary 5.4 (Equivalence of the p-robust estimator efficiency and p-robust solver contraction). Let the assumptions of Theorems 5.1 and 5.3 be satisfied. Then (5.1) holds if and only if (5.2) holds, and

β = √ 1 -α 2 .
Proof. Let u J ∈ V p J be the solution of (2.4), let u i J ∈ V p J be arbitrary, and let u i+1 J ∈ V p J be constructed from u i J by our multilevel solver of Definition 4.4. First, we write the relation between the algebraic errors associated to u i+1 J and u i J . Case ρ i J,alg = 0. Using (4.4) and (4.5), we see

∇(u J -u i+1 J ) 2 = ∇(u J -u i J ) 2 -(η i alg ) 2 .
(5.3)

Case ρ i J,alg = 0. By Definitions 4.4 and 4.1, we have u i+1 J = u i J and η i alg = 0. In particular, this means that ∇(u J -u i+1 J ) = ∇(u J -u i J ) , so that (5.3) still holds. The above observations allow us to write, in any case, starting from (5.2) with 0 < α < 1,

∇(u J -u i+1 J ) 2 ≤ α 2 ∇(u J -u i J ) 2 (5.3) ⇔ ∇(u J -u i J ) 2 -(η i alg ) 2 ≤ α 2 ∇(u J -u i J ) 2 ⇔ ∇(u J -u i J ) 2 (1 -α 2 ) ≤ (η i alg ) 2 ,
which is (5.1) with

β 2 = 1 -α 2 .
In view of Corollary 5.4, we will prove in section 7 below only Theorem 5.1. Importantly, the following also holds.

Corollary 5.5 (Equivalence of vanishing algebraic lifting with the solver reaching the solution). Let the assumptions of Theorems 5.1 and 5.3 be satisfied. Then ρ i J,alg = 0 if and only if u i+1 J = u i J = u J . Finally, by the proofs in section 7, the algebraic error is also equivalent to a localized a posteriori error estimate.

Corollary 5.6 (p-robust localized reliable and efficient a posteriori estimate on the algebraic error). Let the assumptions of Theorem 5.1 be satisfied. Let ρ i J,alg be the algebraic residual lifting constructed in Definition 3.2. Then

∇(u J -u i J ) 2 ≤ C 2 1 ∇ρ i 0 2 + J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s ≤ C 2 2 ∇(u J -u i J ) 2 , (5.4)
where C 2 = 1 β and C 1 is identified in section 7.6. Equivalence (5.4) gives us an idea where the algebraic error is situated levelwise and patchwise. This information can be exploited to tackle problematic areas adaptively, which is the subject of forthcoming works.

Numerical experiments

In this section we report some numerical illustrations of the theoretical results of section 5. In particular, we focus on the p-robustness. In the following tests, we consider the model problem (2.1) with three different choices of the domain Ω ⊂ R 2 and of the exact solution u:

Sine:

u(x, y) := sin(2πx) sin(2πy), Ω := (-1, 1) 2 . (6.1)

Peak: u(x, y) := x(x -1)y(y -1)e -100((x-0.5) 2 -(y-0.117) 2 ) , Ω := (0, 1) 2 . (6.2)

L-shape: u(r, θ) := r 2/3 sin(2θ/3), Ω := (-1, 1)

2 \([0, 1] × [-1, 0]). (6.3)
For the L-shape problem (6.3), we impose an inhomogeneous Dirichlet boundary condition corresponding to the exact solution, which is expressed here in polar coordinates. For each of the test cases, we start with an initial Delaunay triangulation of Ω. Then we consider J uniform refinements where all triangles are decomposed into four congruent subtriangles. Implementation-wise, we opt for Lagrange basis functions with nonuniformly distributed nodes because of their better behavior with respect to high-order methods; see Warburton [START_REF] Warburton | An explicit construction of interpolation nodes on the simplex[END_REF]. Recall that this choice has no influence on the theoretical results of section 5 as well as presented numerical results (in exact arithmetic). Though it is not the focus of this work, we also remark that our solver can be implemented in a matrix-free way and can also be parallelized.

The contraction factor of the solver of Definition 4.4 on each step i is given by ∇(u J -u i+1 J ) / ∇(u Ju i J ) , and, as stated in Corollary 5.4, it reveals the efficiency of the a posteriori estimator η i alg of Definition 4.1. Keeping this in mind, we only focus on the solver and the contraction factor. We will follow a common choice for the stopping criterion, with the notation of section 2.3:

F J -A J U is J F J ≤ 10 -5 F J -A J U 0 J F J . (6.4)
We also introduce the average error contraction factor

ᾱ := 1 i s is-1 i=0 ∇(u J -u i+1 J ) ∇(u J -u i J ) . (6.5)
We expect a p-robust solver to converge in a similar number of iterations and have similar error contraction factors at all iterations for different polynomial degrees p. The tests below cover different numbers of mesh levels J = 3, 4, 5, polynomial degrees p = 1, 3, 6, 9, and the "small" as well as the "large" patches as in Figure 2.

Performance of the damped additive Schwarz (dAS) construction of the solver

A crucial component in the definition of our a posteriori estimator and multilevel solver is the construction of the residual lifting ρ i J,alg of Definition 3.2, where we have used damped additive Schwarz (dAS) to cope with overlapping:

dAS: ρ i J,alg := J j=0 ρ i j and ρ i j := 1 w 1 a∈Vj-s ρ i j,a , 1 ≤ j ≤ J, (6.6)
(∇ρ i j,a , ∇v j,a ) ω a j,s = (f, v j,a ) ω a j,s -(∇u i J , ∇v j,a ) ω a j,s -

1 w 2 j-1 k=0 (∇ρ i k , ∇v j,a ) ω a j,s .
For the three test cases we consider three different choices of the damping weights which satisfy condition (3.5) (see Remark 3.4):

for problem (6.1): w 1 = J(d + 1) and w 2 = 1;

for problem (6.2):

w 1 = 4 √ J and w 2 = ∞;
for problem (6.3): w 1 = d + 1 and w 2 = J.

Recall that the choice w 2 = ∞ means that the construction of the lifting ρ i J,alg can be implemented completely in parallel, levelwise as well as patchwise.

The results are presented in Figures 345and in Table 1. They confirm the expected complete independence of the polynomial degree p for our multilevel solver which uses the construction dAS (6.6) of the lifting. Actually, we observe better contraction factors for higher polynomial degrees. ) for p = p in (2.10a), "small" (left) and "large" (right) patches, and stopping criterion (6.4). Top: error contraction factors ∇(u J -u i+1 J ) / ∇(u J -u i J ) . Bottom: relative algebraic error ∇(u J -u i J ) / ∇u J . ) for p = p in (2.10a), "small" (left) and "large" (right) patches, and stopping criterion (6.4). Top: error contraction factors ∇(u J -u i+1 J ) / ∇(u J -u i J ) . Bottom: relative algebraic error ∇(u J -u i J ) / ∇u J .

An inferior quality of the contraction factors for the case of p = 1 and the use of damping factors w 1 = J(d + 1) and w 2 = 1 appears. This is in line with some precedents in literature, where numerically probust solvers also perform worse for order 1 approximations; we mention, for example, Griebel, Oswald, and Sine problem (6.1) Peak problem (6.2) L-shape problem (6.3) w1 = J(d + 1), w2 = 1 w1 = 4 √ J, w2 = ∞ w1 = d + 1, w2 = J "small" "large" "small" "large" "small" "large" 1]. Recall that we consider no presmoothing and only one postsmoothing step; an important drop in the number of iterations appears if more smoothing steps are employed, which will be explored below. Another observation is that the number of iterations depends on the number of mesh levels J, in accordance with the theoretical result of section 7, even though rather mildly. ) for p = p in (2.10a), "small" (left) and "large" (right) patches, and stopping criterion (6.4). Top: error contraction factors ∇(u J -u i+1 J ) / ∇(u J -u i J ) . Bottom: relative algebraic error ∇(u J -u i J ) / ∇u J .

J p DoF is ᾱ is ᾱ is ᾱ is ᾱ is ᾱ is ᾱ 3 1 5e 3 48
The behavior of the contraction factor in each iteration in Figures 345appears quite different. This seems to be related partly to the smoothness of the problem and partly to choice of the damping weights. We explore this in more detail in Figure 6 by using different choices of the weights and number of postsmoothing steps ν. In particular the degradation of the contraction factors observed in Figure 3 disappears when employing more smoothing steps. : Sine problem (6.1), "small" patches, p = p: study of the contraction factor behavior with respect to the number of postsmoothing steps and damping weights for the solver (4.6).

Performance of the weighted restrictive additive Schwarz (wRAS) construction of the solver

As observed in the literature, replacing the damping with parameter w 1 in (3.7) by hat function weighting via a restrictive additive Schwarz often performs better; cf. Cai and Sarkis [START_REF] Cai | A restricted additive Schwarz preconditioner for general sparse linear systems[END_REF], Efstathiou and Gander [START_REF] Efstathiou | Why restricted additive Schwarz converges faster than additive Schwarz[END_REF], or Loisel, Nabben, and Szyld [START_REF] Loisel | On hybrid multigrid-Schwarz algorithms[END_REF]. Thus, in addition to the dAS construction (6.6), we now numerically also explore the weighted restricted additive Schwarz (wRAS) construction of the lifting ρ i J,alg :

wRAS: ρ i J,alg := J j=0 ρ i j and ρ i j := a∈Vj-s I q j (ψ a j-s ρ i j,a ), 1 ≤ j ≤ J, (6.7)

(∇ρ i j,a , ∇v j,a ) ω a j,s = (f, v j,a ) ω a j,s -(∇u i J , ∇v j,a ) ω a j,s - j-1 k=0 (∇ρ i k , ∇v j,a ) ω a j,s
with q = p except for j = J where q = p, we denote by I q j the P q Lagrange interpolation operator on the mesh level j, i.e., I q j : C 0 (Ω) → V q j , I q j (v) preserves the values of v in the nodes corresponding to the Lagrange degrees of freedom. No damping weights are to be chosen here.

We summarize the results obtained for each of the problems (6.1)-( 6.3) in Table 2. In addition to one postsmoothing step, ν = 1, we also present the results for ν = 3 postsmoothing steps. In both cases, no presmoothing has been employed. In the last two columns for each problem, we present a comparison of our solver of Definition 4.4 employing (6.7) with two standard smoothers for multigrid, namely the Jacobi (J) and the Gauss-Seidel (GS) ones. Here, we employ no presmoothing step, one postsmoothing step, and a coarse solve with polynomials of order 1 as in (3.3) to compare with our approach.

The results using the wRAS (6.7) construction of the lifting indicate an improvement in the error contraction factors with respect to dAS (6.6) of section 6.1 and, moreover, present a complete numerical independence of the number of levels J. Furthermore, the iteration numbers drop by at least half when three postsmoothing steps are employed. In contrast to these results, we see that the multigrid with standard smoothers degrades violently with respect to the polynomial degree p. Note in this respect that for p = 1, the only difference between wRAS of (6.7) with small patches and ν = 1 and standard Jacobi lies in the optimally chosen step-size of Lemma 4.3. This gives a spectacular gain in the number of iterations, and makes the method convergent even when the standard Jacobi fails.

Sine problem (6.1) Peak problem (6.2) L-shape problem (6.3) wRAS MG wRAS MG wRAS MG "small" "large" ν = 1 "small" "large" ν = 1 "small" "large" Table 2: Number of iterations needed to reach the stopping criterion (6.4): wRAS construction (6.7), problems (6.1)-( 6.3), p = p in (2.10a), "small" and "large" patches, ν postsmoothing steps, and standard multigrid method with piecewise affine coarse solve (3.3), initialized by the coarse grid solution, no presmoothing, one postsmoothing step, and Jacobi (J) and Gauss-Seidel (GS) smoothers.

ν = 1 J p ν = 1 ν = 3 ν = 1 ν = 3 J GS ν = 1 ν = 3 ν = 1 ν = 3 J GS ν = 1 ν = 3 ν = 1 ν = 3 J GS

Comparison with other multilevel solvers

Some recent comparisons of state-of-the-art solvers for Poisson problems with multigrid methods in the highorder setting include Gholami et al. [START_REF] Gholami | or multigrid? A comparative study of state-of-the-art Poisson solvers for uniform and nonuniform grids in the unit cube[END_REF], Sundar, Stadler, and Biros [START_REF] Sundar | Comparison of multigrid algorithms for high-order continuous finite element discretizations[END_REF], and Kronbichler and Wall [START_REF] Kronbichler | A performance comparison of continuous and discontinuous Galerkin methods with fast multigrid solvers[END_REF].

In [START_REF] Sundar | Comparison of multigrid algorithms for high-order continuous finite element discretizations[END_REF], it was in particular reported that none of the methods considered behaves fully independently of the polynomial degree. In this subsection, we compare our developments with four well-established methods. We focus on the number of iterations, but we also indicate CPU times of our vectorized Matlab implementation1 , trusting the reader to understand the trickiness inherent in such implementation-and machine-dependent measurements. The timings below involve the solution time only; i.e., they do not include the assembly time of the matrices. The methods we consider for the comparison are as follows:

wRAS, ν = 1 (∼MG(0,1)-bJ): Definition 4.4, small patches, p = p in (2.10a) (to illustrate the associated space hierarchy, we write "1, p → p"), wRAS construction (6.7).

wRAS, ν = 3 (∼MG(0,3)-bJ): Definition 4.4, small patches, p = p in (2.10a) ("1, p → p"), wRAS construction (6.7), three postsmoothing steps employed.

wRAS, ν = 1 (∼MG(0,1)-bJ): Definition 4.4, small patches, p = 1 in (2.10a) ("1 → 1, p"), wRAS construction (6.7).

wRAS, ν = 3 (∼MG(0,3)-bJ): Definition 4.4, small patches, p = 1 in (2.10a) ("1 → 1, p"), wRAS construction (6.7), three postsmoothing steps employed.

PCG(MG(3,3)-bJ): Preconditioned conjugate gradient solver; the preconditioner is multigrid V-cycle [START_REF] Antonietti | A uniform additive Schwarz preconditioner for high-order discontinuous Galerkin approximations of elliptic problems[END_REF][START_REF] Antonietti | A uniform additive Schwarz preconditioner for high-order discontinuous Galerkin approximations of elliptic problems[END_REF] with weighted restrictive additive Schwarz (block Jacobi) smoother associated to small patches; the space hierarchy relies on order p discretization, including the coarsest space ("p → p"); the iterations start with the zero vector. This choice of solver is motivated by Antonietti and Pennesi [START_REF] Antonietti | V -cycle multigrid algorithms for discontinuous Galerkin methods on non-nested polytopic meshes[END_REF], adapted to the conforming finite elements setting.

MG(1,1)-PCG(iChol):

Multigrid solver V-cycle(1,1); the smoother is PCG with incomplete zero level fill-in Cholesky preconditioner; the space hierarchy is of increasing order: from order 1 for the coarsest level to order p for the finest level ("1 p"); the iterations start with the zero vector. This choice of solver is motivated by Botti, Colombo, and Bassi [START_REF] Botti | h-multigrid agglomeration based solution strategies for discontinuous Galerkin discretizations of incompressible flow problems[END_REF], adapted for a symmetric setting. MG(0,1)-bGS: Multigrid solver V-cycle(0,1); the smoother is block Gauss-Seidel associated to small patches; the space hierarchy consists of order 1 for all levels except the finest level, which is of order p ("1 → 1, p"), i.e., as in (2.10a) with p =1; the iterations start with the zero vector. This choice of the solver is motivated by NGSolve [START_REF] Schöberl | C++11 Implementation of Finite Elements in NGsolve[END_REF]; however, the multigrid is used here as a solver instead of a preconditioner.

MG(0,3)-bGS: Multigrid solver analogous to MG(0,1)-bGS where now three postsmoothing steps are employed. MG(3,3)-GS: Multigrid solver V-cycle(3,3); the smoother is standard Gauss-Seidel; the space hierarchy is of increasing order: from order 1 for coarse level to order p for the finest level ("1 p"); the iterations start with the zero vector. Table 3: Comparison of various multilevel solvers (described in section 6.3) for the L-shape case (6.3), i s is the number of iterations to reach the stopping criterion (6.4).

wRAS wRAS wRAS wRAS PCG(MG MG(1,1)-MG(0,1)-MG(0,3)-MG(3,3)- ν = 1 ν = 3 ν = 1 ν = 3 (3,3)-bJ) PCG(iChol) bGS bGS GS 1, p → p 1, p → p 1 → 1, p 1 → 1, p p → p 1 p 1 → 1, p 1 → 1, p 1 p J p i s time i s time i s time i s time i s time i s time i s time i s time i s time 3
As one can see from Table 3, the presented methods split into two groups: numerically p-robust (wRAS, PCG(MG-bJ), MG-bGS) and not (MG-PCG(iChol), MG-GS). Note that the choice of three pre-and three postsmoothing steps makes every iteration of the methods PCG(MG(3,3)-bJ) and MG(3,3)-GS considerably more expensive than those of the methods wRAS and MG-bGS with ν = 1, where the minimalist (0,1) choice is sufficient. The variants wRAS and MG-bGS with ν = 3 are also cheaper. In addition, in PCG(MG(3,3)-bJ), the coarse grid correction is more expensive as it uses order p approximations. The inversion of the Jacobi blocks in PCG(MG(3,3)-bJ) on the finest level J, corresponds to solving the patch problems of order p as in (3.8), so that its cost is the same as for the local problems of wRAS. As for MG(1,1)-PCG(iChol), we find the method to be quite satisfactory for lower-order approximations and small J, but as soon as p and J increase, the number of iterations degrades considerably. In contrast to wRAS, MG-bGS is a multiplicative Schwarz method and is thus less suitable for parallelization. Finally, the classical MG(3,3)-GS is a combination of h-and p-multigrid and gives the best timings in our experiments. The numbers of pre-and postsmoothing steps, however, remain parameters, and their tuning might not be straightforward in order to get an efficient and numerically robust multigrid solver in general (cf. the poor results of the very similar-up to the different number of pre-and postsmoothing steps and a stronger hierarchy-MG(0,1)-GS version in Table 2). The Gauss-Seidel smoother used therein again makes the method harder to parallelize. Lemma 7.1 (Estimate on ∇ρ i J,alg and ∇ρ i j by patchwise contributions). Let ρ i J,alg and ρ i j for j ∈ {1, . . . , J} be given by Definition 3.2. Then

∇ρ i j 2 ≤ d + 1 w 2 1 a∈Vj-s ∇ρ i j,a 2 ω a j,s , (7.1)
∇ρ i J,alg 2 ≤ C 2 max (w 1 ) ∇ρ i 0 2 + J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s , (7.2)
where 2 ≤ C 2 max (w 1 ) := 2 max 1,

J(d + 1) w 2 1 ≤ 2J(d + 1). (7.3)
Proof. Definition 3.2 and inequality

| d+1 k=1 a k | 2 ≤ (d + 1) d+1 k=1 |a k | 2 lead to ∇ρ i j 2 = K∈Tj-s ∇ρ i j 2 K = K∈Tj-s 1 w 1 a∈V K ∇ρ i j,a 2 K ≤ d + 1 w 2 1 K∈Tj-s a∈V K ∇ρ i j,a 2
K = d + 1 w 2 1 a∈Vj-s ∇ρ i j,a 2 ω a j,s .
Note that this allows us to write

J j=1 ∇ρ i j 2 ≤ J J j=1 ∇ρ i j 2 ≤ J(d + 1) w 2 1 J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s . (7.4)
This property together with some simple manipulations, gives the second estimate:

∇ρ i J,alg 2 ≤ 2 ∇ρ i 0 2 + 2 J j=1 ∇ρ i j 2 (7.4) ≤ 2 ∇ρ i 0 2 + 2J(d + 1) w 2 1 J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s ≤ 2 max 1, J(d + 1) w 2 1 ∇ρ i 0 2 + J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s
.

The bounds (7.3) on C 2 max (w 1 ) are easily obtained by using w 1 ≥ 1 requested in Definition 3.2.

J j=1 (f, ρ i j ) -(∇u i J , ∇ρ i j ) (3.3) (3.7)
= ∇ρ i 0 2 + 1 w 1 J j=1 a∈Vj-s (f, ρ i j,a ) ω a j,s -(∇u i J , ∇ρ i j,a ) ω a j,s (3.8)
= ∇ρ i 0 2 + 1 w 1 J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s + 1 w 2 j-1 k=0 (∇ρ i k , ∇ρ i j,a ) ω a j,s (3.7)
= ∇ρ i (∇ρ i k , ∇ρ i j ).

The first two terms above are of the right form to prove the result, but one needs to be a bit more careful with the third one. We estimate it using Young's inequality and the sum interchange

J j=2 j-1 k=1 = J-1 k=1
J j=k+1 ; Young's parameter µ is picked later to control the dependence on J of the final estimate. We have

1 w 2 J j=1 j-1 k=0 (∇ρ i k , ∇ρ i j ) = 1 w 2   J j=2 j-1 k=1 (∇ρ i k , ∇ρ i j ) + J j=1 (∇ρ i 0 , ∇ρ i j )   ≥ 1 w 2 J j=2 j-1 k=1 - 1 2 ∇ρ i k 2 - 1 2 ∇ρ i j 2 + 1 w 2 J j=1 - 1 2µ ∇ρ i 0 2 - µ 2 ∇ρ i j 2 = - 1 2w 2 J j=1 (J -j) ∇ρ i j 2 - 1 2w 2 J j=1 (j -1) ∇ρ i j 2 - J 2µw 2 ∇ρ i 0 2 - µ 2w 2 J j=1 ∇ρ i j 2 ,
where we added the terms in the sum corresponding to k = J and j = 1 since they are zero, and then renamed the summation index when there is no confusion. Picking Young's inequality parameter µ = 2J 3w2 , a few more manipulations on the right-hand side give us

1 w 2 J j=1 j-1 k=0 (∇ρ i k , ∇ρ i j ) ≥ - 3 4 ∇ρ i 0 2 - J -1 + 2J 3w2 2w 2 J j=1 ∇ρ i j 2 (7.1) ≥ - 3 4 ∇ρ i 0 2 - J(1 + 2 3w2 )(d + 1) 2w 2 w 2 1 J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s . 1 4 ∇ρ i 0 2 + 1 w 1 - J(1 + 2 3w2 )(d + 1) 2w 2 w 2 1 J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s (7.6) ≥ C 2 min (w 1 , w 2 )   ∇ρ i 0 2 + J j=1 a∈Vj-s ∇ρ i j,a 2 ω a j,s   . (7.7)
The upper bound on C 2 min (w 1 , w 2 ) in (7.6) is immediate due to the minimum in its expression, while the lower bound is obtained by rewriting condition (3.5) on w 2

w 2 ≥ 5J 2 (d + 1) 2 w 1 6J(d + 1) 1 -w1 6J(d+1) ⇔ 1 w 1 - 5 6
J(d + 1) w 2 w 2 1 ≥ 1 6J(d + 1) . (7.8)
As also w 2 ≥ 1,

1 w 1 - J(1 + 2 3w2 )(d + 1) 2w 2 w 2 1 ≥ 1 w 1 - J(1 + 2 3 )(d + 1) 2w 2 w 2 1 = 1 w 1 - 5
= v # J + b∈V J v p J,b , where v # J ∈ V 1 J and v p J,b ∈ V b J,0 , b ∈ V J ,
and this decomposition is stable in the sense ∇v

# J 2 + b∈V J ∇v p J,b 2 ω b J,0 ≤ C 2 SD ∇v J 2 , (7.9)
where C SD ≥ 1 only depends on the mesh shape regularity parameter κ T and space dimension d.

Similarly to [50, Lemma 3.1 and Theorem 3.1], in the case of quasi-uniform meshes with bounded refinement strength, we have the following. Lemma 7.4 (P 1 -multilevel stable decomposition for quasi-uniform meshes). For all v # J ∈ V 1 J , there exists a multilevel piecewise affine decomposition v

# J = v 1 0 + J j=1 b∈Vj v 1 j,b with v 1 0 ∈ V 1 0 and v 1 j,b ∈ V b j,0 = P 1 (T b j,0 ) ∩ H 1 0 (ω b j,0
). Under Assumption 2.2, this decomposition is stable as

∇v 1 0 2 + J j=1 b∈Vj ∇v 1 j,b 2 ω b j,0 ≤ C 2 MD ∇v # J 2 , (7.10)
where C MD ≥ 1 only depends on the space dimension d, the mesh shape regularity parameter κ T , the maximum strength of refinement parameter C ref , and the quasi-uniformity parameter C qu .

Proof.

Let v # J ∈ V 1 J .
We first apply a levelwise decomposition that follows from Xu, Chen, and Nochetto in [START_REF] Xu | Optimal multilevel methods for H(grad), H(curl), and H(div) systems on graded and unstructured grids[END_REF]Lemma 3.1] by keeping the gradient on level zero. This gives us v

# J = J j=0 v 1 j with v 1 j ∈ V 1 j such that ∇v 1 0 2 + J j=1 h -2 j v 1 j 2 ≤ C 2 ml ∇v # J 2 , (7.11)
where C ml ≥ 1 has the same dependencies as C MD .

We further decompose each of the above v 1 j ∈ V 1 j , j ≥ 1, into patchwise components. For this purpose, we use the standard nodal decomposition v 1 j = b∈Vj v 1 j,b , where v 1 j,b = v 1 j (b)ψ j,b belongs to the local space V b j,0 for p = p = 1. By stability of the P 1 nodal decomposition,

b∈Vj v 1 j,b 2 ω b j,0 ≤ C 2 nd v 1 j 2 , (7.12)
where C nd only depends on the space dimension d and the mesh shape regularity parameter κ T . This can, for instance, be shown by considering a patch ω b j,0 and an element K contained in the patch. Since

v 1 j,b = v 1 j (b)ψ j,b ∈ V b j,0
and by equivalence of norms in finite dimension, we have

v 1 j (b)ψ j,b ω b j,0 ≈ v 1 j (b)ψ j,b K ≤ v 1 j (b)ψ j,b ∞,K |K| 1 2 ≤ b∈V K v 1 j (b)ψ j,b ∞,K |K| 1 2 = v 1 j ∞,K |K| 1 2 v 1 j K . (7.13)
The result (7.12) is obtained by summing both sides over all vertices. Now, the claim (7.10) follows by using an inverse inequality on patches, the quasi-uniformity of the meshes, the above decompositions as

J j=1 b∈Vj ∇v 1 j,b 2 ω b j,0 ≤ C 2 inv J j=1 b∈Vj h -2 ω b j,0 v 1 j,b 2 ω b j,0 (2.8) ≤ C 2 inv C -2 qu J j=1 h -2 j b∈Vj v 1 j,b 2 ω b j,0 (7.12)
≤ C 2 inv C -2 qu C 2 nd J j=1 h -2 j v 1 j 2 (7.11) ≤ C 2 inv C -2 qu C 2 nd C 2 ml ∇v # J 2
and by summing the left-hand side with ∇v 1 0 2 , which satisfies a similar bound from (7.11). We set

C 2 MD := C 2 inv C -2 qu C 2 nd C 2 ml + C 2 ml to obtain the result.
By [START_REF] Xu | Optimal multilevel methods for H(grad), H(curl), and H(div) systems on graded and unstructured grids[END_REF]Theorem 4.3], in the case of graded meshes, we have the following.

Lemma 7.5 (P 1 -multilevel stable decomposition for graded meshes). For all v # J ∈ V 1 J , there exists a multilevel piecewise affine decomposition v

# J = v 1 0 + J j=1 b∈Vj v 1 j,b with v 1 0 ∈ V 1 0 and v 1 j,b ∈ V b j,0 = P 1 (T b j,0 )∩ H 1 0 (ω b j,0
). Under Assumption 2.3, this decomposition is stable as

∇v 1 0 2 + J j=1 b∈Vj ∇v 1 j,b 2 ω b j,0 ≤ C 2 MD ∇v # J 2 , (7.14)
where C MD ≥ 1 only depends on the space dimension d, the mesh shape regularity parameter κ T , the coarse mesh quasi-uniformity parameter C 0 qu , and the local quasi-uniformity parameter C loc qu .

Proof. Let v # J ∈ V 1 J . We apply the results on stable decomposition on graded meshes of Xu, Chen, and Nochetto in [START_REF] Xu | Optimal multilevel methods for H(grad), H(curl), and H(div) systems on graded and unstructured grids[END_REF]Theorem 4.3] adjusted to our setting. On the one hand, this gives us the decomposition

v # J = v 1 0 + J j=1 b∈Bj v 1 j,b + b∈V J v 1 J,b , where v 1 0 ∈ V 1 0 , ∀b ∈ B j , v 1 j,b ∈ V b j,0 for p = p = 1, and ∀b ∈ B J , v 1 J,b ∈ V b J,0 for p = 1.
On the other hand, the result also gives us the following stability inequality:

∇v 1 0 2 + J j=1 h -2 Bj b∈Bj v 1 j,b 2 + b∈V J h -2 ω b J,0 v 1 J,b 2 ω b J,0 ≤ C 2 gra ∇v # J 2 , (7.15)
where C gra has the same dependencies as C MD . First, since the mesh hierarchy is created via bisections, we have the local quasi-uniformity property (2.9). This, together with an inverse inequality and L 2 -stability as in (7.13), gives us

J j=1 b∈Bj ∇v 1 j,b 2 ω b j,0 ≤ C 2 inv J j=1 b∈Bj h -2 ω b j,0 v 1 j,b 2 ω b j,0 (2.9) ≤ C 2 inv (C loc qu ) -2 C 2 nd J j=1 h -2 Bj b∈Bj v 1 j,b 2
.

Second, we only need an inverse inequality to obtain

b∈V J ∇v 1 J,b 2 ω b J,0 ≤ C 2 inv b∈V J h -2 ω b J,0 v 1 J,b 2 ω b J,0
.

Third, we can sum together the above estimations and use (7.15) to obtain

∇v 1 0 2 + J j=1 b∈Bj ∇v 1 j,b 2 ω b j,0 + b∈V J ∇v 1 J,b 2 ω b J,0 ≤ C 2 MD ∇v # J 2 ,
where

C MD := C gra • max(C inv , C inv C nd (C loc qu ) -1 , 1) . Finally, since B j ⊂ V j , we can set v 1 j,b := 0 for b ∈ V j \B j and have a new decomposition v # J = v 1 0 + J j=1 b∈Vj v 1 j,
v J = v 1 0 + J j=1 b∈Vj v j,b , v 1 0 ∈ V 1 0 , v j,b ∈ V b j,0 (7.16)
stable as

∇v 1 0 2 + J j=1 b∈Vj ∇v j,b 2 ω b j,0 ≤ C 2 SMD ∇v J 2 , (7.17)
where C SMD := √ 2C SD C MD ≥ 1 only depends on the space dimension d, the mesh shape regularity parameter κ T , and, depending on whether Assumption 2.2 or 2.3 is satisfied, on either the quasi-uniformity parameter C qu and the maximum strength of refinement parameter C ref or the coarse mesh and local quasi-uniformity parameters C 0 qu ,C loc qu , respectively.

Proof. Let v J ∈ V p J and let us begin by applying the decomposition of Lemma 7.3. This gives v

J = v # J + b∈V J v p J,b with v # J ∈ V 1 J and v p J,b ∈ V b J,0 ∀b ∈ V J .
Then, we further decompose v # J using either Lemma 7.4 or 7.5, depending on whether Assumption 2.2 or 2.3 is satisfied. We obtain v

# J = v 1 0 + J j=1 b∈Vj v 1 j,b , where v 1 0 ∈ V 1 0 and v 1 j,b ∈ V b j,0 (actually V b j,0 ∩ P 1 (T b j,0 )). We set v j,b := v 1 j,b ∀b ∈ V j , j ∈ {1, . . . , J -1} and v J,b := v 1 J,b + v p J,b . Thus, we have v J = v 1 0 + J j=1 b∈Vj v j,b with v 1 0 ∈ V 1 0 and v j,b ∈ V b j,0 .
The stable decomposition results presented in the previous lemmas allow us to write

∇v 1 0 2 + J j=1 b∈Vj ∇v j,b 2 ω b j,0 ≤ 2 ∇v 1 0 2 + J j=1 b∈Vj ∇v 1 j,b 2 ω b j,0 + b∈V J ∇v p J,b 2 ω b J,0 (7.10) or (7.14) ≤ 2C 2 MD ∇v # J 2 + b∈V J ∇v p J,b 2 ω b J,0 (7.9) ≤ 2C 2 SD C 2 MD ∇v J 2 .

Upper bound on ∇ρ i J,alg

Recall that ρi J,alg , introduced in (3.1), is the unknown exact algebraic error. We now estimate ∇ρ i J,alg from above. We introduce some helpful notation first. For all a ∈ V j-s , let I a ⊂ V j be a set containing (fine-mesh) vertices of the interior of the patch ω a j,s such that {I a } a∈Vj-s cover V j and are mutually disjoint; if s = 0, we have I a = {a}. This allows us to write b∈Vj = a∈Vj-s b∈Ia . Moreover, since the indices of I a are localized in the interior of the patch ω a j,s , we have b∈Ia v j,b ∈ V a j,s when v j,b ∈ V b j,0 . Writing it this way will help us to apply the results on the p-robust stable decomposition of Lemma 7.6 given for "small" patches only to the "large"-patch setting as well. Proof. The main ingredient of the proof is to replace locally the uncomputable ρi J,alg = u J -u i J by the constructed local contributions ρ i j,a using the problems they solve on patches. We begin by using Proposition 7.6 applied to ρi J,alg ∈ V p J , writing ρi J,alg = e 0 + We will now estimate each of the above three terms using Young's inequality and patch overlap arguments as done in the proof of Lemma 7.1. First, we have

∇ρ i 0 , b∈V0 ∇e 0,b ≤ C 2 SMD 2 ∇ρ i 0 2 + 1 2C 2 SMD ∇e 0 2 .
For the second term, we similarly obtain

Proof of Theorem 5.1

The results of the previous subsections allow us now to give a concise proof of Theorem 5.1.

Proof of Theorem 5.1. Case ρ i J,alg = 0. By Definition 4.1 this means η i alg = 0, so that it suffices to show that u J = u i J in this case. We do this by using Lemmas 7.2 and 7.7, which lead to ∇(u J -u i J )

Figure 2 :

 2 Figure 2: Illustration of degrees of freedom (p = p = 2) for the space V bj,0 associated to the "small" patch T b j,0 (left) and for the space V a j,1 associated to the "large" patch T a j,1 (right). The mesh T j-1 and its refinement T j are defined in bold and dotted lines, respectively.

Lemma 4 . 3 (

 43 Optimal step-size). Consider a solver of the form (4.3) and suppose ρ i J,alg = 0. Then the choice λ

Figure 3 :

 3 Figure3: Sine problem (6.1), w 1 = J(d + 1), w 2 = 1: results of the solver (4.6) for p = p in (2.10a), "small" (left) and "large" (right) patches, and stopping criterion (6.4). Top: error contraction factors ∇(u J -u i+1 J ) / ∇(u J -u i J ) . Bottom: relative algebraic error ∇(u J -u i J ) / ∇u J .

Figure 4 :

 4 Figure 4: Peak problem (6.2), w 1 = 4√ J, w 2 = ∞: results of the solver (4.6) for p = p in (2.10a), "small" (left) and "large" (right) patches, and stopping criterion (6.4). Top: error contraction factors ∇(u J -u i+1 J ) / ∇(u J -u i J ) . Bottom: relative algebraic error ∇(u J -u i J ) / ∇u J .

Figure 5 :

 5 Figure 5: L-shape problem(6.3), w 1 = d + 1, w 2 = J: results of the solver (4.6) for p = p in (2.10a), "small" (left) and "large" (right) patches, and stopping criterion (6.4). Top: error contraction factors ∇(u J -u i+1 J ) / ∇(u J -u i J ) . Bottom: relative algebraic error ∇(u J -u i J ) / ∇u J .

Figure 6

 6 Figure6: Sine problem (6.1), "small" patches, p = p: study of the contraction factor behavior with respect to the number of postsmoothing steps and damping weights for the solver (4.6).

  3m 3 2.1m 12 2.3m 5 2.5m 2 3.1m +200 +1.0 h 6 2.2m 3 2.7m 8 1.3m

Lemma

.

  Finally, for the third term we additionally use the property w 2 ≥ 1 and rename summation indices when there is no confusion:

J 2 w2 ≤ J 2 and d+1 w 2 1 ≤

 21 d + 1, leading to the desired result:C 2 max (w 1 , w 2 ) ≤ 2(d + 1)C 2 SMD 1 + J 2 (d + 1) ≤ 2(d + 1)C 2 SMD 2J 2 (d + 1) .

Example 7 . 8 ( 1 12C 1 12C SMD J 2 ( 4 √ J and w 2 = ∞ : 1 8C

 7811241 Specific choices of weights). We illustrate here a bound on the efficiency factor β in (5.1) for different choices of the damping weights satisfying the compatibility condition (3.5) from Remark 3.4:w 1 = J(d + 1) and w 2 = 1 : SMD J 2 2(d + 1) 3 ≤ β.w 1 = d + 1 and w 2 = J : SMD J(d + 1) ≤ β.

7. 6

 6 Proof ofCorollary 5.6 If ρ i J,alg = 0, as a result of (7.21), we have ∇(u J -u i J ) = ∇ρ i = 0. Otherwise by(7.18) we set C 1 := C max (w 1 , w 2 ) and by (4.2) and(7.22), C 2 = 1 β = C max (w 1 )C max (w 1 , w 2 )C 2 min (w 1 , w 2 ) gives the result.

Table 1 :

 1

		0.79 34 0.70	74 0.85 43 0.75 38 0.75 20 0.56
	3 4e 4	23 0.63 24 0.59	60 0.83 36 0.70 28 0.68 18 0.53
	6 2e 5	23 0.63 22 0.55	58 0.82 34 0.68 27 0.69 16 0.49
	9 4e 5	23 0.63 19 0.50	58 0.82 31 0.65 25 0.69 14 0.46
	4 1 2e 4	52 0.80 40 0.74	87 0.87 48 0.77 39 0.76 23 0.60
	3 2e 5	27 0.68 26 0.60	66 0.84 41 0.72 28 0.70 23 0.60
	6 1e 6	26 0.66 24 0.57	68 0.84 38 0.70 29 0.72 20 0.58
	9 2e 5	26 0.67 21 0.53	70 0.84 33 0.67 28 0.72 18 0.55
	5 1 1e 5	56 0.81 43 0.75	97 0.88 52 0.78 40 0.76 25 0.63
	3 1e 6	32 0.73 28 0.61	72 0.85 44 0.74 30 0.72 27 0.65
	6 3e 6	29 0.71 26 0.58	78 0.86 42 0.72 31 0.74 25 0.63
	9 6e 6	29 0.71 24 0.56	84 0.86 38 0.71 30 0.74 21 0.59

dAS construction (6.6): problems (6.1)-(6.3), p = p in (2.10a), "small" and "large" patches. i s : the number of iterations needed to reach the stopping criterion (6.4). ᾱ: average error contraction factor given by (6.5).

Schweitzer

[23,

Table

1

] and Kronbichler and Wall

[START_REF] Kronbichler | A performance comparison of continuous and discontinuous Galerkin methods with fast multigrid solvers[END_REF] Table

  Now, we devise a p-robust multilevel stable decomposition. This decomposition relies on the one level p-robust stable decomposition given in Schöberl et al.[43, Proof of Theorem 2.1] and the piecewise affine multilevel decomposition in the spirit of Xu, Chen, and Nochetto [50, Theorems 3.1 and 4.3]. These results are presented below in the form of lemmas. Note that in the decomposition, only "small" patches are used, which will be sufficient for our purposes. Recall also the definition of the local spaces (2.12), which will be useful below. Hereafter, we always assume that Assumption 2.1 is satisfied.By [43, Proof of Theorem 2.1], we have the following.

	6	J(d + 1) w 2 w 2 1	(7.8) ≥	1 6J(d + 1)	.
	7.3 Polynomial-degree-robust multilevel stable decomposition
	Lemma 7.3 (One-level p-robust stable decomposition). For all v J ∈ V p J , there exists a finest-level decom-
	position v J

  b (reusing the notation) such that (7.14) holds. Proposition 7.6 (p-robust multilevel stable decomposition). Let v J ∈ V p J . Under either Assumption 2.2 or 2.3, there exists a decomposition

  7.7 (Estimating ∇ρ i J,alg by local contributions). Let ρi J,alg ∈ V p J be defined by (3.1). We have

					J				
	∇ρ i J,alg	2 ≤ C 2 max (w 1 , w 2 )	 ∇ρ i 0	2 +		∇ρ i j,a	2 ω a j,s	 ,	(7.18)
					j=1 a∈Vj-s
	where
	2(d + 1) ≤ C 2 max (w 1 , w 2 ) := 2(d + 1)C 2 SMD 1 +	J 2 w 2	max 1,	d + 1 1 w 2	≤ 4(d + 1) 2 C 2 SMD J 2 .	(7.19)

  Summing these components together, we can now pursue our main estimate: on both sides, we finally obtain the desired result.The lower bound on C 2 max (w 1 , w 2 ) in (7.[START_REF] Efstathiou | Why restricted additive Schwarz converges faster than additive Schwarz[END_REF]) is obtained by using C SMD ≥ 1 from Proposition 7.6. To derive the upper bound, we use the fact that weights of Definition 3.2 satisfy w 1 ≥ 1, w 2 ≥ 1. This gives

			2 SMD J 2 w 2	∇ρ i 0	2 +	d + 1 w 2 1	J j=1 a∈Vj-s	∇ρ i j,a	2 ω a j,s	+	1 4C 2 SMD	J j=1 b∈Vj	∇e j,b	2 ω b j,0	.
	∇ρ i J,alg	2	(7.19) ≤	C 2 max (w 1 , w 2 ) 2	∇ρ i 0	2 +	J j=1 a∈Vj-s		∇ρ i j,a	2 ω a j,s +	∇e 0	2 +	J j=1 b∈Vj 2C 2 SMD	∇e j,b	2 ω b j,0
			(7.17) ≤	C 2 max (w 1 , w 2 ) 2	∇ρ i 0	2 +	J j=1 a∈Vj-s	∇ρ i j,a	2 ω a j,s	+	1 2	∇ρ i J,alg	2 .
	After subtracting 1 2 ∇ρ i J,alg	2

  Case ρ i J,alg = 0. In this case, we combine the results of Lemmas 7.1, 7.2, and 7.7 , w 2 )C max (w 1 )C max (w 1 , w 2 ) , w 2 ) C max (w 1 )C max (w 1 , w 2 )

				2 (3.2) = ∇ρ i J,alg	2 (7.18) ≤ C 2 max (w 1 , w 2 )	  ∇ρ i 0	2 +	J	∇ρ i j,a	 j,s ω a 2 
												j=1 a∈Vj-s
				(7.5) ≤	C 2 max (w 1 , w 2 ) C 2 min (w 1 , w 2 )	(f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) = 0.	(7.20)
	η i alg =	(f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) ∇ρ i J,alg	(7.2) (7.7) ≥	C 2 C max (w 1 ) min (w 1 , w 2 )	 ∇ρ i 0 	2 +	j=1 a∈Vj-s J	∇ρ i j,a	j,s  ω a 2 1  2
		(7.18) ≥	C 2 min (w 1 ∇ρ i J,alg	(3.2) = β ∇(u J -u i J )	(7.21)
	for
		12	√	1 2C SMD J	5 2 (d + 1)	5 2	≤ β :=	C 2 min (w 1 ≤	8	√	1 d + 1	.

The bounds on β follow from (7.3),

(7.6)

, and

(7.19)

.

The codes were prepared to benefit as much as possible from Matlab's fast operations on matrices and vectors. The experiments were run on one Dell C6220 dual-Xeon E5-2650 node of Inria Sophia Antipolis -Méditerranée "NEF" computation cluster, in a sequential Matlab script.

Conclusions and outlookIn this work, we presented a hierarchical construction of the algebraic residual lifting in the spirit of Papež et al.[START_REF] Papež | Sharp algebraic and total a posteriori error bounds for h and p finite elements via a multilevel approach[END_REF]. This lifting approximates the algebraic error by one iteration of a V-cycle multigrid with no presmoothing step, a single damped additive Schwarz postsmoothing step, and a coarse solve of the lowest polynomial degree. The lifting leads us to an a posteriori estimator on the algebraic error and to a linear iterative solver. We showed that the two following results are equivalent: the (reliable) a posteriori estimator is p-robustly efficient, and the solver contracts p-robustly the error at each iteration. The provided numerical tests agree with these theoretical findings. Moreover, we also presented numerical results for a modified solver corresponding to a weighted restricted additive Schwarz smoothing. In accordance with the literature, this modified solver provides a further speed-up compared to the damped Schwarz smoothing. Although we currently cannot show that our p-robust theoretical result also applies to this construction, the use of high-degree polynomials does not seem to cause a degradation of the solver. So far, our theory involves estimates depending algebraically on the number of mesh levels J, which we do not observe in the numerical results for the weighted restricted variant. In forthcoming works, we plan to develop adaptivity based on the property (5.4), i.e., a computable splitting equivalent to the error and localized not only levelwise but also patchwise. Applications to more involved problems are also on our work list.

We return to the main estimate and obtain the result by using definition (7.6) (f, ρ i J,alg )-(∇u i J , ∇ρ i J,alg )≥

* This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR). The work of Jan Papež was supported by the NLAFET project as part of European Union's Horizon 2020 research and innovation program under grant 671633. The authors are grateful to Inria Sophia Antipolis -Méditerranée "NEF" computation cluster for providing resources and support.

Proofs of the main results

As shown in Corollary 5.4, the results of Theorems 5.1 and 5.3 are equivalent. Therefore it suffices to prove the first one. Our approach to proving Theorem 5.1 consists in studying the uncomputable exact residual lifting ρi J,alg given by (3.1) and its approximation ρ i J,alg given by Definition 3.2. In particular, we will estimate p-robustly the quantities ∇ρ i J,alg , ∇ρ i J,alg , and (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) by local contributions ρ i j,a of (3.8) used to construct ρ i J,alg , and we show that

when ρ i J,alg = 0, ρi J,alg = 0 when ρ i J,alg = 0, which also establishes Corollary 5.5.

Upper bound on ∇ρ i J,alg

We present here properties of the constructed residual lifting ρ i J,alg and its levelwise components ρ i j , where 1 ≤ j ≤ J.

7.2 Lower bound on (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg )

While studying the term (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ), the interaction of different level contributions ρ i j of the lifting ρ i J,alg arises naturally. In order to estimate these terms, the damping parameters w 1 , w 2 used in the construction (3.7) of our lifting prove to be essential. Lemma 7.2 (Estimate on (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) from below by patchwise contributions). Let ρ i J,alg be given by Definition 3.2. Then

where

Proof. We begin by using the construction of ρ i J,alg given in Definition 3.2 to write (f, ρ i J,alg ) -(∇u i J , ∇ρ i J,alg ) = (f, ρ i 0 ) -(∇u i J , ∇ρ i 0 ) +