Guaranteed contraction of adaptive inexact hp -refinement strategies with realistic stopping criteria

,


Introduction
The adaptive finite element method has been used in practice and theoretically studied for several decades.Its recent development has been catalysed by the extensive research dedicated to efficient and reliable a posteriori error estimates, see, e.g., the pioneering works of Babuška and Miller [3] and Dörfler [29] and the survey books by Ainsworth and Oden [1] and Verfürth [56].In the h-adaptive strategy for elliptic problems, Morin et al. [43] provided a plain convergence result, whereas Binev et al. [10] modified the method from [43] to prove not only optimal convergence rate but also optimal computational complexity.Other important results are to be found in Morin et al. [44,45], Cascón et al. [21], Stevenson [54,55], and Carstensen et al. [19].Most of the convergence results were stated for methods driven by residual-type a posteriori error estimates; the works addressing other types of estimators include Kreuzer and Siebert [39] and Cascón and Nochetto [22].In contrast, convergence of hp-adaptive strategies has been addressed only recently in Dörfler and Heuveline [30], Bürg and Dörfler [15], and Bank, Parsania, and Sauter [6].To our knowledge, the state-of-the-art optimality result is by Canuto et al. [16], hinging on an coarsening module due to Binev [9]; a polynomial-degree-robust (p-robust) contraction has been shown in Canuto et al. [17], based on a local saturation assumption.
In a common adaptive finite element method, a (larger and larger) system of algebraic equations needs to be solved at each step.This can be practically quite costly and is actually not necessary upon employing the idea to balance the algebraic error to the level of the discretization one at each step of the adaptive procedure, as suggested in, e.g., [7,10,54,55,8,2,33,50] and the references therein for linear elliptic problems and Holst et al. [37], Carstensen et al. [19], and Gantner et al. [36] for nonlinear elliptic problems.In the seminal works, Binev et al. [10] and Stevenson [54,55] present their h-adaptive finite element algorithms in an abstract setting and the inexact approximations are assumed to be sufficiently/arbitrarily close to the exact ones, without, however, a practical estimate on the algebraic error.The algebraic error can be estimated for a specific solver under appropriate assumptions as in Arioli et al. [2], where a convergence of the inexact h-adaptive finite element algorithm employing the conjugate gradient method as the algebraic solver has been shown, provided a good estimate on the smallest eigenvalue of the finite element system matrix is available.Another approach is based on an a priori argument which requires that the employed iterative solver contracts the algebraic error at least by a factor ρ it < 1 on each iteration.In particular, in such a situation, Becker and Mao [8] design a convergent and quasi-optimal conforming h-adaptive finite element algorithm.We believe, though, that an optimal way to proceed is via a dedicated a posteriori error estimate on the algebraic error.This does not rely on a specific solver or assumptions and may be much more precise, leading to quite sharp identification of all total, algebraic, and discretization error components, see [38,47,46] and the references therein.
The goal of this work, stemming from Chapter 3 of the Ph.D. dissertation [24], is to complete our recently proposed hp-adaptive refinement strategy with computable guaranteed bound on the error reduction factor [25] and its counterpart with an inexact algebraic solver setting [26] by a rigorous convergence proof.Similarly to [25,26], we examine the pure diffusion equation with homogeneous Dirichlet boundary conditions.Let Ω ⊂ R d , d = 2, 3, be a polygonal/polyhedral domain (open, bounded, and connected set) with a Lipschitz boundary ∂Ω, and let H 1 0 (Ω) denote the Sobolev space of all functions in L 2 (Ω) which have all their first-order weak derivatives in L 2 (Ω) and a vanishing trace on ∂Ω.Let f ∈ L 2 (Ω) be a given source term and let A ∈ [L ∞ (Ω)] d×d be a given diffusion tensor that we additionally suppose symmetric and positive definite; at some occasions below, we will additionally suppose it piecewise constant with respect to the computational mesh.The model problem in its weak form reads: find u ∈ H 1 0 (Ω) such that (A∇u, ∇v) = (f, v) ∀v ∈ H 1 0 (Ω), (1.1) where (•, •) stands for the L 2 (Ω) or L 2 (Ω) d inner product.The conforming hp-finite element method is used to discretize the model problem (1.1).Therein, we consider only matching simplicial meshes without hanging nodes.
Scheme 1: Paradigm of an hp-adaptive algorithm employing an exact algebraic solver.
The first part of this manuscript (Sections 2-6) is dedicated to the study of the hp-adaptive algorithm in the exact solver setting, i.e. we assume exact (up to machine precision) solution of all the resulting linear algebraic problems.Typically, such algorithms follow the well-established paradigm presented in Scheme 1. Particularly, in [25], we showed that on step of the adaptive loop from Scheme 1, it is possible to compute explicitly a real number 0 ≤ C ,red ≤ 1 such that where u ex and u ex +1 denote the exact finite element solutions from the respective iterations and + 1 of the adaptive loop.Even though we provided numerical evidence in [25, Section 6], we did not prove that the reduction factor C ,red of [25, Section 5] is bounded by a generic constant strictly smaller than one; this is the subject of our study here.In order to achieve it, in contrast to [25], we introduce two slight modifications of the modules MARK and REFINE of Scheme 1 in Section 4. First, we extend the marked region in module MARK by one layer of elements, see Section 4.2.And second, we either introduce an additional computational check on discrete stability of the local equilibrated fluxes in each patch marked for refinement within module REFINE, see condition (4.14) in Algorithm 1, or we use stronger local h-and p-refinements.
For a piecewise polynomial right-hand side f , a piecewise constant diffusion coefficient A, and when the maximal polynomial degree is bounded by some p max , we prove that estimate (5.16), with C lb possibly depending on p max , is satisfied when each patch marked for refinement is h-refined such that element and face interior nodes arise and/or p-refined by essentially a polynomial increase by d − 1.We will call this below the "interior node, p + d − 1" strategy.In the general case, we prove that the check (4.14) will be satisfied with a constant C lb independent of p upon strong-enough h-and/or p-refinements in each marked patch, where we propose to iterate one newest-vertex bisection or polynomial increase by one.We call this "adaptive residual lifting" strategy below.The requested bound on the maximal polynomial degree in the first result above is restrictive from the theoretical viewpoint (it only gives a h-convergence result with prefinements included but limited), though reasonable in practice.The second result then gives both a h-and p-robust contraction (1.2), but the numbers of h-and/or p-refinement iterations in each marked patch are not necessarily limited; in our numerical experiments, though, only one such iteration was always sufficient, just as in [25,26].The key ingredient of the proof is the stability of the local equilibrated fluxes: non p-robust via the bubble function technique in the fist case and p-robust via polynomial extension operators in the second case.We summarize these developments in Section 5 and state the guaranteed contraction result (1.2) in details in Section 6.
Scheme 2: Paradigm of an adaptive loop employing an inexact algebraic solver.
In the second part of this manuscript (starting with Section 7), we extend our results concerning the convergence (guaranteed contraction) of the hp-adaptive algorithm from the exact solver setting to the inexact solver setting.Following [26], we incorporate the use of an arbitrary algebraic solver within the framework of adaptive algorithm of Scheme 1 by replacing the module SOLVE by an adaptive sub-loop consisting of modules ONE SOLVER STEP and ESTIMATE, see Scheme 2. In contrast to [26], we need to increase the equilibration polynomial degree by one.We treat carefully our choice of a practical adaptive stopping criterion for the algebraic solver in order to ensure the convergence of the adaptive algorithm.In contrast to [10,54,55,8,2], the particular advantage of the present work is that no request on the algebraic solver is made and our bounds on algebraic, discretization, and total errors do not contain any generic constant and remain fully computable, see Section 7.1.Consequently, the values of the stopping parameter γ , which expresses the percentage of the algebraic error with respect to the total error, do not need to theoretically take excessively small (and unknown) values.The proper choice of stopping criterion together with the modifications of modules MARK and REFINE introduced in the exact solver setting then allow us to show, in extension of property (1.2), the error reduction between two consecutive inexact approximations u and u +1 produced by the adaptive loop of Scheme 2 with a fully computable factor 0 ≤ C ,red ≤ C < 1.This theoretical result improves in particular the developments of [26,Theorem 5.4], where the reduction factor was not shown to be bounded by a generic constant strictly smaller than one.In particular, the reduction property (1.3), showed in Section 9 here, implies the convergence of the adaptive algorithm prescribed by Scheme 2. The quick numerical assessment of Section 10, which complements the extensive numerical studies from [25] and [26], illustrates that the generated sequences of meshes and polynomial degree distributions still lead to asymptotic exponential convergence, and this for the stopping parameter γ with the very-reasonable-in-practice value around 0.05 in the inexact solver case.Moreover, condition (4.14) in Algorithm 1, with C lb independent of p which ensures p-robust contraction, was here satisfied already for one newest-vertex bisection step and/or increase of the polynomial degree by one.

Road map
To help the reader orientate in our paper, we provide here a quick "road map" of our developments to come; let temporarily the source term f be a piecewise polynomial and the diffusion coefficient A be piecewise constant, so that no "data oscillations" arise, and suppose an exact algebraic solver.Let u be the (unknown) weak solution from (1.1), let u ex be the (known) finite element approximation on mesh T and space V , and let u ex +1 be the (unknown) finite element approximation on the forthcoming mesh T +1 and space V +1 (containing V ) that has been designed by our hp-refinement strategy on the current step .The Galerkin orthogonality of u ex +1 yields In our approach, we use an equilibrated flux reconstruction σ a on each patch subdomain ω a associated with a vertex a ∈ V ; σ a is a Raviart-Thomas piecewise polynomial with no flow through the boundary of ω a on interior vertices.We then set σ = a∈V σ a , which gives the guaranteed error upper bound We use A 1 2 ∇u ex + A − 1 2 σ K for marking of mesh elements K to h-or p-refinement; more precisely, the set of marked elements M θ has to contain a θ-fraction of the overall estimated error via the bulk-chasing (Dörfler) criterion Moreover, on each patch ω a around a vertex a ∈ V marked for refinement, we construct a residual lifting r a,hp , which is a conforming piecewise polynomial taking zero value on the boundary of ω a .Then we have the guaranteed discrete lower bound Now, we can already bound the error A 1 2 ∇(u−u ex +1 ) , with respect to A 1 2 ∇(u−u ex ) , prior to the computation of u ex +1 .Indeed, using respectively the guaranteed discrete lower bound (2.4), Dörfler's marking (2.3), and the guaranteed upper bound (2.2), we have for the computable constant By (2.6), we can monitor whether a contraction in the finite element approximation is happening, but an important question is whether C ,red is uniformly bounded by a constant C red strictly smaller than one.This is what we establish here.A crucial ingredient for this is that we design the above residual liftings r a,hp such that the local stability holds for a generic constant C lb .Then, overlaps between patches lead to and patches overlaps again, together with the local stabilities (2.7), imply the global stability on the refined elements collected in M θ in the form (2.9) From (2.9), we conclude that indeed

Framework and notation
This section details the notation employed throughout the manuscript.Within the adaptive loops of Schemes 1 and 2, a sequence {(T , p )} ≥0 , with ≥ 0 the iteration counter, is generated.Each pair (T , p ) consists of a matching simplicial mesh T of the computational domain Ω, i.e., a finite collection of (closed) simplices K ∈ T covering Ω and such that the intersection of two different simplices is either empty or their d -dimensional common face, 0 ≤ d ≤ d−1, and of a polynomial-degree distribution vector p := {p ,K } K∈T which assigns a degree p ,K ∈ N ≥1 to each simplex K ∈ T .Moreover, each pair (T , p ) prescribes a discrete finite-dimensional space V , defined as where P p (T ) denotes the space of piecewise polynomials of total degree at most p ,K on each simplex K ∈ T .Let us denote by N the dimension of the -th level space V .Note that we enforce the H 1 0 (Ω)conformity of the spaces V for all ≥ 0. In addition, we make the following nestedness assumption: The initial pair (T 0 , p 0 ) is assumed to be given.Then, the purpose of each step ≥ 0 of the adaptive loops of Schemes 1 and 2 is to determine the next pair (T +1 , p +1 ).The nestedness property (3.1) gives us two crucial restrictions on the meshes and polynomial-degree distributions defining the spaces V : (i) the sequence of meshes {T } ≥0 needs to be hierarchically nested, i.e., for all ≥ 1 the mesh T is a refinement of T −1 such that for all K ∈ T , there is a unique simplex K ∈ T −1 , called the parent of K, satisfying K ⊆ K; (ii) the local polynomial degree is locally increasing, i.e., for all ≥ 1 and all K ∈ T , p ,K ≥ p −1, K , where K ∈ T −1 is the parent of K.Moreover, we assume the following standard shape-regularity property: there exists a constant κ T > 0 such that max K∈T h K /ρ K ≤ κ T for all ≥ 0, where h K is the diameter of K and ρ K is the diameter of the largest ball inscribed in K.This can in practice be ensured by using e.g. the newest-vertex bisection mesh refinement algorithm [53,41,55].
We denote by V (respectively F ) the set of vertices (respectively (d − 1)-dimensional faces) of the mesh T .These are decomposed into interior vertices V int ((d − 1)-dimensional faces F int ) and boundary vertices V ext ((d − 1)-dimensional faces F ext ).For each vertex a ∈ V , ≥ 0, the so-called hat function ψ a is the continuous, piecewise affine function that takes the value 1 at the vertex a and the value 0 at all the other vertices of V ; the function ψ a is contained in the space V for all vertices a ∈ V int and all polynomial degrees p .Furthermore, we consider the simplex patch T a ⊂ T which is the collection of the simplices sharing the vertex a ∈ V , with ω a the corresponding open subdomain of Ω, coinciding with the support of ψ a .Let F a ⊂ F denote the set of all the (d − 1)-dimensional faces in the patch T a .This set can be further decomposed into F a,int , the subset of faces from F a that share the vertex a and are shared by two distinct simplices in T a , and F a,ext , the faces from F a lying in ∂ω a , so that F a = F a,int ∪ F a,ext .For any face F ∈ F , n F stands for a unit vector normal to F with an arbitrary but fixed orientation.The operator [[•]] yields the jump, in the direction of n F , of the traces of the argument from the two simplices that share F ∈ F int , and the actual trace for F ∈ F ext .Finally, for each simplex K ∈ T , V ,K denotes the set of vertices of K and F ,K denotes the set of (d − 1)-dimensional faces of element K ∈ T .

The hp-adaptive algorithm with exact solver
In this section we first recall the modules SOLVE and ESTIMATE as defined in [25].Afterwards, we introduce the slightly modified versions of the remaining modules MARK and REFINE from the adaptive loop of Scheme 1, which will allow us to prove the boundedness of the contraction factor and thus the convergence of such an adaptive algorithm.

The modules SOLVE and ESTIMATE
Let ≥ 0 denote the current iteration counter.The module SOLVE takes as input the current finite element space V ⊂ H 1 0 (Ω) and outputs the Galerkin approximation u ex ∈ V of the weak solution u of (1.1) defined as the unique solution of (A∇u ex , ∇v Note that obtaining u ex is equivalent to computing the exact solution of the system of linear algebraic equations where {ψ n } 1≤n≤N is a basis of the -th level space V , so that Following [27,13,34,28], see also the references therein, the module ESTIMATE relies on an equilibrated flux a posteriori error estimate on the energy error A The equilibrated flux is constructed locally on the simplex patches T a attached to each vertex a ∈ V .For this construction, we consider as in [28,25] the local polynomial degree p est a := max K∈T a p ,K (any other choice so that p est a ≥ max K∈T a p ,K can also be employed).For a fixed vertex a ∈ V , let the broken space where RTN p (K) := [P p (K)] d + P p (K)x is the usual p-th order Raviart-Thomas-Nédélec space (cf.[14,51]) on a simplex K ∈ T .Then, the patchwise normal-trace-continuous spaces V a with homogeneous Neumann boundary conditions in which the local equilibration will be performed are defined by 3) Let P p est a (T a ) be composed of elementwise polynomials of degree p est a on the patch T a .We will also need the L 2 (ω a )-orthogonal projection onto P p est a (T a ) that we denote as Π a .Definition 4.1 (Equilibrated flux σ by local minimizations).Let u ex be the solution of (4.1).For each vertex a ∈ V , let the local equilibrated flux σ a ∈ V a be defined by the local minimization problem Then, extending each local contribution σ a by zero outside of the patch domain ω a , the global H(div, Ω)conforming equilibrated flux σ is constructed as It will turn useful to denote Then, the Neumann compatibility condition ω a g a = ω a ψ a f − ∇ψ a •A∇u ex = 0 for the problem (4.4a) is satisfied for all a ∈ V int as a direct consequence of (4.1) (consider ψ a as a test function in (4.1)).The global H(div, Ω)-conformity of σ follows from imposing the zero normal trace of functions in the local spaces V a , cf. the definition (4.3).
Then, as stated in [27,13], see also [28,Theorem 3.3] and [23,Theorem 4.5], the following guaranteed upper bound on the energy error holds true: ) where c A,K is the smallest eigenvalue of the diffusion tensor A on the simplex K.When f is a piecewise polynomial in the space P p −1 (T ) and when A is piecewise constant, A ∈ P 0 (T ) d×d , the so-called data A,K ) vanishes in the local estimator η K .This follows from (4.8) below.It would also vanish if f ∈ P p (T ) for uniform polynomial degrees p = const, see (6.8) and the discussion in Remark 6.7 below.

The module MARK
The module MARK takes as input the local error estimators from (4.6) and proceeds in two phases.
The first phase corresponds to the module MARK used in [25].We select the smallest subset of marked vertices V θ ⊂ V using a bulk-chasing criterion inspired by the well-known Dörfler's criterion (cf.[29]) where θ ∈ (0, 1] is a fixed threshold parameter, and where, for a subset S ⊂ T , we adopt the notation η(u ex , S) := K∈S η K (u ex ) 2 1/2 .Then, letting T a ⊂ T be the collection of all the simplices that belong to a patch associated with a marked vertex, we observe that (4.7) means that η(u ex , M θ ) ≥ θ η(u ex , T ).Let us denote by ω := a∈V θ ω a the open subdomain corresponding to the set of marked elements M θ .To select a set V θ of minimal cardinality, the mesh vertices in V are sorted by comparing the vertex-based error estimators η(u , T a ) for all a ∈ V , and a greedy algorithm is employed to build the set V θ .Possible reductions of the computational cost of such an algorithm are proposed in [29, Section 5.2] and in Stevenson [55, Section 5], whereas a construction of the minimal set at a linear cost is identified in Pfeiler and Praetorius [48].
In a second phase, we define an extended set of marked vertices V with the corresponding set of marked elements M in the following way In other words, we extend the set M θ by one more layer of neighbouring elements in contact with the boundary ∂ω , see Figure 1 for an illustration.In addition, we define by ω := a∈V ω a the open subdomain corresponding to the set of marked elements M .This extension is motivated by the structure of the error estimate (4.6), stemming from equilibrated flux σ composed of local patchwise contributions σ a .It plays a particular role later in our convergence proof, when we decompose the error estimate η(u ex , M θ ) into a sum of local patchwise contributions in (6.6) with their corresponding domains of definitions possibly exceeding the original marked subdomain ω , but always included in the extended marked subdomain ω .The present theoretical analysis requires such extension of the marked region; this is the price we pay for the precision of our estimates on the error reduction factor C ,red : Remark 4.2 (Marking and analysis without the additional layer).Instead of the elementwise estimate (4.6), one could employ the following patchwise form of the error estimate, cf.[34,Lemma 3.22], [17, since Π a projects to P p est a (T a ), which contains the polynomials of order p ,K in each mesh element K ∈ T a , and recalling (4.5), the term ∇ψ a •A∇u ex − Π a (∇ψ a •A∇u ex ) vanishes if the diffusion tensor is piecewise constant, A ∈ P 0 (T ) d×d , whereas the term ψ a f − Π a (ψ a f ) vanishes if the source term is a piecewise polynomial, f ∈ P p −1 (T ).This would result in no need to extend the obtained set of marked vertices V θ by the additional layer, as well as in an overall a simpler hp-adaptive algorithm based on the patchwise estimators η a (u ex ).However, unless a combination with (4.6) was used, the presence of the constant √ d + 1 already in the error estimate (4.8) would lead to a deterioration of the bound on the error reduction factor C ,red , which was designed as sharp as possible in [25,26].

The module REFINE
The module REFINE takes as input the extended set of marked vertices V and outputs the mesh T +1 and the polynomial-degree distribution p +1 to be used at the next iteration of the adaptive loop from Scheme 1.
An example of a marked patch T a together with its polynomial-degree distribution p a (left), its hrefinement (center ), and p-refinement (right), the "adaptive residual lifting" strategy of Section 4.3.2.

h and p patch refinements and the hp-decision
The idea employed in [25,26] is to emulate separately the effects of h-and p-refinement on each patch T a assigned to a marked vertex a ∈ V using two distinct local primal solves.For this purpose, we will need a matching simplicial refinement T a,h of T a (see the center panel of Figure 2 for an example), as well as an increased polynomial-degree distribution vector p a,p (see the right panel of Figure 2 for an example); let also p a,h be the (original) polynomial-degree distribution in the patch T a given by the vector p a := {p ,K } K∈T a and T a,p = T a .For each vertex a from the extended set of marked vertices V , consider the local patch-based spaces V a,h and V a,p given by In the spirit of Babuška and Miller [3], and similarly to the element-patch problems in Bürg and Dörfler [15, Section 3.2.1],we let r a,h ∈ V a,h and r a,p ∈ V a,p be the h-and p-refinement residual liftings respectively given by the two above local primal problems use Dirichlet boundary conditions, while we refer to [26] for the possible use of Neumann boundary conditions.The hp-decision is then made based on the comparison of "attainable error decreases" A 1 2 ∇r a,h ω a and A 1 2 ∇r a,p ω a and is described on line 9 of Algorithm 1, also relying on the details of the two following refinement strategies.We would like to confess immediately that, actually, this hp-decision is currently not exploited in the proof of the contraction property (6.1) of Theorem 6.1, see Remark 6.6 below for a discussion.

"Adaptive residual lifting" refinement strategy; hp residual lifting
For each vertex a from the extended set of marked vertices V , this strategy corresponds to iteratively applying the original strategy from [25,26].Let C lb be a positive constant; theoretically, this should be taken as (5.9) from Section 5.1, whereas in practice, we use C lb = 10.
For h-refinement testing, we define T a,h by diving each simplex K ∈ T a into at least two children simplices simply by applying the selected mesh refinement algorithm, e.g., the newest-vertex bisection [53,41,55], see Figure 2, center.The polynomial-degree distribution p a,h corresponding to T a,h is then obtained from p a by assigning to each newly-created simplex the same polynomial degree as that of its parent.For p-refinement testing, we set T a,p = T a and the polynomial-degree distribution p a,p is obtained from p a by assigning to each simplex K ∈ T a,p the polynomial degree p ,K + δ a K where, as in [25,26], This is illustrated in Figure 2, right, and only increases the smallest polynomial degree in the patch T a .In difference with respect to [25,26], though, we also immediately consider the suggested new pair (T +1 , p +1 ) identified by the hp-decision on all vertices a ∈ V after step 20 of Algorithm 1, leading to the suggested finite element space V +1 := P p +1 (T +1 ) ∩ H 1 0 (Ω).Then we consider local subspaces of V +1 , still on ω a , given by and construct the patch residual lifting r a,hp ∈ V a,hp by solving Input: extended set of marked vertices V

3:
Output: new pair (T +1 , p +1 ) 4: for all a ∈ V do 5: Prepare the patch h-refinement space V a,h and the patch p-refinement space V a,p from (4.9) in depen-

6:
dence on the chosen strategy of either Section 4.3.2 or Section 4.3.3

14:
for all K ∈ T +1 with its parent element K ∈ T do 15: if "adaptive residual lifting" strategy of Section 4.3.2then 16: else if "interior node, p + d − 1" strategy of Section 4.3.3then for all a ∈ V do

23:
Compute the hp residual lifting r a,hp given by (4.13).

24:
if the condition is not satisfied then 25: re-assign T := T +1 and p := p +1 and go back to line 4.

27:
end for 28: end if 29: end module Now, the idea is to check the condition (4.14), similar to the "local saturation property" (iii) from [17,Section 4].If it is satisfied (which was always the case in the numerical experiments in Section 10 below), then we accept the suggested new pair (T +1 , p +1 ).If (4.14) is not satisfied, we continue working on ω a by further diving each simplex K ∈ T a to create a more refined T a,h and by increasing p a,p , until (4.14) is satisfied (which will be reached by virtue of Proposition 5.1 below, at least for small data oscillation).

"Interior node, p + d − 1" hp refinement strategy
We will see below that the strategy from Section 4.3.2 is p-robust, but it may lead to an iteration loop for finding (T +1 , p +1 ).Alternatively, when the polynomial degrees in p can be supposed bounded by p max (which is restrictive from the theoretical point of view, though acceptable in practice), we may identify a sufficient mesh refinement and polynomial degree increase in each marked patch ω a , a ∈ V , to automatically satisfy a condition of the form (4.14), without any iteration, as follows.Let for simplicity the datum f be a piecewise polynomial of variable degree at most p f with respect to the coarsest partition T 0 , f ∈ P p f (T 0 ) for p f ≤ p 0 − 1, as well as A ∈ P 0 (T 0 ) d×d piecewise constant.
For each marked vertex a ∈ V , we let T a,h be obtained from T a by dividing each simplex K ∈ T a into at least six children simplices such that a node interior to all the elements K ∈ T a as well as all the faces F ∈ F a,int is created, (4.15) cf. the requirement of having an interior node and its implementation in two and three dimensions in [43].This requirement can actually be relaxed, as done previously in, e.g., [31,Lemma 11].In particular, when the space dimension d = 2 and either f = 0 or p f,K ≤ p ,K − 2 for all K ∈ T a , then only the faces from F a,int need to have an interior node but not the elements from T a .The polynomial-degree distribution p a,h corresponding to T a,h is then copied from p a , with a possible increase to match the polynomial degree of its marked neighbors.More precisely, for each K ∈ T a , we set p ,K as the maximum of the polynomial degrees p ,K and all p ,K where K ∈ T a shares a face with K.An illustration can be found in the right panel of Figure 3. Next, for p-refinement testing on each marked patch, we let T a,p := T a and define the polynomial-degree distribution vector p a,p by assigning to each simplex K ∈ T a,p the polynomial degree with p f,K the local polynomial degree of f on K, see line 18 of Algorithm 1 and Figure 3, right, for illustration.
In this strategy, r a,hp from (4.13) is only to be computed at a later stage in Section 4.3.4below.Congruently, there is no need to check condition (4.14) on line 24 of Algorithm 1, since it is automatically satisfied in view of (5.16) from Proposition 5.2 below, though possibly with a non p-robust constant C lb .Figure 3: An example of a marked patch T a together with its polynomial-degree distribution p a (left), its hrefinement (center ), and p-refinement (right), the "interior node, p + d − 1" strategy of Section 4.3.3.

Computable discrete lower bound
Once the new pair (T +1 , p +1 ) is determined using one of the above strategies, the finite element space V +1 := P p +1 (T +1 ) ∩ H 1 0 (Ω)) to be used on iteration ( + 1) of the adaptive loop of Scheme 1 is at our disposal.We now proceed by extending the discrete lower bound of [25,Lemma 5.1] to the present setting, which relies on constructing the hp residual lifting r a,hp by (4.13).After extending r a,hp from (4.13) by zero outside ω a , we have for the current (available) approximation u ex ∈ V and the next level's (unavailable) Moreover, this lower bound can be further localized using the fact that each simplex has (d + 1) vertices as Indeed, this can be seen from 5 Discrete stability of equilibrated fluxes in an exact solver setting We prove in this section that condition (4.14) of Algorithm 1, requested for the "adaptive residual lifting" strategy of Section 4.3.2,can be satisfied.Section 5.1 does it in a p-robust way, for a sufficient h-and/or p-refinement of each marked patch ω a .In Section 5.2, we then show that this condition automatically holds for the "interior node, p + d − 1" strategy of Section 4.3.3 with a fixed h-and/or p-refinement, albeit in a non-p-robust way.

p-robust stability for the "adaptive residual lifting" strategy
We prove here condition (4.14) for the "adaptive residual lifting" strategy of Section 4.3.2, with the constant C lb independent of the polynomial degree p.Let us consider the infinite-dimensional local H(div) space with the normal trace understood as usual by duality, and the infinite-dimensional constrained minimization for a constant C st that only depends on the space dimension d, the mesh shape-regularity κ T , and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient A on ω a .Computable upper bounds on C st are discussed in [34,Lemma 3.23 and Remark 3.24].
Let also (5.4) where the constant C cont,PF only depends on the space dimension d, the mesh shape-regularity κ T , and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient A on ω a ; C cont,PF relies on the Poincaré-Friedrichs inequality and is fully computable.Finally, denote η osc,a := where we recall that g a is given by (4.5) and that c A,K is the smallest eigenvalue of the diffusion tensor A on the simplex K. Also recall from Remark 4.2 that ∇ψ a •A∇u ex − Π a (∇ψ a •A∇u ex ) vanishes if the diffusion tensor is piecewise constant, A ∈ P 0 (T a ) d×d , whereas the term ψ a f − Π a (ψ a f ) vanishes if the source term is a piecewise polynomial, f ∈ P p −1 (T a ).
and let the local equilibrated flux σ a be constructed by (4.4a).Let data oscillation do not dominate in that Then, there exists a constant C lb ≥ 1 only depending on the space dimension d, the mesh shape-regularity κ T , and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient A so that, for a sufficiently fine patch space V a,hp from (4.12) and r a,hp given by (4.13), there holds More precisely, it is possible to take where the constants C st and C cont,PF respectively stem from (5.3) and (5.5).
Proof.We proceed in several steps.
Step 1. Remark that the left-hand side of (5.8) can be bounded by the triangle inequality as (5.10) Step 2. In addition to σ a , let us denote by σ a the solution to problem (5.2) where the polynomial projection Π a (ψ a f −∇ψ a •A∇u ex )= Π a (g a ) is replaced by ψ a f −∇ψ a •A∇u ex = g a .Treating the mismatch between these two terms as in [34, proof of Theorem 3.17], one uncovers from (5.3) that Moreover, the primal-dual equivalence as in [34,Remark 3.15] or [35,Corollary 3.6] shows that (5.12) Step 3. From (5.5), we can rewrite and bound the right-hand side of (5.12) as max where r a ∈ H 1 0 (ω a ) is given by (5.14) Step 4. From the above steps, we have From here, condition (5.7) on η osc,a yields Thus, Step 5. Now, recall that the patch ω a and the data f and u ex are fixed, whereas r a and r a,hp are prescribed therefrom by respectively (5.14) and (4.13).Moreover, the space V a,hp from (4.12) for computing r a,hp (local to the patch ω a and discrete) is created by applying h-and/or p-refinements in the patch ω a .Thus, there holds A and there exists a sufficiently fine patch space V a,hp such that A ∇r a ω a /2, so that (5.8) holds with the constant C lb given by (5.9).

Non p-robust stability for the "interior node, p + d − 1" strategy
We prove here condition of the form (4.14) for the "interior node, p + d − 1" strategy of Section 4.3.3,with the constant C lb possibly dependent on the polynomial degree p.We suppose here f piecewise (p 0 − 1)polynomial and A piecewise constant; otherwise the bubble function construction below would require additional polynomial degree enrichment.Proposition 5.2 (Discrete local non p-robust stability of the flux equilibration).Let f ∈ P p f (T 0 ) fulfill p f ≤ p 0 − 1 and let the diffusion tensor A be piecewise constant with respect to T 0 , A ∈ P 0 (T 0 ) d×d .Let u ex ∈ V satisfy the hat function orthogonality and let the local equilibrated flux σ a be constructed by (4.4a).Finally, let the patch space V a,hp from (4.12) involve the interior nodes (4.15) or the polynomial degree increase (4.16) and let r a,hp be given by (4.13).
Then there exists a constant C lb ≥ 1 only depending on the space dimension d, the mesh shape-regularity κ T , the ratio of the largest and the smallest eigenvalue of the diffusion coefficient A, and the maximal polynomial degree p max such that To prove Proposition 5.2, we will rely on two auxiliary results employing the bubble function technique, cf.[56].We present the forthcoming developments with A being an identity matrix for simplicity; extension to piecewise constant A is immediate.From now on, we use the shorthand notation x 1 x 2 when there exists a generic positive constant C that only depends on the space dimension d, the shape-regularity κ T of the underlying hierarchy of meshes, and the polynomial degree of approximation employed locally such that x 1 ≤ Cx 2 .
Lemma 5.3 (Discrete stability of the element residuals).Let f ∈ P p f (T 0 ) with p f ≤ p 0 − 1 and A = I.Let a vertex a ∈ V , the patch space V a,hp from (4.12), involving the interior nodes (4.15) or the polynomial degree increase (4.16) of Section 4.3.3,and the corresponding r a,hp from (4.13) be given.Then (5.17) Proof.Fix an element K ∈ T a and set We note that v K ∈ P max{p ,K −2,p f,K } (K), using the assumption f ∈ P p f (T 0 ).Let us define ψ K , a bubble function on element K, depending on which refinement has been applied within the REFINE module in Section 4.3.3 for the vertex a, to be a function more precisely, in the first case, ψ K is a piecewise affine function taking nonzero values in the(/all) interior nodes in K, whereas in the second case, ψ K is the usual element bubble function.We crucially remark that ψ K v K extended by zero is contained in the local hp-refinement space V a,hp due to (4.12) with either (4.15) or (4.16).In particular, either h-refinement has been employed and then T +1 | K contains an interior node, or p-refinement has been applied and then the polynomial degree has been increased to at least max{p ,K + d − 1, p f,K + d + 1}.Now, the equivalence of norms on finite-dimensional spaces gives where the hidden constant depends on the polynomial degrees.This in combination with (5.18) yields Then, noting that (∇u ex )| K ∈ H(div, K) and ψ K v K ∈ H 1 0 (K), employing the Green theorem and using (4.13) with A = I and the Cauchy-Schwarz inequality, we obtain (5.21) Then, using the inverse inequality (cf., e.g.[49, Proposition 6.3.2]),where again the hidden constant depends on the polynomial degrees, we have Moreover, from the definition of the bubble function ψ K , there holds Finally, (5.21) with (5.22) and (5.23) lead to the assertion of the lemma.
For the relaxation mentioned in Section 4.3.3,where for h-refinement in the case d = 2 only the faces from F a,int have an interior node and either f = 0 or p f,K ≤ p ,K − 2 for all K ∈ T a , one rather uses , the edge bubble function of a newly-added edge, in place of the first choice in (5.19).
Lemma 5.4 (Discrete stability of the face residuals).Let f ∈ P p f (T 0 ) with p f ≤ p 0 − 1 and A = I.Let a vertex a ∈ V , the patch space V a,hp from (4.12), involving the interior nodes (4.15) or the polynomial degree increase (4.16) of Section 4.3.3,and the corresponding r a,hp from (4.13) be given.Let T F denote the set of elements containing the two elements K ∈ T a that share a face F ∈ F a,int , and let ω F be the corresponding open subdomain.Then Proof.Fix an edge F ∈ F a,int .This time, we set and note that v F ∈ P max{p ,K −1,p ,K −1} , where K, K ∈ T a share the face F ; this is the motivation for the local polynomial degree increase in the REFINE module from Section 4.3.3.Let ψ F be the bubble function on T F , depending on whether h-or p-refinement was applied within the REFINE module in Section 4.3.3;namely, for the vertex a ∈ V , (5.25) similarly to (5.19), ψ K is respectively a piecewise affine function taking nonzero values in the(/all) interior nodes in F , or the usual face bubble function.Again, the face interior node property obtained from (4.15) or the polynomial degree increase obtained from (4.16) ensure that ψ F v F extended by zero lies in the local hp-refinement space V a,hp (4.12), which will be crucial below.By equivalence of norms on finite-dimensional spaces, similarly to (5.20), there holds, with a p-dependent constant, Let us keep the same notation for the extension of the function v F , with its original domain of definition being only the face F , to a function defined on the two simplices of T F .The extension is done by constant values in the direction of the barycenter of the face towards the vertex opposite to F .Then, there holds Employing (5.26), (5.24), and expanding the jump term with T F := {K, K } and the convention that n F points from K to K , we have We have also used the fact that (ψ F v F ) | ∂ω F = 0 due to the definition of the bubble function ψ F .Moreover, for each simplex K ∈ T F , (ψ F v F ) | K ∈ H 1 (K) and (∇u ex )| K ∈ H(div, K), so we are able to apply the Green theorem for both terms on the right-hand side of (5.28).Adding and subtracting (f, ψ F v F ) ω F , using that ψ F v F extended by zero lies in V a,hp , and recalling the discrete problem (4.13) with A = I, we finally obtain From the definition of the bubble function ψ F , there holds while the inverse inequality and the shape-regularity of the mesh T give, again with a p-dependent constant, Then, the following chain of inequalities holds true due to the Cauchy-Schwarz inequality, (5.29), (5.30), (5.31), and (5.27): Now, the assertion of the lemma follows from the above result combined with (5.17) and the shape-regularity of the mesh T .
Proof of Proposition 5.2.Let a ∈ V be fixed and recall we consider here A = I and f ∈ P p f (T 0 ) with p f ≤ p 0 − 1, so that Π a (g a ) =Π a (ψ a f −∇ψ a •A∇u ex ) = ψ a f −∇ψ a •A∇u ex = g a .We will crucially employ two reformulations of the local minimization problem (4.4a).Let the patchwise discontinuous piecewise polynomial spaces Q a be defined by (5.32) Then the Euler-Lagrange conditions for (4.4a), imposing the divergence constraint via a Lagrange multiplier, amount to solving the local Neumann mixed finite element problem: seek the pair (σ a , δ a ) ∈ V a × Q a such that (5.33b) Next, imposing the normal trace continuity constraint on faces from F a,int and the no-flux condition on all faces from F a,ext for a ∈ V int and only faces from F a,ext not sharing the vertex a ∈ V ext via a Lagrange multiplier leads to the hybridized formulation: seek σ a in the broken space RTN p est a (T a ), δ a ∈ Q a , and where we used the notation (5.35) We note that ψ a ∇u ex +σ a ∈ RTN p est a (T a ), as (4.3) uses the maximal local polynomial degree p est a .Thus, we can use it as a test function v in (5.34a).Afterwards, employing δ a as a test function q in (5.34b) and λ F as ξ in (5.34c), summing (5.34c) over all F ∈ F a \ F a,∂Ω , and finally summing (5.34a)-(5.34c),we obtain For the estimate on the right-hand side of (5.36), we have also used the fact that the hat function ψ a | F = 0 ∀F ∈ F a,ext for vertices a ∈ V int , so that the boundary faces are excluded from the sum, and ψ a | F = 0 ∀F ∈ F a,ext \ F a,∂Ω for a ∈ V ext .We now proceed by bounding the terms in the estimate of (5.36).Firstly, by the inf-sup stability of the mixed discretization (5.33), see e.g.[14] or [57, Theorem 5.9], we have δ a ω a h ω a ψ a ∇u ex + σ a ω a . (5.37) Secondly, ψ a ∞,K ≤ 1 together with Lemma 5.3 yield the estimate (5.38) Combining (5.36) with (5.38), the Cauchy-Schwarz inequality, the shape-regularity leading to h K ≈ h ω a , and (5.37), we obtain the estimate on the first term ∇r a,hp ω a ψ a ∇u ex + σ a ω a . (5.39) We continue by bounding the second term in the estimate of (5.36).By the characterization of the degrees of freedom in the Raviart-Thomas-Nédélec spaces, if F ∈ F a,int is a face of an element K ∈ T a , we have Fix v ∈ RTN p est a (K) with the constraints as in (5.40).Using this v as a test function in (5.34a) and the Cauchy-Schwarz inequality, we obtain (5.41) We now treat the two terms of (5.41) separately.For the first term, we start by employing the inverse inequality ∇•v K h −1 K v K .Furthermore, by scaling arguments, under the constraints of (5.40), we have (5.42) The first term of (5.41) can then be bounded as follows: where we have also used the bound (5.37) and the mesh shape-regularity yielding h ω a ≈ h K ≈ h F .On the other hand, (5.42) leads to the following bound on the second term of (5.41) In order to bound the supremum in (5.40), we combine (5.43), (5.44), and (5.41), to get F ψ a ∇u ex + σ a ω a . (5.45) Afterwards, ψ a F ≤ 1 and Lemma 5.4 yield ∇r a,hp ω F .

Proof of guaranteed contraction with an exact solver
Our first main result concerns the contraction of our hp adaptive algorithm with an exact algebraic solver: Theorem 6.1 (Guaranteed contraction, "adaptive residual lifting" strategy of Section 4.3.2).Let u be the weak solution of (1.1) and let u ex be the approximate solution of (4.1) on hp iteration ≥ 0. Let the pair (T +1 , p +1 ) be obtained by the module REFINE of Section 4.3 with the "adaptive residual lifting" strategy of Section 4.3.2, and let V +1 = P p +1 (T +1 ) ∩ H 1 0 (Ω) be the finite element space to be used on iteration ( + 1) of the adaptive algorithm prescribed by Scheme 1.Let also η M be the computable discrete lower bound defined by (4.17).Let finally the data oscillation do not dominate in that (5.7) holds.Then, two options arise.Either the a posteriori error estimate η(u ex , T ) from (4.6) is zero, in which case u ex = u, and the adaptive loop terminates.Or the new numerical solution so that C ,red is a fully computable bound on the error reduction factor and the generic constant C red only depends on the marking parameter θ, the space dimension d, the mesh shape-regularity κ T , and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient A; here, C lb is the constant from (5.8) that can theoretically be chosen as (5.9).
Proof.Let us assume that η(u ex , M θ ) = 0, the other case being trivial.The proof proceeds in two steps.
Step 1.We follow the proof of [25, Theorem 5.2], see also proofs of [43, Theorem 3.1] and [55, Theorem 5.3].To start with, from (1.1) and (4.1) on mesh level ( + 1), the Pythagorean relation reads By the computable lower bound on the incremental error (4.17) and the marking criterion (4.7), we have Hence, combining (6.3) with (6.4) and using the error estimate (4.6), we infer Step 2. We show that the reduction factor C ,red is indeed bounded by the positive constant C red strictly smaller than one, again defined in (6.2).We decompose the elementwise estimators from (4.6) using the partition of unity a∈V K ψ a | K = 1, the triangle inequality, the Cauchy-Schwarz inequality, and the fact that each simplex has (d + 1) vertices, as We note that in (6.6), we have crucially employed the extended set of marked vertices V of Section 4.2, as well as the divergence constraint from (4.4a) together with a∈V ,K ∇ψ a •A∇u ex = 0 and the notation (4.5).Next, (5.8) from Proposition 5.1 applied on each patchwise contribution of the sum in (6.6) and the localization (4.18) of η M given by (4.17) yield Finally, employing (6.7) to bound C ,red in (6.5) from above and defining the constant C red as in (6.2) finishes the proof.
Assuming the linear systems (4.2) being solved exactly and all the algebraic computations being performed in exact arithmetic, Theorem 6.1, in particular, implies: Corollary 6.2 (Convergence of the hp adaptive algorithm).Let the Assumptions of Theorem 6.1 be satisfied.Then lim We finish this section by several remarks.

Remark 6.3 (p-robust computable reduction factors C ,red
).We remark that the computable reduction factors C ,red from (6.2) are bounded by the constant C red that is independent of the polynomial degrees p .Moreover, in contrast to the available results, see, e.g., the reduction properties in [43,55,22,19,16,17], the reduction factor C ,red in (6.1) remains fully computable as in [25].
Remark 6.4 (Guaranteed contraction for the "interior node, p+d−1" hp refinement strategy).The "adaptive residual lifting" strategy of Section 4.3.2 that we put forward in Theorem 6.1 relies on a sufficiently fine h-and/or p-refinement of each marked patch ω a , a ∈ V , in order to satisfy condition (5.8).Unfortunately, the proof of Proposition 5.1 does not give an idea of how fine they need to be, neither how to obtain them in an optimal way.We thus do not know how many "one bisection/increase of p be one" iterations are possibly needed in the strategy of Section 4.3.2.In the numerical experiments in Section 10 below, though, only one iteration in each marked patch was always enough.The "interior node, p + d − 1" refinement strategy of Section 4.3.3(at least when there are no data oscillations in that f is a piecewise polynomial and A is piecewise constant), in turn, gives an a priori quantification of how fine the patch h-and p-refinements need to be: the interior node property gives rise to a contraction factor C red as in (6.2), where, however, C lb is now the constant from (5.16) that additionally possibly depends on the polynomial degrees p .Thus, the strategy of Section 4.3.3 is not proven p-robust, and, moreover, theoretically requests p to be bounded by a maximal polynomial degree p max .We also remark that Propositions 5.1 and 5.2 might be seen as two possible proofs of assumption ( 8) from Bürg and Dörfler [15] in the current context.Remark 6.5 (Patch p-refinement).A p-robust contraction is discussed in Canuto et al. [17], for a prefinement strategy under an assumption of a local saturation property on each marked patch.This assumption is similar to (5.8); disregarding data oscillation, the right-hand side of property (iii) in [17,Section 4] is actually equivalent to the left-hand side of (5.8), but the left-hand side of property (iii) is stronger (smaller) than the right-hand side of (5.8), working with H 1 * (ω a ) from (5.4) in place of H 1 0 (ω a ).It is shown in [17,Theorem 6.3] that this local saturation property is a consequence of three elementwise saturation properties.It is then conjectured and numerically illustrated in [17] on a reference triangle, and proven for one of these properties in [18] on a reference square, that these three elementwise properties request a p-refinement of a form p +1,K = p ,K + λp ,K for a parameter λ.The polynomial degree increase has thus theoretically be proportional to p ,K .This stands in contrast to the (usual) increment of the local polynomial degree by a constant factor, as we for instance use in the "interior node, p + d − 1" refinement strategy of Section 4.3.3.The "sufficient p-refinement" of each marked patch requested in theory in our "adaptive residual lifting" strategy of Section 4.3.2is, in turn, related to this.We, however, do not observe a need for such an aggressive p-refinement in our experiments in Section 10 below, where p +1,K = p ,K + 1 is numerically sufficient to satisfy (5.8) and thus for a p-robust contraction factor C ,red .Thus, the theoretical state-of-the-art strong p-refinement of the form p +1,K = p ,K + λp ,K may not be necessary in practice.Remark 6.6 (hp-decision).We confess that the hp-decision on line 9 of Algorithm 1, based on r a,h and r a,p from (4.10), is currently not exploited in the proof of the contraction property (6.1) of Theorem 6.1.On the other hand, this heuristic choice does have an important motivation in that the proof of Proposition 5.1 hinges on the limit (5.15),where clearly the choice of the bigger of the two quantities A indicates a good direction.Currently, our hp-decision does not reflect the cost of h and p variants, which might be included following the ideas in Bürg and Dörfler [15].The hp-decision on line 9 of Algorithm 1 is also purely local and based on residual lifting, whereas some popular hp-decision strategies would rather try to estimate the local regularity of the solution and ensure a steepest descent globally.
Remark 6.7 (Data oscillation).Whenever the datum f is elementwise smooth, the data oscillation terms A,K ) can converge two orders of magnitude faster than the error and estimators, cf.[34, Remark 3.6] or [28].Indeed, recall that p est a is the maximum of the polynomial degrees p ,K over all elements sharing the vertex a.Then, from [28,Lemma 3.6], where the minimum goes over all vertices V ,K of the element K. Thus, f minus the divergence of σ is orthogonal to at least the polynomials of order p ,K in each mesh element K ∈ T .Then, at least when the polynomial degrees are uniform around K, for a regular-enough right-hand side f (one additional division by p ,K is lost in (4.6)where only the division by the Poincaré constant π is used for the purpose of obtaining of a computable constant).As for η osc,a from (5.6), we similarly have, at least when A is piecewise constant, which still converges one order of magnitude faster (in terms of h K /p est a ) than the error and estimators.Consequently, condition (5.7) is likely to be satisfied in practice, see also the discussion in Canuto et al. [17,Remark 4.2].Possibilities for avoiding a condition of the form (5.7) in convergence and optimality proofs are discussed in, e.g., [55,21,19], including a choice of separate marking for data oscillation.However, such options do not necessarily lead to contraction on each step that we achieve here in Theorem 6.1.For general less regular datum f , down to f ∈ H −1 (Ω), the discussion is much more subtle, because in such a case, f − ∇•σ K or g a − Π a (g a ) K are not even defined.For this reason, f ∈ L 2 (Ω) is typically assumed in a posteriori error analysis.For f ∈ H −1 (Ω) with a known structure, we refer to [32] and the references therein.For an in-depth discussion and remedies for a general right-hand side function f ∈ H −1 (Ω), we refer to Kreuzer and Veeser [40] and the references therein.Remark 6.8 (Optimality).We do not study here any question of optimality, i.e., how fast the error A 1 2 ∇(u − u ex +1 ) goes down in terms of the number of degrees of freedom or actual computational cost.We would like to thank the anonymous referee for pointing to us that in particular a convergence rate exponential in the number of degrees of freedom is incompatible with bulk-chasing (Dörfler) criterion of the form (4.7) and a contraction property such as (6.1).Remark 6.9 (A priori hp-refinement strategies).If the location of (corner and/or edge) singularities is known beforehand, than a priori hp-refinement strategies as in Babuška and Guo [4,5] or Schötzau and Schwab [52] will immediately lead to an exponential convergence with respect to the number of degrees of freedom.

The hp-adaptive algorithm with inexact solver
We now recall the building blocks of an inexact hp-adaptive algorithm of Scheme 2, by summing up the necessary ingredients and notation of modules ONE SOLVER STEP and ESTIMATE from [26, Sections 4.1 and 4.2].A new version of the adaptive stopping criterion controlling the adaptive sub-loop in Scheme 2 is defined in Section 7.2.In the present inexact solver setting, we crucially employ slightly modified modules MARK and REFINE, described in Section 7.3.Finally, Section 7.4 generalizes to the inexact solver setting the contents of Section 4.3.4 and presents a first error reduction estimate which highlights the issues to overcome.

Adaptive sub-loop of ONE SOLVER STEP and ESTIMATE
On the initial step = 0, let the exact algebraic solver and hp-refinement criteria from Section 4 be used.Let henceforth ≥ 1 be the current iteration number.We employ the module ONE SOLVER STEP exactly as in [26], with its output being an approximation u := N n=1 (U ) n ψ n ∈ V of the unavailable exact Galerkin approximation u ex given by (4.1).We recall the notation in (4.2) and let the algebraic residual vector R associated with U be given by R := F − A U .
Moreover, we employ its functional representation r ∈ P p (T ), r For details regarding the construction of such r , we refer to Papež et al. [47,Section 5.1].Note that the property (7.2) together with the algebraic system (4.2) yield the functional equivalent of the algebraic relation that we solely use here, so that the forthcoming developments are independent of the choice of the basis ψ n in (4.2).
Definition 7.1 (Discretization flux reconstruction σ ,dis by local minimization).Let u be the current inexact approximation.For each vertex a ∈ V , let the local Raviart-Thomas space V a be prescribed by (4.3) with, however, p est a increased by one (this is related to p-robustness as discussed in [46,Remark 7.10]).Let the local discretization flux reconstruction σ a ,dis ∈ V a be given by the following local minimization problem Then, the discretization flux reconstruction is defined as σ ,dis := a∈V σ a ,dis , with each local contribution σ a ,dis being extended by zero outside of ω a .Note that the Neumann compatibility condition for the problem (7.4) is a direct consequence of (7.3).
Moreover, using the multilevel construction described in [46, Section 4.1] and [26, Section 4.2.1],with, again, the reconstruction polynomial degrees increased by one, we build the algebraic error flux reconstruction σ ,alg (a piecewise in Raviart-Thomas polynomial) in H(div, Ω), such that, thanks to this increase, Then, as established in [26], cf. also [33,46], the following upper bounds on the total and algebraic error hold: , (7.6a) Next, for each vertex a ∈ V , let V a be a suitable finite-dimensional subspace of H 1 * (ω a ) from (5.4).Following [47, Theorem 2], cf. also [26, Theorem 3.5 and Section 4.2.3],we construct the total residual lifting ρ ,tot := a∈V ψ a ρ a ,tot ∈ H 1 0 (Ω), where each vertex contribution ρ a ,tot ∈ V a solves the local primal finite element problem Then, besides the upper bounds (7.6a) and (7.6b), the following lower bound on the total error is at our disposal, supposing ρ ,tot = 0,

Adaptive stopping criterion for the algebraic solver
The above error estimators computed within module ESTIMATE enable us to stop the adaptive sub-loop between the modules ONE SOLVER STEP and ESTIMATE once the current approximation u is such that i.e. when the algebraic error is γ -times smaller than the discretization error.Recalling the Galerkin orthogonality of u ex , A , requirement (7.8) is ensured via (7.6b) and (7.7) by the following criterion under the condition η alg (u , T ) ≤ µ(u ).As already used in [47, Section 6.3] and [26, Section 4.3], (7.9) is equivalent to criterion η alg (u , T ) ≤ γ µ(u ) (7.10) The choice of the parameter γ in (7.10) will follow from the theoretical analysis to be carried out in Section 9 below.

Modules MARK and REFINE
For the sake of brevity, we keep the notation introduced in Section 4.2 unchanged.The module MARK in the inexact solver setting takes as input an inexact approximation u and the corresponding error indicators computed within module ESTIMATE.It proceeds again in two phases.However, in the marking criterion (4.7) we employ in place of the total error estimator η(u , •) only its component corresponding to the discretization error, so the modified bulk chasing criterion to select a subset of marked vertices V θ ⊂ V reads where , for any subset S ⊂ T , and θ ∈ (0, 1] is a fixed threshold parameter.The rest of the module follows straightforwardly the steps described in Section 4.2 leading to the extended set of marked vertices V .In the module REFINE, the new pair (T +1 , p +1 ) is determined on the basis of the two local problems (4.10), with an approximate solution u in place of the exact Galerkin approximation u ex (recall that u ex is not at our disposal here), by one of the strategies of Sections 4.3.2 or 4.3.3.

Discrete lower bound on the incremental error on marked simplices
Proceeding as in the exact solver setting, for each vertex from the extended set of marked vertices V , let the residual lifting r a,hp ∈ V a,hp be defined as the solution of the local problem (4.13), with u ex replaced by the inexact approximation u : with the local space V a,hp given by (4.12).Extending the residual liftings r a,hp by zero outside ω a , let u ex +1 ∈ V +1 be the (unavailable) exact Galerkin approximation on the next level.By [26,Lemma 5.1] used with the extended set V in place of V θ , cf. (4.17), the following lower bound holds true We note that the localization in the spirit of (4.18) is again available here.
Remark 7.2 (Estimate from [26] and issues to overcome).With the results available above, a computable guaranteed bound on the error reduction factor between the errors in the inexact approximations u +1 and u has been given in [26,Theorem 5.4].However, this reduction factor was only less than or equal to one.
In what follows, we design an analysis giving a error reduction factor uniformly smaller than one, and this under a realistic condition on the stopping parameter γ from (7.10).
8 Discrete stability of equilibrated fluxes in an inexact solver setting In this section, we extend the results of Section 5, namely the local stabilities (5.8) and (5.16), to the inexact solver setting.We use the notation defined therein.We then also prove a global stability of the form (6.7).

Local stability
Recall that η osc,a is defined by (5.6), with g a given by (4.5).The counterpart of Proposition 5.1 is: Proposition 8.1 (Discrete local p-robust stability of the flux equilibration, inexact solver).Let u ∈ V satisfy the hat function orthogonality and let the local equilibrated flux σ a ,dis be constructed by (7.4).Let data oscillation do not dominate in that Let the algebraic error flux reconstruction σ ,alg ∈ H(div, Ω) satisfy (7.5).Then, for a sufficiently fine patch space V a,hp from (4.12) and r a,hp given by (7.12), there holds 2) for all a ∈ V , where in particular C lb can be chosen as the maximum of (5.9) and 1.
Proof.The proof follows that of Proposition 5.1, upon replacing u ex by u and f by f − r , and while working with the algebraic residual representation r as in [46,proof of Theorem 7.7].In particular, step 1 yielding (5.10) stays intact.On step 2, it is crucial that, thanks to the increase of p est a by one in (7.4), Π a (ψ a r ) = ψ a r .Then, on step 3, relying on (7.5), we can proceed as for [46, eq. (B.4)], so that, on step 4, we obtain from where the assertion follows by virtue of assumption (8.1).
Similarly, the counterpart of Proposition 5.2 is: Proposition 8.2 (Discrete local non p-robust stability of the flux equilibration, inexact solver).Let f ∈ P p f (T 0 ) fulfill p f ≤ p 0 − 1 and let the diffusion tensor A be piecewise constant with respect to T 0 , A ∈ P 0 (T 0 ) d×d .Let the approximate solution u ∈ V satisfy the hat function orthogonality Let the local discretization flux reconstruction σ a ,dis be constructed by (7.4) and let the algebraic error flux reconstruction σ ,alg , a piecewise Raviart-Thomas polynomial in H(div, Ω), satisfy (7.5).Finally, let the patch space V a,hp from (4.12) involve the interior nodes (4.15) or the polynomial degree increase (4.16) and let r a,hp be given by (7.12).Then there exists a constant C lb ≥ 1 only depending on the space dimension d, the mesh shape-regularity κ T , the ratio of the largest and the smallest eigenvalue of the diffusion coefficient A, and the maximal polynomial degree p max such that Proof.The proof again follows the proof of Proposition 5.2 and of [46, proof of Theorem 7.7].In particular, in place of (5.38), we estimate and to treat the term r K , we use property (7.5) of the algebraic error flux reconstruction σ ,alg and an inverse inequality yielding Remark 8.4 (Boundedness of ξ from above).Observe that unless η(u , T ) = 0, in which case u = u and the adaptive loop terminates, (7.11) together with (9.7) below implies η dis+osc (u , M θ ) = 0 and thus ξ < ∞; actually the first two inequalities in (9.9) below give θ 2 ξ 2 < 2.
Proof of Lemma 8.3.We proceed similarly as in step 2 of the proof of Theorem 6.1.Remark from (7.4), (7.5), and (4.5) that on each element K ∈ T , Thus, employing the notation from (7.6a) and that below (7.11), we obtain as in (6.6) which yields the second inequality in the assertion (8.7) and concludes the proof.

Proof of guaranteed contraction with an inexact solver
The second main result of this paper, completing Theorem 6.1 in the inexact solver setting, is summarized in the following.Recall that θ ∈ (0, 1] is the Dörfler marking parameter from (7.11) and the notation ξ and ξ from (8.7).
Let u be the weak solution of (1.1) and let the current inexact approximation u , ≥ 0, be obtained by the adaptive sub-loop of Section 7.1 employing the stopping criterion (7.10) with the parameter γ satisfying Owing to the exact initial solve, set γ = 0 and σ ,alg = 0 for = 0. Let the pair (T +1 , p +1 ) be obtained by the module REFINE of Section 7.3 with the "adaptive residual lifting" strategy of Section 4.3.2 and let η * M be the computable discrete lower bound defined by (7.13).Moreover, let u +1 ∈ V +1 be the inexact finite element approximation on iteration + 1 of the adaptive loop prescribed by Scheme 2, satisfying the stopping criterion (7.10) and conditions (9.1) with the iteration counter + 1 in place of .Let finally the data oscillation do not dominate in that (8.1) holds.Then, two options arise.Either the total a posteriori error estimate η(u , T ) from (7.6a) is zero, in which case u = u, and the adaptive loop terminates.Or so that C ,red is a fully computable bound on the error reduction factor and the generic constant C # red only depends on the marking parameter θ, the space dimension d, the mesh shape-regularity κ T , and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient A; here, C lb is the constant from (8.2) that can theoretically be chosen as the maximum of (5.9) and 1.
Proof.Similarly the proof of [26,Theorem 5.4] but using the Galerkin orthogonality, we have We proceed in five steps.
Step 1 (bound on A 1 2 ∇(u−u ex +1 ) ).Let ≥ 1 and consider only the non trivial case when η(u , T ) = 0. Note first that the Galerkin orthogonality of u ex +1 gives Employing the discrete lower bound (7.13), adding and subtracting A − 1 2 σ ,alg 2 , using the adaptive stopping criterion (7.10), and relying on the total error bounds (7.7) and (7.6a), we obtain As we suppose an exact solve on the initial mesh, the claim (9.6) for = 0 immediately holds with γ = 0 and σ ,alg = 0.
Thus, from (9.9), we infer where we have also employed the bound (8.7) from Lemma 8.3.
In extension of Corollary 6.2, we observe that Theorem 9.1 implies: Corollary 9.2 (Convergence of the adaptive algorithm with inexact solvers).Let the assumptions of Theorem 9.1 be satisfied.Then lim Remark 9.3 (Consistency with Theorem 6.1).Theorem 9.1 is a consistent extension of Theorem 6.1, in that for an exact solver, where γ = γ +1 = 0 and σ ,alg = 0, they give the same contraction factor C ,red when the Dörfler criterion (4.7) is employed with equality, i.e., η(u ex , M θ ) = θ η(u ex , T ). Remark This would give a valid (theoretical, since ξ from (8.7) is generally unknown) choice for γ , but we rather stress that (9.10) shows that picking, e.g., γ as a fixed fraction of the min in (9.1) is uniformly bounded from below in all the iteration counter , the mesh sizes, and the polynomial degrees.
Remark 9.6 (Contraction factors C red , C * red , and C # red ).We remark that the contraction factors C red from (6.2), C * red from (9.9), and C # red from (9.3) have all the same structure, upon the coefficients 1, 4, and 8; this was actually our intention beyond the design of condition (9.1).A crucial consequence is that we do not loose (p-robust) contraction with inexact solvers.

Numerical illustration
We illustrate here the proposed exact and inexact hp-adaptive algorithms on a two-dimensional test case with a singular weak solution.The setting is taken from [25,26], where ample numerical illustrations were supplied for the hp-adaptive strategies that we call henceforth "practical".Our goal is to show that their theory-prone modifications, namely the "adaptive residual lifting" hp strategy of Section 4.3.2 and the "interior node, p + d − 1" hp strategy of Section 4.3.3,lead to structurally the same behavior, though with typically worse ratio price/outcome.
Consider the classical re-entrant corner problem, cf.[42,25,26], posed on the L-shape domain Ω We start on a coarse criss-cross grid T 0 with max K∈T0 h K = 0.25 and all the polynomial degrees set to 1.We first consider an exact algebraic solver and assess the overall quality of the "adaptive residual lifting" hp strategy of Section 4.3.2, as well as of the "interior node, p + d − 1" hp strategy of Section 4.3.3.Figure 4, left, traces the decrease of the relative energy error ∇(u − u ex ) / ∇u in function of the number of degrees of freedom.Asymptotic exponential convergence rates are observed.In comparison with the "practical" strategy of [25, Section 6.2], however, worse constants and a worse ratio price/outcome are observed.In the "practical" strategy, there in particular holds ∇ with C 1 = 4.73 and C 2 = 0.69, whereas in the "adaptive residual lifting" and "interior node, p + d − 1" settings, the (second) constants respectively worsen to C 1 = 0.51, C 2 = 0.30 and C 1 = 0.64, C 2 = 0.31.This is explained by the present additional assumptions on the h-and p-refinements, namely the extension of the marked region, the interior node property, and the steeper polynomial degree increase.The last adaptively refined mesh T 25 29 and the corresponding distribution of the polynomial degrees for the "adaptive residual lifting" strategy are then presented in Figure 4, right.Unfortunately, an excessive p-refinement close to the re-entrant corner is produced.On the positive side, with C lb = 10, criterion (4.14) for the "adaptive residual lifting" strategy of Section 4.3.2 is always satisfied, so that no iteration until (4.14) is satisfied is necessary in Algorithm 1.
As in [25], we next investigate the quality of C ,red , the guaranteed bound on the error reduction factor from Theorem 6.1.We do so on all hp-refinement iterations in terms of the effectivity index defined as . (10.1) We also verify the sharpness of the underlying discrete lower bound η M given by (4.17) in terms of the effectivity index defined as The results for both strategies of Sections 4.3.2 and 4.3.3 are reported in Figure 5, where we also plot the corresponding "practical" effectivity indices taken from [25,Figure 15].Very similar and close-to-optimal results can be observed at all iterations.Moreover, in contrast to above, the two theory-prone strategies lead to a slightly better overall performance.We next investigate criterion (4.14) for the "adaptive residual lifting" hp strategy of Section 4.3.2 in more details.Similarly as in (8.6) in the inexact solver setting (and recalling that f = 0 and A = I here, so that no data oscillations arise from f and A), let us define Figure 5: Effectivity index for the error reduction factor estimate C ,red of Theorem 6.1 given by (10.1) (left); effectivity index for the discrete lower bound η M of (4.17) given by (10.2) (right).The "adaptive residual lifting" hp strategy of Section 4.3.2, the "interior node, p+d−1" hp strategy of Section 4.3.3,and the "practical" hp strategy from [25,Figure 15].Exact algebraic solver.both C lb and C lb in Figure 6, left (recall from Figure 4 that the maximal polynomial degree produced by the strategy on the last mesh T 25 equals 8).We see that C lb ≤ 2.2 < 10, which explains why no iteration on h-or p-refinements appears above for our choice C lb = 10.Moreover, in Figure 6, right, we test the (artificial) case where we always increase the polynomial degrees to p + 2 in place of p + 1, which leads to the maximal polynomial degree 15 on the last mesh.We see that the values of C lb and C lb essentially do not change with respect to the left part of Figure 6, with C lb ≤ 2.25 < 10.Thus, though the p-robust guaranteed contraction of Theorem 6.1 a priori needs more h-and/or p-refinements in each marked patch to satisfy criterion (4.14), one minimal h-and/or p-refinement is sufficient in the present case with C lb = 10.
To support this numerical observation further, we plot in Figure 7, left, C lb and C lb from (10.3) also for the "sharp Gaussian" smooth solution from [25, Section 6.1], which shows a similar behavior.
We finally quickly illustrate the inexact solver case.Similar results as above for the exact solver are observed in terms of exponential relative energy error decay, hp-meshes, the effectivity indices of C ,red and of η M , and the values of C lb as well as of C lb .We thus only focus on the stopping parameter γ that we choose following (9.1) in order to satisfy the theoretical assumption needed in the contraction of Theorem 9.1.In Figure 7, right, we plot the evolution of γ during the hp adaptation iterations for the "interior node, p + d − 1" strategy of Section 4.3.3.The theoretical requirement (9.1) is more stringent than the "practical" requirement of [26, equation (5.12)], also plotted in Figure 7, right (data taken from [26, Figure 20 (right)]).Regardless, comparable results in both cases are observed.In particular, γ satisfying (9.1) still takes the very reasonable (in practice) values around 0.05, which leads to requesting at most 5% algebraic error in the total error, recalling (7.10).given by C lb and C lb from (10.3), the "adaptive residual lifting" hp strategy of Section 4.3.2 for the "sharp Gaussian" smooth solution from [25,Section 6.1].Exact algebraic solver (left).Values of the parameter γ , the "interior node, p+d−1" hp strategy of Section 7 compared to the parameter γ as in [26,Theorem 5.4].Inexact algebraic solver (right).

Conclusions and outlook
The contribution of this work is a theoretical study of the adaptive hp-refinement strategies for conforming finite elements proposed in [25] for exact algebraic solvers and in [26] for inexact algebraic solvers.Under some modifications of these two methods (extension of the marked patches by one layer, increase of the equilibration polynomial degree by one for inexact solvers) and under the criterion (4.14) that the algorithm possibly practically checks on each marked patch, we have in particular shown a strict error reduction on each step of the adaptive loop, and consequently convergence of the resulting adaptive algorithms.Abstracting from data oscillations, we have theoretically proven that criterion (4.14) is satisfied 1) uniformly in the iteration index but not uniformly in the polynomial degree p when interior nodes have been introduced and the polynomial degree has been increased to p+d−1 (by the bubble function technique); 2) uniformly in both and p when each marked patch is sufficiently h-and p-refined (via polynomial extension operators).In the numerical experiments we have carried out, though, criterion (4.14) is p-robustly satisfied (and thus p-robust error reduction appears) already for one newest-vertex bisection step and/or increase of the polynomial degree by one, as in [25,26].In the case of an inexact algebraic solver, we have identified that the stopping coefficients γ have to theoretically verify condition (9.1).This turns out not to request an excessively small tolerance and leads in our numerical experiments to values of γ around 0.05, meaning that the algebraic iterations can be stopped when the algebraic error is roughly 20-times smaller than the total one.Some important limitations of the proposed algorithms and analysis remain.The hp-decision on line 9 of Algorithm 1 is not exploited in the theoretical proof of the contraction, we have not been able to identify the best possible way to satisfy criterion (4.14), and neither were we able to prove any optimality result.Moreover, even though we numerically observe exponential convergence rates with respect to the number of degrees of freedom, our hp decision criterion is not performing to our entire satisfaction, leading to an increased p-refinement towards the singularities.

1 2 ∇
(u − u ex ) .The module ESTIMATE takes as input the finite element solution u ex and outputs a collection of local error indicators {η K } K∈T .

Figure 1 :
Figure 1: An example of local error estimators ηK (u ex ) from (4.6) (left) and illustration of the corresponding set of marked vertices V θ = {a1} and its extension V = {a1, a2, . . ., a8}, θ = 0.5 (right).The region highlighted in the red color corresponds to the original marked subdomain ω , while its union with the yellow region amounts to the extended marked subdomain ω .

Proposition 5 . 1 (
Discrete local p-robust stability of the flux equilibration).Let u ex ∈ V satisfy the hat function orthogonality (4.15)  and a typical polynomial degree increase by d − 1, see (4.16), are sufficient for guaranteed contraction.Upon relying on Proposition 5.2 in place of Proposition 5.1, it and the weak solution, in polar coordinates, u(r, ϕ) = r

Figure 4 :
Figure 4: Relative energy error ∇(u − u ex ) / ∇u as a function of DoF 1 3 , obtained with the present "adaptive residual lifting" hp strategy of Section 4.3.2, the "interior node, p + d − 1" hp strategy of Section 4.3.3,and the "practical" hp strategy from [25].Also plotted are the pure h-adaptive mesh refinement and hp-adaptation exploiting the a priori knowledge of the weak solution (left).The last adaptively refined mesh T25 and the corresponding distribution of the polynomial degrees (whole domain and [−1e − 2, 1e − 2] 2 crop), the "adaptive residual lifting" hp strategy of Section 4.3.2(right).Exact algebraic solver.

1 2
∇r a,hp ω a ∀a ∈ V , C lb := max a∈V C a lb , and C lb := min a∈V one newest-vertex bisection step or polynomial increase to p + 1 in each marked patch ω a .In other words, we just investigate the ratio of the left-hand side to the right-hand side in (4.14) or (5.8) after the first h-or p-refinement iteration of the "adaptive residual lifting" hp strategy of Section 4.3.2.We plot (practical) hp-adaptivity (adaptive resid.lifting) hp-adaptivity (interior node, p+d-1) (practical) hp-adaptivity (adaptive resid.lifting) hp-adaptivity (interior node, p+d-1)

Figure 6 :
Figure 6: Maximal and minimal values over all marked vertices of the local stability constants C a lb given by C lb and C lb from (10.3).The "adaptive residual lifting" hp strategy of Section 4.3.2corresponding to Figure 4 with maximal polynomial degree 8 (left), maximal polynomial degree 15 (right).Exact algebraic solver.

Figure 7 :
Figure 7: Maximal and minimal values over all marked vertices of the local stability constants C a lb given by C lb and C lb from (10.3), the "adaptive residual lifting" hp strategy of Section 4.3.2 for the "sharp Gaussian" smooth solution from [25, Section 6.1].Exact algebraic solver (left).Values of the parameter γ , the "interior node, p+d−1" hp strategy of Section 7 compared to the parameter γ as in [26, Theorem 5.4].Inexact algebraic solver (right).
Lemma 8.3 (Discrete global stability of the flux equilibration, inexact solver).Let ≥ 1, let η dis+osc (u , M θ ) be given by (7.11), let η * M be the computable discrete lower bound defined by (7.13), and consider either the setting of Proposition 8.1, or that of Proposition 8.2.Then, there holds . Recall now that C lb from (8.6) are bounded by C lb by Propositions 8.1 or 8.2, which yields the first inequality in the assertion (8.7).Then, relying on (8.6) in (8.8) and using the inequality (a + b) 2 ≤ 2a 2 + 2b 2 together with the localization of the lower bound η * M in the spirit of (4.18) leads to 9.4 (The "interior node, p + d − 1" hp refinement strategy).As in Remark 6.4, the result of Theorem 9.1 holds true for the "interior node, p + d − 1" strategy of Section 4.3.3when the polynomial degrees are bounded by p max , with the non p-robust constant C lb from (8.3) in place of the p-robust constant C lb from (8.2).Remark 9.5 (Stopping parameter γ ).Our analysis shows that choosing the stopping parameter γ following (9.1)(practically, since both θ and ξ from (8.7) are available) is sufficient to obtain the uniform contraction (9.2).Criterion (9.1) does not seem to request excessively small values of γ , leading in particular to the reasonable choice of γ ≈ 0.05 in the numerical experiments in Section 10 below, in contrast to (to our knowledge) all results available in the literature so far.Moreover, from (8.7), there holds