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

Patrik Daniel, Martin Vohralík

<figure><img src="image_1.png" /></figure>

## To cite this version:

Patrik Daniel, Martin Vohralík. Guaranteed contraction of adaptive inexact hp-refinement strategies with realistic stopping criteria. ESAIM: Mathematical Modelling and Numerical Analysis, 2023, 57 (1), pp.329 - 366. 10.1051/m2an/2022082. hal-02486433v3

HAL Id: hal-02486433

https://inria.hal.science/hal-02486433v3

Submitted on 23 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Guaranteed contraction of adaptive inexact hp-refinement strategies with realistic stopping criteria*

Patrik Daniel†‡

Martin Vohralík†‡

July 23, 2022

## Abstract

The purpose of this contribution is to theoretically analyze the adaptive refinement strategies for conforming hp-finite element approximations of elliptic problems proposed for exact algebraic solvers in [Daniel, Ern, Smears, Vohralík, Comput. Math. Appl. 76 (2018), 967–983] and for inexact algebraic solvers in [Daniel, Ern, Vohralík, Comput. Methods Appl. Mech. Engrg. 359 (2020), 112607]. Both of these strategies are driven by guaranteed equilibrated flux energy error estimators. The employed hp-refinement criterion stems from solving two separate local residual problems posed only on the patches of elements around marked vertices selected by a bulk-chasing criterion. In the above references, we have derived a fully computable guaranteed bound on the ratio of the error on two successive steps of the hp-adaptive loop. Here, our focus is to prove that this ratio is uniformly smaller than one, and thus the convergence of the adaptive and adaptive inexact hp-refinement strategies. To be able to achieve this goal, we have to introduce some additional assumptions on the h- and p-refinements, namely an extension of the marked region, as well as a sufficient h- or p-refinement of each marked patch. We investigate two such strategies, where one ensures a polynomial-degree-robust guaranteed contraction. In the inexact case, a sufficiently precise stopping criterion for the algebraic solver is requested, but this criterion remains fully computable and also realistic in the sense that in our numerical experiments, it does not request the algebraic error to be excessively small in comparison with the total error.

**Key words:** elliptic problem, finite element method, a posteriori error estimate, equilibrated flux, hp-adaptivity, inexact solver, error reduction, convergence analysis

## 1 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.

<footnote>*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).</footnote>

<footnote>†Inria, 2 rue Simone Iff, 75589 Paris, France</footnote>

<footnote>‡CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée, France</footnote>
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 $\rho_{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 $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, be a polygonal/polyhedral domain (open, bounded, and connected set) with a Lipschitz boundary $\partial\Omega$, and let $H_0^1(\Omega)$ denote the Sobolev space of all functions in $L^2(\Omega)$ which have all their first-order weak derivatives in $L^2(\Omega)$ and a vanishing trace on $\partial\Omega$. Let $f \in L^2(\Omega)$ be a given source term and let $A \in [L^\infty(\Omega)]^{d \times 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 \in H_0^1(\Omega)$ such that

$$ (A\nabla u, \nabla v) = (f, v) \quad \forall v \in H_0^1(\Omega), \qquad (1.1) $$

where $(\cdot, \cdot)$ stands for the $L^2(\Omega)$ or $[L^2(\Omega)]^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.

<figure><img src="image_3.png" /><figcaption>**Scheme 1:** Paradigm of an hp-adaptive algorithm employing an exact algebraic solver.</figcaption></figure>

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 $\ell$ of the adaptive loop from Scheme 1, it is possible to compute explicitly a real number $0 \le C_{\ell,\text{red}} \le 1$ such that

$$ \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell+1}^{\text{ex}}) \| \le C_{\ell,\text{red}} \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}^{\text{ex}}) \|, \qquad (1.2) $$

where $u_{\ell}^{\text{ex}}$ and $u_{\ell+1}^{\text{ex}}$ denote the exact finite element solutions from the respective iterations $\ell$ and $\ell+1$ of the adaptive loop. Even though we provided numerical evidence in [25, Section 6], we did not prove that the reduction factor $C_{\ell,\text{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 $p$-refinements 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 first 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.

<figure><img src="image_4.png" /><figcaption>**Scheme 2:** Paradigm of an adaptive loop employing an inexact algebraic solver.</figcaption></figure>

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 $\gamma_\ell$, 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_\ell$ and $u_{\ell+1}$ produced by the adaptive loop of Scheme 2

$$ \| \mathbf{A}^{\frac{1}{2}} \nabla(u - u_{\ell+1}) \| \le C_{\ell, \text{red}} \| \mathbf{A}^{\frac{1}{2}} \nabla(u - u_{\ell}) \|, \quad (1.3) $$

with a fully computable factor $0 \le C_{\ell, \text{red}} \le 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 $\gamma_\ell$ 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.
## 2 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_{\ell}^{\text{ex}}$ be the (known) finite element approximation on mesh $\mathcal{T}_{\ell}$ and space $V_{\ell}$, and let $u_{\ell+1}^{\text{ex}}$ be the (unknown) finite element approximation on the forthcoming mesh $\mathcal{T}_{\ell+1}$ and space $V_{\ell+1}$ (containing $V_{\ell}$) that has been designed by our hp-refinement strategy on the current step $\ell$. The Galerkin orthogonality of $u_{\ell+1}^{\text{ex}}$ yields

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\|^2 = \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\|^2 - \|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}})\|^2. \quad (2.1) $$

In our approach, we use an equilibrated flux reconstruction $\sigma_{\ell}^{\text{a}}$ on each patch subdomain $\omega_{\ell}^{\text{a}}$ associated with a vertex $\mathbf{a} \in \mathcal{V}_{\ell}$; $\sigma_{\ell}^{\text{a}}$ is a Raviart-Thomas piecewise polynomial with no flow through the boundary of $\omega_{\ell}^{\text{a}}$ on interior vertices. We then set $\sigma_{\ell} = \sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \sigma_{\ell}^{\text{a}}$, which gives the guaranteed error upper bound

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\| \le \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}) := \|\mathbf{A}^{\frac{1}{2}}\nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}}\sigma_{\ell}\|. \quad (2.2) $$

We use $\|\mathbf{A}^{\frac{1}{2}}\nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}}\sigma_{\ell}\|_{K}$ for marking of mesh elements $K$ to h- or p- refinement; more precisely, the set of marked elements $\mathcal{M}_{\ell}^{\theta}$ has to contain a $\theta$-fraction of the overall estimated error via the bulk-chasing (Dörfler) criterion

$$ \eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta}) \ge \theta \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}). \quad (2.3) $$

Moreover, on each patch $\omega_{\ell}^{\text{a}}$ around a vertex $\mathbf{a} \in \mathcal{V}_{\ell}^{\text{a}}$ marked for refinement, we construct a residual lifting $r^{\text{a},hp}$, which is a conforming piecewise polynomial taking zero value on the boundary of $\omega_{\ell}^{\text{a}}$. Then we have the guaranteed discrete lower bound

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}})\| \ge \eta_{\mathcal{M}_{\ell}^{\text{a}}} = \frac{\sum_{\mathbf{a} \in \mathcal{V}_{\ell}^{\text{a}}} \|\mathbf{A}^{\frac{1}{2}}\nabla r^{\text{a},hp}\|_{\omega_{\ell}^{\text{a}}}^{2}}{\|\mathbf{A}^{\frac{1}{2}}\nabla\left(\sum_{\mathbf{a} \in \mathcal{V}_{\ell}^{\text{a}}} r^{\text{a},hp}\right)\|_{\omega_{\ell}^{\text{a}}}}. \quad (2.4) $$

Now, we can already bound the error $\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\|$, with respect to $\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\|$, prior to the computation of $u_{\ell+1}^{\text{ex}}$. 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

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}})\| \stackrel{(2.4)}{\ge} \eta_{\mathcal{M}_{\ell}^{\text{a}}} = \frac{\eta_{\mathcal{M}_{\ell}^{\text{a}}}}{\eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell})} \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}) \stackrel{(2.3)}{\ge} \theta \frac{\eta_{\mathcal{M}_{\ell}^{\text{a}}}}{\eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})} \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}) \quad (2.5a) $$

$$ \stackrel{(2.2)}{\ge} \theta \frac{\eta_{\mathcal{M}_{\ell}^{\text{a}}}}{\eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})} \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\|. \quad (2.5b) $$

Thus

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\| \le C_{\ell,\text{red}} \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\| \quad (2.6) $$

for the computable constant

$$ 0 \le C_{\ell,\text{red}} := \sqrt{1 - \theta^2 \frac{\eta_{\mathcal{M}_{\ell}^{\text{a}}}^2}{\eta^2(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})}} \le 1. $$

By (2.6), we can monitor whether a contraction in the finite element approximation is happening, but an important question is whether $C_{\ell,\text{red}}$ is uniformly bounded by a constant $C_{\text{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^{\text{a},hp}$ such that the local stability

$$ \|\psi_{\ell}^{\text{a}}\mathbf{A}^{\frac{1}{2}}\nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}}\sigma_{\ell}^{\text{a}}\|_{\omega_{\ell}^{\text{a}}} \le C_{\text{lb}}\|\mathbf{A}^{\frac{1}{2}}\nabla r^{\text{a},hp}\|_{\omega_{\ell}^{\text{a}}} \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\text{a}} \quad (2.7) $$

holds for a generic constant $C_{\text{lb}}$. Then, overlaps between patches lead to

$$ \eta_{\mathcal{M}_{\ell}^{\text{a}}} \ge \frac{\{\sum_{\mathbf{a} \in \mathcal{V}_{\ell}^{\text{a}}} \|\mathbf{A}^{\frac{1}{2}}\nabla r^{\text{a},hp}\|_{\omega_{\ell}^{\text{a}}}^{2}\}^{1/2}}{\sqrt{d+1}}, \quad (2.8) $$
and patches overlaps again, together with the local stabilities (2.7), imply the global stability on the refined elements collected in $\mathcal{M}_{\ell}^{\theta}$ in the form

$$ \eta^2(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta}) \le (d+1) \sum_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}^{2} \stackrel{(2.7)}{\le} (d+1) \sum_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} C_{\text{lb}}^{2} \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}, \text{hp}}\|_{\omega_{\ell}^{\mathbf{a}}}^{2} \stackrel{(2.8)}{\le} (d+1)^2 C_{\text{lb}}^{2} \eta_{\mathcal{M}_{\ell}^{\sharp}}^{2}. \quad (2.9) $$

From (2.9), we conclude that indeed

$$ 0 \le C_{\ell, \text{red}} = \sqrt{1 - \theta^2 \frac{\eta_{\mathcal{M}_{\ell}^{\sharp}}^2}{\eta^2(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})}} \le C_{\text{red}} := \sqrt{1 - \frac{\theta^2}{(d+1)^2 C_{\text{lb}}^2}} < 1. $$

## 3 Framework and notation

This section details the notation employed throughout the manuscript. Within the adaptive loops of Schemes 1 and 2, a sequence $\{(\mathcal{T}_\ell, \mathbf{p}_\ell)\}_{\ell \ge 0}$, with $\ell \ge 0$ the iteration counter, is generated. Each pair $(\mathcal{T}_\ell, \mathbf{p}_\ell)$ consists of a matching simplicial mesh $\mathcal{T}_\ell$ of the computational domain $\Omega$, i.e., a finite collection of (closed) simplices $K \in \mathcal{T}_\ell$ covering $\bar{\Omega}$ and such that the intersection of two different simplices is either empty or their $d'$-dimensional common face, $0 \le d' \le d-1$, and of a polynomial-degree distribution vector $\mathbf{p}_\ell := \{p_{\ell,K}\}_{K \in \mathcal{T}_\ell}$ which assigns a degree $p_{\ell,K} \in \mathbb{N}_{\ge 1}$ to each simplex $K \in \mathcal{T}_\ell$. Moreover, each pair $(\mathcal{T}_\ell, \mathbf{p}_\ell)$ prescribes a discrete finite-dimensional space $V_\ell$, defined as

$$ V_{\ell} := \mathbb{P}_{\mathbf{p}_{\ell}}(\mathcal{T}_{\ell}) \cap H_{0}^{1}(\Omega), \quad \forall \ell \ge 0, $$

where $\mathbb{P}_{\mathbf{p}_\ell}(\mathcal{T}_\ell)$ denotes the space of piecewise polynomials of total degree at most $p_{\ell,K}$ on each simplex $K \in \mathcal{T}_\ell$. Let us denote by $N_\ell$ the dimension of the $\ell$-th level space $V_\ell$. Note that we enforce the $H_0^1(\Omega)$-conformity of the spaces $V_\ell$ for all $\ell \ge 0$. In addition, we make the following nestedness assumption:

$$ V_{\ell} \subset V_{\ell+1}, \quad \forall \ell \ge 0. \quad (3.1) $$

The initial pair $(\mathcal{T}_0, \mathbf{p}_0)$ is assumed to be given. Then, the purpose of each step $\ell \ge 0$ of the adaptive loops of Schemes 1 and 2 is to determine the next pair $(\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+1})$. The nestedness property (3.1) gives us two crucial restrictions on the meshes and polynomial-degree distributions defining the spaces $V_\ell$: (i) the sequence of meshes $\{\mathcal{T}_\ell\}_{\ell \ge 0}$ needs to be hierarchically nested, i.e., for all $\ell \ge 1$ the mesh $\mathcal{T}_\ell$ is a refinement of $\mathcal{T}_{\ell-1}$ such that for all $K \in \mathcal{T}_\ell$, there is a unique simplex $\tilde{K} \in \mathcal{T}_{\ell-1}$, called the parent of $K$, satisfying $K \subseteq \tilde{K}$; (ii) the local polynomial degree is locally increasing, i.e., for all $\ell \ge 1$ and all $K \in \mathcal{T}_\ell$, $p_{\ell,K} \ge p_{\ell-1, \tilde{K}}$, where $\tilde{K} \in \mathcal{T}_{\ell-1}$ is the parent of $K$. Moreover, we assume the following standard shape-regularity property: there exists a constant $\kappa_T > 0$ such that $\max_{K \in \mathcal{T}_\ell} h_K / \rho_K \le \kappa_T$ for all $\ell \ge 0$, where $h_K$ is the diameter of $K$ and $\rho_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 $\mathcal{V}_\ell$ (respectively $\mathcal{F}_\ell$) the set of vertices (respectively $(d-1)$-dimensional faces) of the mesh $\mathcal{T}_\ell$. These are decomposed into interior vertices $\mathcal{V}_\ell^{\text{int}}$ ($(d-1)$-dimensional faces $\mathcal{F}_\ell^{\text{int}}$) and boundary vertices $\mathcal{V}_\ell^{\text{ext}}$ ($(d-1)$-dimensional faces $\mathcal{F}_\ell^{\text{ext}}$). For each vertex $\mathbf{a} \in \mathcal{V}_\ell$, $\ell \ge 0$, the so-called hat function $\psi_\ell^\mathbf{a}$ is the continuous, piecewise affine function that takes the value 1 at the vertex $\mathbf{a}$ and the value 0 at all the other vertices of $\mathcal{V}_\ell$; the function $\psi_\ell^\mathbf{a}$ is contained in the space $V_\ell$ for all vertices $\mathbf{a} \in \mathcal{V}_\ell^{\text{int}}$ and all polynomial degrees $\mathbf{p}_\ell$. Furthermore, we consider the simplex patch $\mathcal{T}_\ell^\mathbf{a} \subset \mathcal{T}_\ell$ which is the collection of the simplices sharing the vertex $\mathbf{a} \in \mathcal{V}_\ell$, with $\omega_\ell^\mathbf{a}$ the corresponding open subdomain of $\Omega$, coinciding with the support of $\psi_\ell^\mathbf{a}$. Let $\mathcal{F}_\ell^\mathbf{a} \subset \mathcal{F}_\ell$ denote the set of all the $(d-1)$-dimensional faces in the patch $\mathcal{T}_\ell^\mathbf{a}$. This set can be further decomposed into $\mathcal{F}_\ell^{\mathbf{a},\text{int}}$, the subset of faces from $\mathcal{F}_\ell^\mathbf{a}$ that share the vertex $\mathbf{a}$ and are shared by two distinct simplices in $\mathcal{T}_\ell^\mathbf{a}$, and $\mathcal{F}_\ell^{\mathbf{a},\text{ext}}$, the faces from $\mathcal{F}_\ell^\mathbf{a}$ lying in $\partial\omega_\ell^\mathbf{a}$, so that $\mathcal{F}_\ell^\mathbf{a} = \mathcal{F}_\ell^{\mathbf{a},\text{int}} \cup \mathcal{F}_\ell^{\mathbf{a},\text{ext}}$. For any face $F \in \mathcal{F}_\ell$, $\mathbf{n}_F$ stands for a unit vector normal to $F$ with an arbitrary but fixed orientation. The operator $[\cdot]$ yields the jump, in the direction of $\mathbf{n}_F$, of the traces of the argument from the two simplices that share $F \in \mathcal{T}_\ell^{\text{int}}$, and the actual trace for $F \in \mathcal{F}_\ell^{\text{ext}}$. Finally, for each simplex $K \in \mathcal{T}_\ell$, $\mathcal{V}_{\ell,K}$ denotes the set of vertices of $K$ and $\mathcal{F}_{\ell,K}$ denotes the set of $(d-1)$-dimensional faces of element $K \in \mathcal{T}_\ell$.
# 4 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.

## 4.1 The modules SOLVE and ESTIMATE

Let $\ell \ge 0$ denote the current iteration counter. The module SOLVE takes as input the current finite element space $V_\ell \subset H_0^1(\Omega)$ and outputs the Galerkin approximation $u_\ell^{\text{ex}} \in V_\ell$ of the weak solution $u$ of (1.1) defined as the unique solution of

$$ (\mathbf{A}\nabla u_{\ell}^{\text{ex}}, \nabla v_{\ell}) = (f, v_{\ell}) \quad \forall v_{\ell} \in V_{\ell}. \qquad (4.1) $$

Note that obtaining $u_\ell^{\text{ex}}$ is equivalent to computing the exact solution of the system of linear algebraic equations

$$ \mathbf{A}_{\ell} \mathbf{U}_{\ell}^{\text{ex}} = \mathbf{F}_{\ell}, \qquad (4.2) $$

where $\{\psi_\ell^n\}_{1 \le n \le N_\ell}$ is a basis of the $\ell$-th level space $V_\ell$, so that $u_\ell^{\text{ex}} := \sum_{n=1}^{N_\ell} (U_\ell^{\text{ex}})_n \psi_\ell^n$, $(\mathbf{A}_\ell)_{k,m} = (\mathbf{A}\nabla\psi_\ell^m, \nabla\psi_\ell^k)$, and $(\mathbf{F}_\ell)_k = (f, \psi_\ell^k)$, $1 \le k, m \le N_\ell$.

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 $\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_\ell^{\text{ex}})\|$. The module ESTIMATE takes as input the finite element solution $u_\ell^{\text{ex}}$ and outputs a collection of local error indicators $\{\eta_K\}_{K \in \mathcal{T}_\ell}$.

The equilibrated flux is constructed locally on the simplex patches $\mathcal{T}_\ell^a$ attached to each vertex $\mathbf{a} \in V_\ell$. For this construction, we consider as in [28, 25] the local polynomial degree $p_\ell^{\text{est}} := \max_{K \in \mathcal{T}_\ell^a} p_{\ell,K}$ (any other choice so that $p_\ell^{\text{est}} \ge \max_{K \in \mathcal{T}_\ell^a} p_{\ell,K}$ can also be employed). For a fixed vertex $\mathbf{a} \in V_\ell$, let the broken space

$$ \mathbf{R}\mathbf{T}\mathbf{N}_p(\mathcal{T}_\ell^a) := \{\mathbf{v}_\ell \in [L^2(\omega_\ell^a)]^d; \mathbf{v}_\ell|_K \in \mathbf{R}\mathbf{T}\mathbf{N}_p(K), \forall K \in \mathcal{T}_\ell^a\}, $$

where $\mathbf{R}\mathbf{T}\mathbf{N}_p(K) := [\mathbb{P}_p(K)]^d + \mathbb{P}_p(K)\mathbf{x}$ is the usual $p$-th order Raviart-Thomas-Nédélec space (cf. [14, 51]) on a simplex $K \in \mathcal{T}_\ell$. Then, the patchwise normal-trace-continuous spaces $\mathbf{V}_\ell^a$ with homogeneous Neumann boundary conditions in which the local equilibration will be performed are defined by

$$ \mathbf{V}_\ell^a := \begin{cases} \{\mathbf{v}_\ell \in \mathbf{R}\mathbf{T}\mathbf{N}_p^{\text{est}}(\mathcal{T}_\ell^a) \cap \mathbf{H}(\text{div}, \omega_\ell^a); \mathbf{v}_\ell \cdot \mathbf{n}_{\omega_\ell^a} = 0 \text{ on } \partial\omega_\ell^a\} & \text{if } \mathbf{a} \in \mathcal{V}_\ell^{\text{int}}, \\ \{\mathbf{v}_\ell \in \mathbf{R}\mathbf{T}\mathbf{N}_p^{\text{est}}(\mathcal{T}_\ell^a) \cap \mathbf{H}(\text{div}, \omega_\ell^a); \mathbf{v}_\ell \cdot \mathbf{n}_{\omega_\ell^a} = 0 \text{ on } \partial\omega_\ell^a \setminus \text{faces in } \partial\omega_\ell^a \cap \partial\Omega \text{ that share } \mathbf{a}\} & \text{if } \mathbf{a} \in \mathcal{V}_\ell^{\text{ext}}. \end{cases} \qquad (4.3) $$

Let $\mathbb{P}_p^{\text{est}}(\mathcal{T}_\ell^a)$ be composed of elementwise polynomials of degree $p_\ell^{\text{est}}$ on the patch $\mathcal{T}_\ell^a$. We will also need the $L^2(\omega_\ell^a)$-orthogonal projection onto $\mathbb{P}_p^{\text{est}}(\mathcal{T}_\ell^a)$ that we denote as $\Pi_\ell^a$.

**Definition 4.1** (Equilibrated flux $\sigma_\ell$ by local minimizations). Let $u_\ell^{\text{ex}}$ be the solution of (4.1). For each vertex $\mathbf{a} \in V_\ell$, let the local equilibrated flux $\sigma_\ell^a \in \mathbf{V}_\ell^a$ be defined by the local minimization problem

$$ \sigma_{\ell}^{\mathbf{a}} := \arg \min_{\mathbf{v}_{\ell} \in \mathbf{V}_{\ell}^{\mathbf{a}}} \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \mathbf{v}_{\ell}\|_{\omega_{\ell}^{\mathbf{a}}}. \qquad (4.4a) $$
$$ \nabla \cdot \mathbf{v}_{\ell} = \Pi_{\ell}^{\mathbf{a}} (\psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}}) $$

Then, extending each local contribution $\sigma_\ell^a$ by zero outside of the patch domain $\omega_\ell^a$, the global $\mathbf{H}(\text{div}, \Omega)$-conforming equilibrated flux $\sigma_\ell$ is constructed as

$$ \sigma_{\ell} := \sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \sigma_{\ell}^{\mathbf{a}}. \qquad (4.4b) $$

It will turn useful to denote

$$ g_{\ell}^{\mathbf{a}} := \psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}}. \qquad (4.5) $$

Then, the Neumann compatibility condition $\int_{\omega_\ell^a} g_\ell^\mathbf{a} = \int_{\omega_\ell^a} \psi_\ell^\mathbf{a} f - \nabla \psi_\ell^\mathbf{a} \cdot \mathbf{A} \nabla u_\ell^\text{ex} = 0$ for the problem (4.4a) is satisfied for all $\mathbf{a} \in \mathcal{V}_\ell^{\text{int}}$ as a direct consequence of (4.1) (consider $\psi_\ell^\mathbf{a}$ as a test function in (4.1)). The global $\mathbf{H}(\text{div}, \Omega)$-conformity of $\sigma_\ell$ follows from imposing the zero normal trace of functions in the local spaces $\mathbf{V}_\ell^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:

$$ \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}^{\text{ex}}) \| \leq \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}) := \left\{ \sum_{K \in \mathcal{T}_{\ell}} \eta_{K}^{2}(u_{\ell}^{\text{ex}}) \right\}^{\frac{1}{2}}, \quad \eta_{K} := \| \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell} \|_{K} + \frac{h_{K}}{\pi c_{\mathbf{A},K}^{1/2}} \| f - \nabla \cdot \boldsymbol{\sigma}_{\ell} \|_{K}, \quad (4.6) $$

where $c_{\mathbf{A},K}$ is the smallest eigenvalue of the diffusion tensor $\mathbf{A}$ on the simplex $K$. When $f$ is a piecewise polynomial in the space $\mathbb{P}_{\mathbf{p}_{\ell}-1}(\mathcal{T}_{\ell})$ and when $\mathbf{A}$ is piecewise constant, $\mathbf{A} \in \mathbb{P}_0(\mathcal{T}_{\ell})^{d \times d}$, the so-called data oscillation term $h_K \|f - \nabla \cdot \boldsymbol{\sigma}_{\ell}\|_K / (\pi c_{\mathbf{A},K}^{1/2})$ vanishes in the local estimator $\eta_K$. This follows from (4.8) below. It would also vanish if $f \in \mathbb{P}_{\mathbf{p}_{\ell}}(\mathcal{T}_{\ell})$ for uniform polynomial degrees $\mathbf{p}_{\ell} = \text{const}$, see (6.8) and the discussion in Remark 6.7 below.

## 4.2 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 $\mathcal{V}_{\ell}^{\theta} \subset \mathcal{V}_{\ell}$ using a bulk-chasing criterion inspired by the well-known Dörfler's criterion (cf. [29])

$$ \eta \left( u_{\ell}^{\text{ex}}, \bigcup_{\mathbf{a} \in \mathcal{V}_{\ell}^{\theta}} \mathcal{T}_{\ell}^{\mathbf{a}} \right) \geq \theta \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}), \quad (4.7) $$

where $\theta \in (0, 1]$ is a fixed threshold parameter, and where, for a subset $S \subset \mathcal{T}_{\ell}$, we adopt the notation $\eta(u_{\ell}^{\text{ex}}, S) := \{\sum_{K \in S} \eta_K(u_{\ell}^{\text{ex}})^2\}^{1/2}$. Then, letting

$$ \mathcal{M}_{\ell}^{\theta} := \bigcup_{\mathbf{a} \in \mathcal{V}_{\ell}^{\theta}} \mathcal{T}_{\ell}^{\mathbf{a}} \subset \mathcal{T}_{\ell} $$

be the collection of all the simplices that belong to a patch associated with a marked vertex, we observe that (4.7) means that $\eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta}) \geq \theta \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell})$. Let us denote by $\omega_{\ell} := \bigcup_{\mathbf{a} \in \mathcal{V}_{\ell}^{\theta}} \omega_{\ell}^{\mathbf{a}}$ the open subdomain corresponding to the set of marked elements $\mathcal{M}_{\ell}^{\theta}$. To select a set $\mathcal{V}_{\ell}^{\theta}$ of minimal cardinality, the mesh vertices in $\mathcal{V}_{\ell}$ are sorted by comparing the vertex-based error estimators $\eta(u_{\ell}, \mathcal{T}_{\ell}^{\mathbf{a}})$ for all $\mathbf{a} \in \mathcal{V}_{\ell}$, and a greedy algorithm is employed to build the set $\mathcal{V}_{\ell}^{\theta}$. 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 $\mathcal{V}_{\ell}^{\sharp}$ with the corresponding set of marked elements $\mathcal{M}_{\ell}^{\sharp}$ in the following way

$$ \mathcal{V}_{\ell}^{\sharp} := \bigcup_{K \in \mathcal{M}_{\ell}^{\theta}} \mathcal{V}_{\ell, K} \quad \text{and} \quad \mathcal{M}_{\ell}^{\sharp} := \bigcup_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} \mathcal{T}_{\ell}^{\mathbf{a}}. $$

In other words, we extend the set $\mathcal{M}_{\ell}^{\theta}$ by one more layer of neighbouring elements in contact with the boundary $\partial \omega_{\ell}$, see Figure 1 for an illustration. In addition, we define by $\omega_{\ell}^{\sharp} := \bigcup_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} \omega_{\ell}^{\mathbf{a}}$ the open subdomain corresponding to the set of marked elements $\mathcal{M}_{\ell}^{\sharp}$. This extension is motivated by the structure of the error estimate (4.6), stemming from equilibrated flux $\boldsymbol{\sigma}_{\ell}$ composed of local patchwise contributions $\boldsymbol{\sigma}_{\ell}^{\mathbf{a}}$. It plays a particular role later in our convergence proof, when we decompose the error estimate $\eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})$ into a sum of local patchwise contributions in (6.6) with their corresponding domains of definitions possibly exceeding the original marked subdomain $\omega_{\ell}$, but always included in the extended marked subdomain $\omega_{\ell}^{\sharp}$. 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_{\ell, \text{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, Proposition 4.2].
<figure><img src="image_5.png" /><figcaption>**Figure 1:** An example of local error estimators $\eta_K(u_\ell^{\text{ex}})$ from (4.6) (left) and illustration of the corresponding set of marked vertices $V_\ell^\theta = \{a_1\}$ and its extension $V_\ell^\sharp = \{a_1, a_2, \dots, a_8\}$, $\theta = 0.5$ (right). The region highlighted in the red color corresponds to the original marked subdomain $\omega_\ell$, while its union with the yellow region amounts to the extended marked subdomain $\omega_\ell^\sharp$.</figcaption></figure>

tion 3.1], and (6.6) below:

$$ \begin{gathered} \|\mathbf{A}^{\frac{1}{2}} \nabla(u - u_{\ell}^{\text{ex}})\| \le \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}) \le \sqrt{d+1} \left\{ \sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \eta_{\mathbf{a}}(u_{\ell}^{\text{ex}})^2 \right\}^{\frac{1}{2}}, \quad (4.8) \\ \eta_{\mathbf{a}}(u_{\ell}^{\text{ex}})^2 := \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \left[ \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{K} + \frac{h_K}{\pi c_{A,K}^{1/2}} \|g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}})\|_{K} \right]^2; \end{gathered} $$

since $\Pi_\ell^\mathbf{a}$ projects to $\mathbb{P}_{p_\ell, K}^\text{est}(\mathcal{T}_\ell^\mathbf{a})$, which contains the polynomials of order $p_\ell, K$ in each mesh element $K \in \mathcal{T}_\ell^\mathbf{a}$, and recalling (4.5), the term $\nabla\psi_\ell^\mathbf{a} \cdot \mathbf{A} \nabla u_\ell^{\text{ex}} - \Pi_\ell^\mathbf{a}(\nabla\psi_\ell^\mathbf{a} \cdot \mathbf{A} \nabla u_\ell^{\text{ex}})$ vanishes if the diffusion tensor is piecewise constant, $\mathbf{A} \in \mathbb{P}_0(\mathcal{T}_\ell)^{d \times d}$, whereas the term $\psi_\ell^\mathbf{a} f - \Pi_\ell^\mathbf{a}(\psi_\ell^\mathbf{a} f)$ vanishes if the source term is a piecewise polynomial, $f \in \mathbb{P}_{p_\ell-1}(\mathcal{T}_\ell)$. This would result in no need to extend the obtained set of marked vertices $V_\ell^\theta$ by the additional layer, as well as in an overall a simpler hp-adaptive algorithm based on the patchwise estimators $\eta_\mathbf{a}(u_\ell^{\text{ex}})$. However, unless a combination with (4.6) was used, the presence of the constant $\sqrt{d+1}$ already in the error estimate (4.8) would lead to a deterioration of the bound on the error reduction factor $C_{\ell, \text{red}}$, which was designed as sharp as possible in [25, 26].

## 4.3 The module REFINE

The module REFINE takes as input the extended set of marked vertices $V_\ell^\sharp$ and outputs the mesh $\mathcal{T}_{\ell+1}$ and the polynomial-degree distribution $\mathbf{p}_{\ell+1}$ to be used at the next iteration of the adaptive loop from Scheme 1.

<figure><img src="image_6.png" /><figcaption>**Figure 2:** An example of a marked patch $\mathcal{T}_\ell^\mathbf{a}$ together with its polynomial-degree distribution $\mathbf{p}_\ell^\mathbf{a}$ (left), its h-refinement (center), and p-refinement (right), the “adaptive residual lifting” strategy of Section 4.3.2.</figcaption></figure>
### 4.3.1 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 $\mathcal{T}_\ell^a$ assigned to a marked vertex $\mathbf{a} \in \mathcal{V}_\ell^\sharp$ using two distinct local primal solves. For this purpose, we will need a matching simplicial refinement $\mathcal{T}_\ell^{a,h}$ of $\mathcal{T}_\ell^a$ (see the center panel of Figure 2 for an example), as well as an increased polynomial-degree distribution vector $\mathbf{p}_\ell^{a,p}$ (see the right panel of Figure 2 for an example); let also $\mathbf{p}_\ell^{a,h}$ be the (original) polynomial-degree distribution in the patch $\mathcal{T}_\ell^a$ given by the vector $\mathbf{p}_\ell^a := \{p_{\ell,K}\}_{K \in \mathcal{T}_\ell^a}$ and $\mathcal{T}_\ell^{a,p} = \mathcal{T}_\ell^a$. For each vertex $\mathbf{a}$ from the extended set of marked vertices $\mathcal{V}_\ell^\sharp$, consider the local patch-based spaces $V_\ell^{a,h}$ and $V_\ell^{a,p}$ given by

$$V_{\ell}^{a,h} := \mathbb{P}_{\mathbf{p}_{\ell}^{a,h}}(\mathcal{T}_{\ell}^{a,h}) \cap H_{0}^{1}(\omega_{\ell}^{a}), \quad V_{\ell}^{a,p} := \mathbb{P}_{\mathbf{p}_{\ell}^{a,p}}(\mathcal{T}_{\ell}^{a,p}) \cap H_{0}^{1}(\omega_{\ell}^{a}). \qquad (4.9)$$

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} \in V_\ell^{a,h}$ and $r^{a,p} \in V_\ell^{a,p}$ be the h- and p-refinement residual liftings respectively given by

$$(\mathbf{A}\nabla r^{a,h}, \nabla v^{a,h})_{\omega_{\ell}^{a}} = (f, v^{a,h})_{\omega_{\ell}^{a}} - (\mathbf{A}\nabla u_{\ell}^{\text{ex}}, \nabla v^{a,h})_{\omega_{\ell}^{a}} \quad \forall v^{a,h} \in V_{\ell}^{a,h}, \qquad (4.10a)$$

$$(\mathbf{A}\nabla r^{a,p}, \nabla v^{a,p})_{\omega_{\ell}^{a}} = (f, v^{a,p})_{\omega_{\ell}^{a}} - (\mathbf{A}\nabla u_{\ell}^{\text{ex}}, \nabla v^{a,p})_{\omega_{\ell}^{a}} \quad \forall v^{a,p} \in V_{\ell}^{a,p}; \qquad (4.10b)$$

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” $\|\mathbf{A}^{\frac{1}{2}}\nabla r^{a,h}\|_{\omega_{\ell}^{a}}$ and $\|\mathbf{A}^{\frac{1}{2}}\nabla r^{a,p}\|_{\omega_{\ell}^{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.

### 4.3.2 “Adaptive residual lifting” refinement strategy; hp residual lifting

For each vertex $\mathbf{a}$ from the extended set of marked vertices $\mathcal{V}_\ell^\sharp$, 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 $\mathcal{T}_\ell^{a,h}$ by diving each simplex $K \in \mathcal{T}_\ell^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 $\mathbf{p}_\ell^{a,h}$ corresponding to $\mathcal{T}_\ell^{a,h}$ is then obtained from $\mathbf{p}_\ell^a$ by assigning to each newly-created simplex the same polynomial degree as that of its parent. For p-refinement testing, we set $\mathcal{T}_\ell^{a,p} = \mathcal{T}_\ell^a$ and the polynomial-degree distribution $\mathbf{p}_\ell^{a,p}$ is obtained from $\mathbf{p}_\ell^a$ by assigning to each simplex $K \in \mathcal{T}_\ell^{a,p}$ the polynomial degree $p_{\ell,K} + \delta_K^a$ where, as in [25, 26],

$$\delta_K^a := \begin{cases} 1 & \text{if } p_{\ell,K} = \min_{K' \in \mathcal{T}_\ell^a} p_{\ell,K'} \\ 0 & \text{otherwise.} \end{cases} \qquad (4.11)$$

This is illustrated in Figure 2, right, and only increases the smallest polynomial degree in the patch $\mathcal{T}_\ell^a$.

In difference with respect to [25, 26], though, we also immediately consider the suggested new pair $(\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+1})$ identified by the hp-decision on all vertices $\mathbf{a} \in \mathcal{V}_\ell^\sharp$ after step 20 of Algorithm 1, leading to the suggested finite element space $V_{\ell+1} := \mathbb{P}_{\mathbf{p}_{\ell+1}}(\mathcal{T}_{\ell+1}) \cap H_0^1(\Omega)$. Then we consider local subspaces of $V_{\ell+1}$, still on $\omega_\ell^a$, given by

$$V_{\ell}^{a, hp} := V_{\ell+1}|_{\omega_{\ell}^{a}} \cap H_{0}^{1}(\omega_{\ell}^{a}), \qquad (4.12)$$

and construct the patch residual lifting $r^{a, hp} \in V_\ell^{a, hp}$ by solving

$$(\mathbf{A}\nabla r^{a, hp}, \nabla v^{a, hp})_{\omega_{\ell}^{a}} = (f, v^{a, hp})_{\omega_{\ell}^{a}} - (\mathbf{A}\nabla u_{\ell}^{\text{ex}}, \nabla v^{a, hp})_{\omega_{\ell}^{a}} \quad \forall v^{a, hp} \in V_{\ell}^{a, hp}. \qquad (4.13)$$
<caption>**Algorithm 1 (module REFINE)**</caption>

1: **module REFINE** ($V_\ell^\sharp$)
2: &nbsp;&nbsp;&nbsp;▷ **Input:** extended set of marked vertices $V_\ell^\sharp$
3: &nbsp;&nbsp;&nbsp;▷ **Output:** new pair ($\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+1}$)
4: &nbsp;&nbsp;&nbsp;**for all** $\mathbf{a} \in V_\ell^\sharp$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Prepare the patch $h$-refinement space $V_\ell^{a,h}$ and the patch $p$-refinement space $V_\ell^{a,p}$ from (4.9) in dependence on the chosen strategy of either Section 4.3.2 or Section 4.3.3
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute the residual liftings $r^{a,h}, r^{a,p}$ given by (4.10a) and (4.10b).
8: &nbsp;&nbsp;&nbsp;**end for**
9: &nbsp;&nbsp;&nbsp;Set $V_\ell^p := \{\mathbf{a} \in V_\ell^\sharp \mid \|\mathbf{A}^{\frac{1}{2}}\nabla r^{a,h}\|_{\omega_\ell^a} < \|\mathbf{A}^{\frac{1}{2}}\nabla r^{a,p}\|_{\omega_\ell^a}\}$.
10: &nbsp;&nbsp;&nbsp;Set $V_\ell^h := V_\ell^\sharp \setminus V_\ell^p$.
11: &nbsp;&nbsp;&nbsp;Select $M_\ell^h := \{K \in \mathcal{T}_\ell \mid V_{\ell,K} \cap V_\ell^h \neq \emptyset\} \subset M_\ell^\sharp$.
12: &nbsp;&nbsp;&nbsp;Select $M_\ell^p := \{K \in \mathcal{T}_\ell \mid V_{\ell,K} \cap V_\ell^p \neq \emptyset\} \subset M_\ell^\sharp$.
13: &nbsp;&nbsp;&nbsp;Build $\mathcal{T}_{\ell+1}$ from $\mathcal{T}_\ell$ and $M_\ell^h$, such that $\mathcal{T}_\ell^{a,h} \subset \mathcal{T}_{\ell+1} \forall \mathbf{a} \in V_\ell^h$.
14: &nbsp;&nbsp;&nbsp;**for all** $K \in \mathcal{T}_{\ell+1}$ with its parent element $\tilde{K} \in \mathcal{T}_\ell$ **do**
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** “adaptive residual lifting” strategy of Section 4.3.2 **then**
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p_{\ell+1,K} := \begin{cases} p_{\ell,\tilde{K}} & \text{if } \tilde{K} \notin M_\ell^p, \\ p_{\ell,\tilde{K}} + \delta_{\tilde{K}}^a & \text{if } \tilde{K} \in M_\ell^p, \end{cases}$ where $\delta_{\tilde{K}} := \begin{cases} 1 & \text{if } p_{\ell,\tilde{K}} = \min_{\tilde{K}' \in \mathcal{T}_\ell^a} p_{\ell,\tilde{K}'}, \\ 0 & \text{otherwise.} \end{cases}$
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** “interior node, $p+d-1$” strategy of Section 4.3.3 **then**
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p_{\ell+1,K} := \begin{cases} \bar{p}_{\ell,\tilde{K}} & \text{if } \tilde{K} \notin M_\ell^p, \\ \max\{\bar{p}_{\ell,\tilde{K}} + d-1, p_{\ell,\tilde{K}} + d+1\} & \text{if } \tilde{K} \in M_\ell^p, \end{cases}$ where $\bar{p}_{\ell,\tilde{K}} := \max\{p_{\ell,\tilde{K}}, p_{\ell,\text{neigh},\tilde{K}}\}$.
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
20: &nbsp;&nbsp;&nbsp;**end for**
21: &nbsp;&nbsp;&nbsp;**if** “adaptive residual lifting” strategy of Section 4.3.2 **then**
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for all** $\mathbf{a} \in V_\ell^\sharp$ **do**
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute the hp residual lifting $r^{a,hp}$ given by (4.13).
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** the condition
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \left\{ \sum_{K \in \mathcal{T}_\ell^a} \left[ \|\psi_\ell^a \mathbf{A}^{\frac{1}{2}} \nabla u_\ell^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_\ell^a\|_K + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \|g_\ell^a - \Pi_\ell^a(g_\ell^a)\|_K \right]^2 \right\}^{\frac{1}{2}} \le C_{\text{lb}} \|\mathbf{A}^{\frac{1}{2}} \nabla r^{a,hp}\|_{\omega_\ell^a} \quad (4.14) $$ is not satisfied **then**
26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;re-assign $\mathcal{T}_\ell := \mathcal{T}_{\ell+1}$ and $\mathbf{p}_\ell := \mathbf{p}_{\ell+1}$ and go back to line 4.
27: &nbsp;&nbsp;&nbsp;**end for**
28: &nbsp;&nbsp;&nbsp;**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 ($\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+1}$). If (4.14) is not satisfied, we continue working on $\omega_\ell^a$ by further diving each simplex $K \in \mathcal{T}_\ell^a$ to create a more refined $\mathcal{T}_\ell^{a,h}$ and by increasing $\mathbf{p}_\ell^{a,p}$, until (4.14) is satisfied (which will be reached by virtue of Proposition 5.1 below, at least for small data oscillation).
### 4.3.3 “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 ($\mathcal{T}_{\ell+1}$, $\mathbf{p}_{\ell+1}$). Alternatively, when the polynomial degrees in $\mathbf{p}_\ell$ 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 $\omega_\ell^\mathbf{a}$, $\mathbf{a} \in \mathcal{V}_\ell^\mathbf{a}$, 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 $\mathbf{p}_f$ with respect to the coarsest partition $\mathcal{T}_0$, $f \in \mathbb{P}_{\mathbf{p}_f}(\mathcal{T}_0)$ for $\mathbf{p}_f \le \mathbf{p}_0 - 1$, as well as $\mathbf{A} \in \mathbb{P}_0(\mathcal{T}_0)^{d \times d}$ piecewise constant.

For each marked vertex $\mathbf{a} \in \mathcal{V}_\ell^\mathbf{a}$, we let $\mathcal{T}_\ell^{\mathbf{a},h}$ be obtained from $\mathcal{T}_\ell^\mathbf{a}$ by dividing each simplex $K \in \mathcal{T}_\ell^\mathbf{a}$ into at least six children simplices such that

$$ \text{a node interior to all the elements } K \in \mathcal{T}_\ell^\mathbf{a} \text{ as well as all the faces } F \in \mathcal{F}_\ell^{\mathbf{a},\text{int}} \text{ is created,} \quad (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} \le p_{\ell,K} - 2$ for all $K \in \mathcal{T}_\ell^\mathbf{a}$, then only the faces from $\mathcal{F}_\ell^{\mathbf{a},\text{int}}$ need to have an interior node but not the elements from $\mathcal{T}_\ell^\mathbf{a}$. The polynomial-degree distribution $\mathbf{p}_\ell^{\mathbf{a},h}$ corresponding to $\mathcal{T}_\ell^{\mathbf{a},h}$ is then copied from $\mathbf{p}_\ell^\mathbf{a}$, with a possible increase to match the polynomial degree of its marked neighbors. More precisely, for each $K \in \mathcal{T}_\ell^\mathbf{a}$, we set $\bar{p}_{\ell,K}$ as the maximum of the polynomial degrees $p_{\ell,K}$ and all $p_{\ell,K'}$ where $K' \in \mathcal{T}_\ell^\mathbf{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 $\mathcal{T}_\ell^{\mathbf{a},p} := \mathcal{T}_\ell^\mathbf{a}$ and define the polynomial-degree distribution vector $\mathbf{p}_\ell^{\mathbf{a},p}$ by assigning to each simplex $K \in \mathcal{T}_\ell^{\mathbf{a},p}$ the polynomial degree

$$ \max\{\bar{p}_{\ell,K} + d - 1, p_{f,K} + d + 1\}, \quad (4.16) $$

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^{\mathbf{a},hp}$ from (4.13) is only to be computed at a later stage in Section 4.3.4 below. 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><img src="image_7.png" /><figcaption>**Figure 3:** An example of a marked patch $\mathcal{T}_\ell^\mathbf{a}$ together with its polynomial-degree distribution $\mathbf{p}_\ell^\mathbf{a}$ (left), its h-refinement (center), and p-refinement (right), the “interior node, $p+d-1$” strategy of Section 4.3.3.</figcaption></figure>

### 4.3.4 Computable discrete lower bound

Once the new pair $(\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+1})$ is determined using one of the above strategies, the finite element space $V_{\ell+1} := \mathbb{P}_{\mathbf{p}_{\ell+1}}(\mathcal{T}_{\ell+1}) \cap H_0^1(\Omega)$ to be used on iteration $(\ell+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^{\mathbf{a},hp}$ by (4.13). After extending $r^{\mathbf{a},hp}$ from (4.13) by zero outside $\omega_\ell^\mathbf{a}$, we have for the current (available) approximation $u_\ell^{\text{ex}} \in V_\ell$ and the next level's (unavailable)
approximation $u_{\ell+1}^{\text{ex}} \in V_{\ell+1}$, by [25, Lemma 5.1] used with the extended set $V_{\ell}^{\sharp}$ in place of $V_{\ell}^{\theta}$,

$$ \| \mathbf{A}^{\frac{1}{2}} \nabla (u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}}) \|_{\omega_{\ell}^{\sharp}} \geq \underline{\eta}_{\mathcal{M}_{\ell}^{\sharp}}, \quad \underline{\eta}_{\mathcal{M}_{\ell}^{\sharp}} := \begin{cases} \frac{\sum_{a \in V_{\ell}^{\sharp}} \| \mathbf{A}^{\frac{1}{2}} \nabla^{a,hp} \|_{\omega_{\ell}^{a}}^{2}}{\| \mathbf{A}^{\frac{1}{2}} \nabla (\sum_{a \in V_{\ell}^{\sharp} r^{a,hp}) \|_{\omega_{\ell}^{\sharp}}}^{2}} & \text{if } \sum_{a \in V_{\ell}^{\sharp}} r^{a,hp} \neq 0, \\ 0 & \text{otherwise.} \end{cases} \quad (4.17) $$

Moreover, this lower bound can be further localized using the fact that each simplex has $(d+1)$ vertices as

$$ \underline{\eta}_{\mathcal{M}_{\ell}^{\sharp}} \geq \frac{\{ \sum_{a \in V_{\ell}^{\sharp}} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{a,hp} \|_{\omega_{\ell}^{a}}^{2} \}^{1/2}}{\sqrt{d+1}}. \quad (4.18) $$

Indeed, this can be seen from

$$ \left\| \mathbf{A}^{\frac{1}{2}} \nabla \left( \sum_{a \in V_{\ell}^{\sharp}} r^{a,hp} \right) \right\|_{\omega_{\ell}^{\sharp}}^{2} \leq \sum_{K \in \mathcal{M}_{\ell}^{\sharp}} (d+1) \sum_{a \in V_{\ell}^{\sharp} \cap V_K} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{a,hp} \|_{K}^{2} = (d+1) \sum_{a \in V_{\ell}^{\sharp}} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{a,hp} \|_{\omega_{\ell}^{a}}^{2}. $$

# 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 $\omega_{\ell}^{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.

## 5.1 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 $\mathbf{H}(\text{div})$ space

$$ \mathbf{V}^{\text{a}} := \begin{cases} \{\mathbf{v} \in \mathbf{H}(\text{div}, \omega_{\ell}^{\text{a}}); \mathbf{v}_{\ell} \cdot \mathbf{n}_{\omega_{\ell}^{\text{a}}} = 0 \text{ on } \partial \omega_{\ell}^{\text{a}}\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}}, \\ \{\mathbf{v} \in \mathbf{H}(\text{div}, \omega_{\ell}^{\text{a}}); \mathbf{v}_{\ell} \cdot \mathbf{n}_{\omega_{\ell}^{\text{a}}} = 0 \text{ on } \partial \omega_{\ell}^{\text{a}} \setminus \text{faces in } \partial \omega_{\ell}^{\text{a}} \cap \partial \Omega \text{ that share } \mathbf{a}\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{ext}}, \end{cases} \quad (5.1) $$

with the normal trace understood as usual by duality, and the infinite-dimensional constrained minimization

$$ \sigma^{\text{a}} := \arg \min_{\substack{\mathbf{v} \in \mathbf{V}^{\text{a}} \\ \nabla \cdot \mathbf{v} = \Pi_{\ell}(\psi_{\ell}^{\text{a}} f - \nabla \psi_{\ell}^{\text{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}})}} \| \psi_{\ell}^{\text{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \mathbf{v} \|_{\omega_{\ell}^{\text{a}}}; \quad (5.2) $$

note that (4.3)-(4.4a) are discrete versions of (5.1)-(5.2). The seminal p-robustness result by Braess et al. [12, Theorem 7] in two space dimensions and its extension to three space dimensions [35, Corollaries 3.3 and 3.8] show that

$$ \| \psi_{\ell}^{\text{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \sigma^{\text{a}} \|_{\omega_{\ell}^{\text{a}}} \leq C_{\text{st}} \| \psi_{\ell}^{\text{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \sigma^{\text{a}} \|_{\omega_{\ell}^{\text{a}}} \quad (5.3) $$

for a constant $C_{\text{st}}$ that only depends on the space dimension $d$, the mesh shape-regularity $\kappa_T$, and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient $\mathbf{A}$ on $\omega_{\ell}^{a}$. Computable upper bounds on $C_{\text{st}}$ are discussed in [34, Lemma 3.23 and Remark 3.24].

Let also

$$ H_{*}^{1}(\omega_{\ell}^{\text{a}}) := \begin{cases} \{v \in H^{1}(\omega_{\ell}^{\text{a}}), (v, 1)_{\omega_{\ell}^{\text{a}}} = 0\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}}, \\ \{v \in H^{1}(\omega_{\ell}^{\text{a}}), v = 0 \text{ on faces in } \partial \omega_{\ell}^{\text{a}} \cap \partial \Omega \text{ that share } \mathbf{a}\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{ext}}. \end{cases} \quad (5.4) $$

For all $v \in H_{*}^{1}(\omega_{\ell}^{\text{a}})$, there holds $\psi_{\ell}^{\text{a}}v \in H_{0}^{1}(\omega_{\ell}^{\text{a}})$. Moreover, as in Carstensen and Funken [20, Theorem 3.1], Braess et al. [12, Section 3], [34, estimate (3.29)], or [11, Lemma 2.4], there holds

$$ \| \mathbf{A}^{\frac{1}{2}} \nabla (\psi_{\ell}^{\text{a}} v) \|_{\omega_{\ell}^{\text{a}}} \leq C_{\text{cont,PF}} \| \mathbf{A}^{\frac{1}{2}} \nabla v \|_{\omega_{\ell}^{\text{a}}}, \quad (5.5) $$
where the constant $C_{\text{cont,PF}}$ only depends on the space dimension $d$, the mesh shape-regularity $\kappa_{\mathcal{T}}$, and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient $\mathbf{A}$ on $\omega_{\ell}^{\mathbf{a}}$; $C_{\text{cont,PF}}$ relies on the Poincaré–Friedrichs inequality and is fully computable. Finally, denote

$$ \eta_{\text{osc},\mathbf{a}} := \left\{ \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \left[ \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \|g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}})\|_{K} \right]^2 \right\}^{\frac{1}{2}}, \quad (5.6) $$

where we recall that $g_{\ell}^{\mathbf{a}}$ is given by (4.5) and that $c_{\mathbf{A},K}$ is the smallest eigenvalue of the diffusion tensor $\mathbf{A}$ on the simplex $K$. Also recall from Remark 4.2 that $\nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}} - \Pi_{\ell}^{\mathbf{a}}(\nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}})$ vanishes if the diffusion tensor is piecewise constant, $\mathbf{A} \in \mathbb{P}_0 (\mathcal{T}_{\ell}^{\mathbf{a}})^{d \times d}$, whereas the term $\psi_{\ell}^{\mathbf{a}} f - \Pi_{\ell}^{\mathbf{a}}(\psi_{\ell}^{\mathbf{a}} f)$ vanishes if the source term is a piecewise polynomial, $f \in \mathbb{P}_{\mathbf{P}_{\ell}-1}(\mathcal{T}_{\ell}^{\mathbf{a}})$.

**Proposition 5.1** (Discrete local p-robust stability of the flux equilibration). Let $u_{\ell}^{\text{ex}} \in V_{\ell}$ satisfy the hat function orthogonality

$$ (f, \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A} \nabla u_{\ell}^{\text{ex}}, \nabla \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} = 0 \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}} $$

and let the local equilibrated flux $\sigma_{\ell}^{\mathbf{a}}$ be constructed by (4.4a). Let data oscillation do not dominate in that

$$ \eta_{\text{osc},\mathbf{a}} \le \frac{1}{(1+2C_{\text{st}})} \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \sigma_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}. \quad (5.7) $$

Then, there exists a constant $C_{\text{lb}} \ge 1$ only depending on the space dimension $d$, the mesh shape-regularity $\kappa_{\mathcal{T}}$, and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient $\mathbf{A}$ so that, for a sufficiently fine patch space $V_{\ell}^{\mathbf{a},\text{hp}}$ from (4.12) and $r^{\mathbf{a},\text{hp}}$ given by (4.13), there holds

$$ \left\{ \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \left[ \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \sigma_{\ell}^{\mathbf{a}} \|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \| g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}}) \|_{K} \right]^2 \right\}^{\frac{1}{2}} \le C_{\text{lb}} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a},\text{hp}} \|_{\omega_{\ell}^{\mathbf{a}}} \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}. \quad (5.8) $$

More precisely, it is possible to take

$$ C_{\text{lb}} := 4C_{\text{st}}C_{\text{cont,PF}}, \quad (5.9) $$

where the constants $C_{\text{st}}$ and $C_{\text{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

$$ \left\{ \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \left[ \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \sigma_{\ell}^{\mathbf{a}} \|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \| g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}}) \|_{K} \right]^2 \right\}^{\frac{1}{2}} \le \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \sigma_{\ell}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} + \eta_{\text{osc},\mathbf{a}}. \quad (5.10) $$

**Step 2.** In addition to $\sigma^{\mathbf{a}}$, let us denote by $\tilde{\sigma}^{\mathbf{a}}$ the solution to problem (5.2) where the polynomial projection $\Pi_{\ell}^{\mathbf{a}}(\psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}}) = \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}})$ is replaced by $\psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}} = g_{\ell}^{\mathbf{a}}$. Treating the mismatch between these two terms as in [34, proof of Theorem 3.17], one uncovers from (5.3) that

$$ \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \sigma_{\ell}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} \le C_{\text{st}} \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \tilde{\sigma}_{\ell}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} + C_{\text{st}} \eta_{\text{osc},\mathbf{a}}. \quad (5.11) $$

Moreover, the primal-dual equivalence as in [34, Remark 3.15] or [35, Corollary 3.6] shows that

$$ \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \tilde{\sigma}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} = \max_{v \in H_{*}^{1}(\omega_{\ell}^{\mathbf{a}})} \frac{(f, \psi_{\ell}^{\mathbf{a}} v)_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A} \nabla u_{\ell}^{\text{ex}}, \nabla (\psi_{\ell}^{\mathbf{a}} v))_{\omega_{\ell}^{\mathbf{a}}}}{\| \mathbf{A}^{\frac{1}{2}} \nabla v \|_{\omega_{\ell}^{\mathbf{a}}}}. \quad (5.12) $$

**Step 3.** From (5.5), we can rewrite and bound the right-hand side of (5.12) as

$$ \begin{aligned} & \max_{v \in H_{*}^{1}(\omega_{\ell}^{\mathbf{a}})} \frac{(f, \psi_{\ell}^{\mathbf{a}} v)_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A} \nabla u_{\ell}^{\text{ex}}, \nabla (\psi_{\ell}^{\mathbf{a}} v))_{\omega_{\ell}^{\mathbf{a}}}}{\| \mathbf{A}^{\frac{1}{2}} \nabla (\psi_{\ell}^{\mathbf{a}} v) \|_{\omega_{\ell}^{\mathbf{a}}}} \le C_{\text{cont,PF}} \max_{v \in H_{0}^{1}(\omega_{\ell}^{\mathbf{a}})} \frac{(f, v)_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A} \nabla u_{\ell}^{\text{ex}}, \nabla v)_{\omega_{\ell}^{\mathbf{a}}}}{\| \mathbf{A}^{\frac{1}{2}} \nabla v \|_{\omega_{\ell}^{\mathbf{a}}}} \quad (5.13) \\ & = C_{\text{cont,PF}} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}}, \end{aligned} $$
where $r^a \in H_0^1(\omega_\ell^a)$ is given by

$$ (\mathbf{A}\nabla r^a, \nabla v)_{\omega_\ell^a} = (f, v)_{\omega_\ell^a} - (\mathbf{A}\nabla u_\ell^{\text{ex}}, \nabla v)_{\omega_\ell^a} \quad \forall v \in H_0^1(\omega_\ell^a). \quad (5.14) $$

Step 4. From the above steps, we have

$$ \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} \le C_{\text{st}} C_{\text{cont,PF}} \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} + C_{\text{st}} \eta_{\text{osc,a}}. $$

From here, condition (5.7) on $\eta_{\text{osc,a}}$ yields

$$ (1 + C_{\text{st}}) \eta_{\text{osc,a}} \le \frac{(1 + C_{\text{st}})}{(1 + 2C_{\text{st}})} \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} \le C_{\text{st}} C_{\text{cont,PF}} \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}. $$

Thus,

$$ \begin{aligned} \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} + \eta_{\text{osc,a}} &\le C_{\text{st}} C_{\text{cont,PF}} \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} + (1 + C_{\text{st}}) \eta_{\text{osc,a}} \\ &\le 2C_{\text{st}} C_{\text{cont,PF}} \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}. \end{aligned} $$

Step 5. Now, recall that the patch $\omega_\ell^a$ and the data $f$ and $u_\ell^{\text{ex}}$ are fixed, whereas $r^a$ and $r^{a, hp}$ are prescribed therefrom by respectively (5.14) and (4.13). Moreover, the space $V_\ell^{a, hp}$ from (4.12) for computing $r^{a, hp}$ (local to the patch $\omega_\ell^a$ and discrete) is created by applying $h$- and/or $p$-refinements in the patch $\omega_\ell^a$. Thus, there holds

$$ \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}, hp}\|_{\omega_{\ell}^{\mathbf{a}}} \to \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}, \quad (5.15) $$

and there exists a sufficiently fine patch space $V_\ell^{a, hp}$ such that $\|\mathbf{A}^{\frac{1}{2}} \nabla r^{a, hp}\|_{\omega_\ell^a} \ge \|\mathbf{A}^{\frac{1}{2}} \nabla r^a\|_{\omega_\ell^a}/2$, so that (5.8) holds with the constant $C_{\text{lb}}$ given by (5.9). $\square$

## 5.2 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_{\text{lb}}$ possibly dependent on the polynomial degree $p$. We suppose here $f$ piecewise $(\mathbf{p}_0 - 1)$-polynomial and $\mathbf{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 \in \mathbb{P}_{\mathbf{p}_f}(\mathcal{T}_0)$ fulfill $\mathbf{p}_f \le \mathbf{p}_0 - 1$ and let the diffusion tensor $\mathbf{A}$ be piecewise constant with respect to $\mathcal{T}_0$, $\mathbf{A} \in \mathbb{P}_0(\mathcal{T}_0)^{d \times d}$. Let $u_\ell^{\text{ex}} \in V_\ell$ satisfy the hat function orthogonality

$$ (f, \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A} \nabla u_{\ell}^{\text{ex}}, \nabla \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} = 0 \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}} $$

and let the local equilibrated flux $\boldsymbol{\sigma}_\ell^a$ be constructed by (4.4a). Finally, let the patch space $V_\ell^{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_{\text{lb}} \ge 1$ only depending on the space dimension $d$, the mesh shape-regularity $\kappa_\tau$, the ratio of the largest and the smallest eigenvalue of the diffusion coefficient $\mathbf{A}$, and the maximal polynomial degree $p_{\text{max}}$ such that

$$ \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} \le C_{\text{lb}} \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}, hp}\|_{\omega_{\ell}^{\mathbf{a}}} \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}. \quad (5.16) $$

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 $\mathbf{A}$ being an identity matrix for simplicity; extension to piecewise constant $\mathbf{A}$ is immediate. From now on, we use the shorthand notation $x_1 \lesssim x_2$ when there exists a generic positive constant $C$ that only depends on the space dimension $d$, the shape-regularity $\kappa_\tau$ of the underlying hierarchy of meshes, and the polynomial degree of approximation employed locally such that $x_1 \le C x_2$.

**Lemma 5.3** (Discrete stability of the element residuals). Let $f \in \mathbb{P}_{\mathbf{p}_f}(\mathcal{T}_0)$ with $\mathbf{p}_f \le \mathbf{p}_0 - 1$ and $\mathbf{A} = \mathbf{I}$. Let a vertex $\mathbf{a} \in \mathcal{V}_\ell^\sharp$, the patch space $V_\ell^{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

$$ h_K \|f + \Delta u_{\ell}^{\text{ex}}\|_{K} \lesssim \|\nabla r^{a, hp}\|_{K} \quad \forall K \in \mathcal{T}_{\ell}^{a}. \quad (5.17) $$
Proof. Fix an element $K \in \mathcal{T}_{\ell}^{\mathbf{a}}$ and set

$$v_K := (f + \Delta u_{\ell}^{\text{ex}})|_K. \qquad (5.18)$$

We note that $v_K \in \mathbb{P}_{\max\{p_{\ell,K}-2,p_{f,K}\}}(K)$, using the assumption $f \in \mathbb{P}_{p_f}(\mathcal{T}_0)$. Let us define $\psi_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

$$\psi_K \in \begin{cases} \mathbb{P}_1(\mathcal{T}_{\ell+1}|_K) \cap H_0^1(K) & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^h, \\ \mathbb{P}_{d+1}(K) \cap H_0^1(K) & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^p, \end{cases} \quad \|\psi_K\|_{\infty,K} = 1; \qquad (5.19)$$

more precisely, in the first case, $\psi_K$ is a piecewise affine function taking nonzero values in the(/all) interior nodes in $K$, whereas in the second case, $\psi_K$ is the usual element bubble function. We crucially remark that $\psi_K v_K$ extended by zero is contained in the local hp-refinement space $V_{\ell}^{\mathbf{a},hp}$ due to (4.12) with either (4.15) or (4.16). In particular, either h-refinement has been employed and then $\mathcal{T}_{\ell+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_{\ell,K} + d - 1, p_{f,K} + d + 1\}$.

Now, the equivalence of norms on finite-dimensional spaces gives

$$(v_K, v_K)_K \lesssim (v_K, \psi_K v_K)_K, \qquad (5.20)$$

where the hidden constant depends on the polynomial degrees. This in combination with (5.18) yields

$$\|v_K\|_K^2 = (v_K, v_K)_K \lesssim (v_K, \psi_K v_K)_K = (f, \psi_K v_K)_K + (\Delta u_{\ell}^{\text{ex}}, \psi_K v_K)_K.$$

Then, noting that $(\nabla u_{\ell}^{\text{ex}})|_K \in \mathbf{H}(\text{div}, K)$ and $\psi_K v_K \in H_0^1(K)$, employing the Green theorem and using (4.13) with $\mathbf{A} = \mathbf{I}$ and the Cauchy-Schwarz inequality, we obtain

$$\begin{aligned} \|v_K\|_K^2 & \lesssim (f, \psi_K v_K)_K - (\nabla u_{\ell}^{\text{ex}}, \nabla(\psi_K v_K))_K = (\nabla r^{\mathbf{a},hp}, \nabla(\psi_K v_K))_K, \\ & \le \|\nabla r^{\mathbf{a},hp}\|_K \|\nabla(\psi_K v_K)\|_K. \end{aligned} \qquad (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

$$\|\nabla(\psi_K v_K)\|_K \lesssim h_K^{-1} \|\psi_K v_K\|_K. \qquad (5.22)$$

Moreover, from the definition of the bubble function $\psi_K$, there holds

$$\|\psi_K v_K\|_K \le \|\psi_K\|_{\infty,K} \|v_K\|_K = \|v_K\|_K. \qquad (5.23)$$

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 $\mathcal{F}_{\ell}^{\mathbf{a},\text{int}}$ have an interior node and either $f=0$ or $p_{f,K} \le p_{\ell,K} - 2$ for all $K \in \mathcal{T}_{\ell}^{\mathbf{a}}$, one rather uses $\psi_K \in \mathbb{P}_2(\mathcal{T}_{\ell+1}|_K) \cap H_0^1(K)$ with $\|\psi_K\|_{\infty,K} = 1$, the edge bubble function of a newly-added edge, in place of the first choice in (5.19). $\square$

**Lemma 5.4** (Discrete stability of the face residuals). Let $f \in \mathbb{P}_{p_f}(\mathcal{T}_0)$ with $p_f \le p_0 - 1$ and $\mathbf{A} = \mathbf{I}$. Let a vertex $\mathbf{a} \in \mathcal{V}_{\ell}^{\mathbf{a}}$, the patch space $V_{\ell}^{\mathbf{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^{\mathbf{a},hp}$ from (4.13) be given. Let $\mathcal{T}_F$ denote the set of elements containing the two elements $K \in \mathcal{T}_{\ell}^{\mathbf{a}}$ that share a face $F \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}$, and let $\omega_F$ be the corresponding open subdomain. Then

$$h_F^{\frac{1}{2}} \|\nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_F\|_F \lesssim \|\nabla r^{\mathbf{a},hp}\|_{\omega_F} \quad \forall F \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}.$$

Proof. Fix an edge $F \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}$. This time, we set

$$v_F := \|\nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_F\|_F \qquad (5.24)$$
and note that $v_F \in \mathbb{P}_{\max\{p_{\ell,K}-1, p_{\ell,K'}-1\}}$, where $K, K' \in \mathcal{T}_{\ell}^{\mathbf{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 $\psi_F$ be the bubble function on $\mathcal{T}_F$, depending on whether $h$- or $p$-refinement was applied within the **REFINE** module in Section 4.3.3; namely, for the vertex $\mathbf{a} \in \mathcal{V}_{\ell}^{\mathbf{h}}$,

$$ \psi_F \in \left\{ \begin{array}{ll} \mathbb{P}_1 (\mathcal{T}_{\ell+1} | \mathcal{T}_F) \cap H_0^1 (\omega_F) & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\mathbf{h}}, \\ \mathbb{P}_d (\mathcal{T}_F) \cap H_0^1 (\omega_F) & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\mathbf{p}}, \end{array} \right. \quad \|\psi_F\|_{\infty, \omega_F} = 1; \qquad (5.25) $$

similarly to (5.19), $\psi_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 $\psi_F v_F$ extended by zero lies in the local $hp$-refinement space $V_{\ell}^{\mathbf{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,

$$ (v_F, v_F)_F \lesssim (v_F, \psi_F v_F)_F. \qquad (5.26) $$

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 $\mathcal{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

$$ \|v_F\|_{\omega_F} \lesssim h_F^{\frac{1}{2}} \|v_F\|_F. \qquad (5.27) $$

Employing (5.26), (5.24), and expanding the jump term with $\mathcal{T}_F := \{K, K'\}$ and the convention that $\mathbf{n}_F$ points from $K$ to $K'$, we have

$$ \begin{aligned} \|v_F\|_F^2 &= (v_F, v_F)_F \lesssim (v_F, \psi_F v_F)_F = ((\nabla u_{\ell}^{\text{ex}}|_K)|_F \cdot \mathbf{n}_F, \psi_F v_F)_F - ((\nabla u_{\ell}^{\text{ex}}|_{K'})|_F \cdot \mathbf{n}_F, \psi_F v_F)_F \\ &= ((\nabla u_{\ell}^{\text{ex}}|_K)|_{\partial K} \cdot \mathbf{n}_{\partial K}, \psi_F v_F)_{\partial K} + ((\nabla u_{\ell}^{\text{ex}}|_{K'})|_{\partial K' \cdot \mathbf{n}_{\partial K'}}, \psi_F v_F)_{\partial K'}. \end{aligned} \qquad (5.28) $$

We have also used the fact that $(\psi_F v_F)|_{\partial \omega_F} = 0$ due to the definition of the bubble function $\psi_F$. Moreover, for each simplex $K \in \mathcal{T}_F$, $(\psi_F v_F)|_K \in H^1(K)$ and $(\nabla u_{\ell}^{\text{ex}})|_K \in \mathbf{H}(\text{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, \psi_F v_F)_{\omega_F}$, using that $\psi_F v_F$ extended by zero lies in $V_{\ell}^{\mathbf{a},hp}$, and recalling the discrete problem (4.13) with $\mathbf{A} = \mathbf{I}$, we finally obtain

$$ \begin{aligned} \|v_F\|_F^2 &\lesssim (\nabla u_{\ell}^{\text{ex}}, \nabla(\psi_F v_F))_{\omega_F} - (f, \psi_F v_F)_{\omega_F} + (f + \Delta u_{\ell}^{\text{ex}}, \psi_F v_F)_{\omega_F} \\ &= -(\nabla r^{\mathbf{a},hp}, \nabla(\psi_F v_F))_{\omega_F} + (f + \Delta u_{\ell}^{\text{ex}}, \psi_F v_F)_{\omega_F}. \end{aligned} \qquad (5.29) $$

From the definition of the bubble function $\psi_F$, there holds

$$ \|\psi_F v_F\|_{\omega_F} \le \|\psi_F\|_{\infty, \omega_F} \|v_F\|_{\omega_F} = \|v_F\|_{\omega_F}, \qquad (5.30) $$

while the inverse inequality and the shape-regularity of the mesh $\mathcal{T}_\ell$ give, again with a $p$-dependent constant,

$$ \|\nabla(\psi_F v_F)\|_{\omega_F} \lesssim h_F^{-1} \|\psi_F v_F\|_{\omega_F}. \qquad (5.31) $$

Then, the following chain of inequalities holds true due to the Cauchy-Schwarz inequality, (5.29), (5.30), (5.31), and (5.27):

$$ \begin{aligned} \|v_F\|_F^2 &\lesssim \|\nabla r^{\mathbf{a},hp}\|_{\omega_F} \|\nabla(\psi_F v_F)\|_{\omega_F} + \|f + \Delta u_{\ell}^{\text{ex}}\|_{\omega_F} \|\psi_F v_F\|_{\omega_F} \\ &\lesssim \|\nabla r^{\mathbf{a},hp}\|_{\omega_F} h_F^{-1} \|v_F\|_{\omega_F} + \|f + \Delta u_{\ell}^{\text{ex}}\|_{\omega_F} \|v_F\|_{\omega_F} \\ &\lesssim (\|\nabla r^{\mathbf{a},hp}\|_{\omega_F} + h_F \|f + \Delta u_{\ell}^{\text{ex}}\|_{\omega_F}) h_F^{-\frac{1}{2}} \|v_F\|_F. \end{aligned} $$

Now, the assertion of the lemma follows from the above result combined with (5.17) and the shape-regularity of the mesh $\mathcal{T}_\ell$. $\square$
Proof of Proposition 5.2. Let $\mathbf{a} \in \mathcal{V}_\ell^\sharp$ be fixed and recall we consider here $\mathbf{A} = \mathbf{I}$ and $f \in \mathbb{P}_{\mathbf{p}_f}(\mathcal{T}_0)$ with $\mathbf{p}_f \le \mathbf{p}_0 - 1$, so that $\Pi_\ell^\mathbf{a}(g_\ell^\mathbf{a}) = \Pi_\ell^\mathbf{a}(\psi_\ell^\mathbf{a}f - \nabla\psi_\ell^\mathbf{a} \cdot \mathbf{A}\nabla u_\ell^\mathbf{a}) = \psi_\ell^\mathbf{a}f - \nabla\psi_\ell^\mathbf{a} \cdot \mathbf{A}\nabla u_\ell^\mathbf{a} = g_\ell^\mathbf{a}$. We will crucially employ two reformulations of the local minimization problem (4.4a). Let the patchwise discontinuous piecewise polynomial spaces $\mathcal{Q}_\ell^\mathbf{a}$ be defined by

$$Q_\ell^\mathbf{a} := \begin{cases} \{q_\ell \in \mathbb{P}_{p_\ell}^{\text{est}}(\mathcal{T}_\ell^\mathbf{a}); (q_\ell, 1)_{\omega_\ell^\mathbf{a}} = 0\} & \text{if } \mathbf{a} \in \mathcal{V}_\ell^{\text{int}}, \\ \mathbb{P}_{p_\ell}^{\text{est}}(\mathcal{T}_\ell^\mathbf{a}) & \text{if } \mathbf{a} \in \mathcal{V}_\ell^{\text{ext}}. \end{cases} \quad (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 $(\sigma_\ell^\mathbf{a}, \delta_\ell^\mathbf{a}) \in \mathbf{V}_\ell^\mathbf{a} \times \mathcal{Q}_\ell^\mathbf{a}$ such that

$$(\sigma_{\ell}^{\mathbf{a}}, \mathbf{v}_{\ell})_{\omega_{\ell}^{\mathbf{a}}} - (\delta_{\ell}^{\mathbf{a}}, \nabla \cdot \mathbf{v}_{\ell})_{\omega_{\ell}^{\mathbf{a}}} = -(\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}}, \mathbf{v}_{\ell})_{\omega_{\ell}^{\mathbf{a}}} \quad \forall \mathbf{v}_{\ell} \in \mathbf{V}_{\ell}^{\mathbf{a}}, \quad (5.33\text{a})$$

$$(\nabla \cdot \sigma_{\ell}^{\mathbf{a}}, q_{\ell})_{\omega_{\ell}^{\mathbf{a}}} = (\psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \nabla u_{\ell}^{\text{ex}}, q_{\ell})_{\omega_{\ell}^{\mathbf{a}}} \quad \forall q_{\ell} \in \mathcal{Q}_{\ell}^{\mathbf{a}}. \quad (5.33\text{b})$$

Next, imposing the normal trace continuity constraint on faces from $\mathcal{F}_\ell^{\mathbf{a},\text{int}}$ and the no-flux condition on all faces from $\mathcal{F}_\ell^{\mathbf{a},\text{ext}}$ for $\mathbf{a} \in \mathcal{V}_\ell^\text{int}$ and only faces from $\mathcal{F}_\ell^{\mathbf{a},\text{ext}}$ not sharing the vertex $\mathbf{a} \in \mathcal{V}_\ell^\text{ext}$ via a Lagrange multiplier leads to the hybridized formulation: seek $\sigma_\ell^\mathbf{a}$ in the broken space $\mathbf{RTN}_{p_\ell}^{\text{est}}(\mathcal{T}_\ell^\mathbf{a})$, $\delta_\ell^\mathbf{a} \in \mathcal{Q}_\ell^\mathbf{a}$, and $\lambda_\ell^F \in \mathbb{P}_{p_\ell}^{\text{est}}(F)$ for all $F \in \mathcal{F}_\ell^\mathbf{a} \setminus \mathcal{F}_\ell^{\mathbf{a},\partial\Omega}$ such that

$$\sum_{K \in \mathcal{T}_\ell^\mathbf{a}} (\psi_\ell^\mathbf{a} \nabla u_\ell^\text{ex} + \sigma_\ell^\mathbf{a}, \mathbf{v}_\ell)_K - \sum_{K \in \mathcal{T}_\ell^\mathbf{a}} (\nabla \cdot \mathbf{v}_\ell, \delta_\ell^\mathbf{a})_K + \sum_{K \in \mathcal{T}_\ell^\mathbf{a}} \sum_{F \in \mathcal{F}_\ell, K \in \mathcal{F}_\ell^{\mathbf{a},\partial\Omega}} (\mathbf{v}_\ell \cdot \mathbf{n}_K, \lambda_\ell^F)_F = 0 \\ \forall \mathbf{v}_\ell \in \mathbf{RTN}_{p_\ell}^{\text{est}}(\mathcal{T}_\ell^\mathbf{a}), \quad (5.34\text{a})$$

$$\sum_{K \in \mathcal{T}_\ell^\mathbf{a}} (\nabla \cdot \sigma_\ell^\mathbf{a}, q_\ell)_K = \sum_{K \in \mathcal{T}_\ell^\mathbf{a}} (\psi_\ell^\mathbf{a} f - \nabla \psi_\ell^\mathbf{a} \cdot \nabla u_\ell^\text{ex}, q_\ell)_K \quad \forall q_\ell \in \mathcal{Q}_\ell^\mathbf{a}, \quad (5.34\text{b})$$

$$- \sum_{K \in \mathcal{T}_F} (\sigma_{\ell}^{\mathbf{a}} \cdot \mathbf{n}_K, \xi_{\ell})_F = 0 \qquad \forall \xi_{\ell} \in \mathbb{P}_{p_\ell}^{\text{est}}(F), \forall F \in \mathcal{F}_\ell^\mathbf{a} \setminus \mathcal{F}_\ell^{\mathbf{a},\partial\Omega}, \quad (5.34\text{c})$$

where we used the notation

$$\mathcal{F}_{\ell}^{\mathbf{a},\partial\Omega} := \begin{cases} \emptyset & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}}, \\ \{F \in \mathcal{F}_{\ell}^{\mathbf{a},\text{ext}}; F \text{ shares } \mathbf{a}\} & \text{if } \mathbf{a} \in \mathcal{V}_{\ell}^{\text{ext}}. \end{cases} \quad (5.35)$$

We note that $\psi_\ell^\mathbf{a} \nabla u_\ell^\text{ex} + \sigma_\ell^\mathbf{a} \in \mathbf{RTN}_{p_\ell}^{\text{est}}(\mathcal{T}_\ell^\mathbf{a})$, as (4.3) uses the maximal local polynomial degree $p_\ell^{\text{est}}$. Thus, we can use it as a test function $\mathbf{v}_\ell$ in (5.34a). Afterwards, employing $\delta_\ell^\mathbf{a}$ as a test function $q_\ell$ in (5.34b) and $\lambda_\ell^F$ as $\xi_\ell$ in (5.34c), summing (5.34c) over all $F \in \mathcal{F}_\ell^\mathbf{a} \setminus \mathcal{F}_\ell^{\mathbf{a},\partial\Omega}$, and finally summing (5.34a)-(5.34c), we obtain

$$\begin{aligned} \| \psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}}^{2} &= \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} (\psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \nabla u_{\ell}^{\text{ex}} + \nabla \cdot (\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}}), \delta_{\ell}^{\mathbf{a}})_K - \sum_{F \in \mathcal{F}_{\ell}^{\mathbf{a}} \setminus \mathcal{F}_{\ell}^{\mathbf{a},\partial\Omega}} (\llbracket \psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_F \rrbracket, \lambda_{\ell}^{F})_F \\ &= \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} (\psi_{\ell}^{\mathbf{a}}(f + \Delta u_{\ell}^{\text{ex}}), \delta_{\ell}^{\mathbf{a}})_K - \sum_{F \in \mathcal{F}_{\ell}^{\mathbf{a}} \setminus \mathcal{F}_{\ell}^{\mathbf{a},\partial\Omega}} (\psi_{\ell}^{\mathbf{a}}[\nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_F], \lambda_{\ell}^{F})_F \quad (5.36) \\ &\le \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \| \psi_{\ell}^{\mathbf{a}}(f + \Delta u_{\ell}^{\text{ex}}) \|_{K} \| \delta_{\ell}^{\mathbf{a}} \|_{K} + \sum_{F \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}} \| \psi_{\ell}^{\mathbf{a}}[\nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_F] \|_{F} \| \lambda_{\ell}^{F} \|_{F}. \end{aligned}$$

For the estimate on the right-hand side of (5.36), we have also used the fact that the hat function $\psi_\ell^\mathbf{a}|_F = 0$ $\forall F \in \mathcal{F}_\ell^{\mathbf{a},\text{ext}}$ for vertices $\mathbf{a} \in \mathcal{V}_\ell^\text{int}$, so that the boundary faces are excluded from the sum, and $\psi_\ell^\mathbf{a}|_F = 0$ $\forall F \in \mathcal{F}_\ell^{\mathbf{a},\text{ext}} \setminus \mathcal{F}_\ell^{\mathbf{a},\partial\Omega}$ for $\mathbf{a} \in \mathcal{V}_\ell^\text{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

$$\| \delta_{\ell}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} \lesssim h_{\omega_{\ell}^{\mathbf{a}}} \| \psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}}. \quad (5.37)$$
Secondly, $\|\psi_{\ell}^{\mathbf{a}}\|_{\infty,K} \le 1$ together with Lemma 5.3 yield the estimate

$$ \|\psi_{\ell}^{\mathbf{a}}(f + \Delta u_{\ell}^{\text{ex}})\|_{K} \le \|f + \Delta u_{\ell}^{\text{ex}}\|_{K} \lesssim h_{K}^{-1} \|\nabla r^{\mathbf{a},\text{hp}}\|_{K} \quad \forall K \in \mathcal{T}_{\ell}^{\mathbf{a}}. \quad (5.38) $$

Combining (5.36) with (5.38), the Cauchy–Schwarz inequality, the shape-regularity leading to $h_K \approx h_{\omega_\ell^\mathbf{a}}$, and (5.37), we obtain the estimate on the first term

$$ \begin{aligned} \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \|\psi_{\ell}^{\mathbf{a}}(f + \Delta u_{\ell}^{\text{ex}})\|_{K} \|\delta_{\ell}^{\mathbf{a}}\|_{K} & \lesssim \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} h_{K}^{-1} \|\nabla r^{\mathbf{a},\text{hp}}\|_{K} \|\delta_{\ell}^{\mathbf{a}}\|_{K} \\ & \lesssim h_{\omega_{\ell}^{\mathbf{a}}}^{-1} \|\nabla r^{\mathbf{a},\text{hp}}\|_{\omega_{\ell}^{\mathbf{a}}} \|\delta_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} \\ & \lesssim \|\nabla r^{\mathbf{a},\text{hp}}\|_{\omega_{\ell}^{\mathbf{a}}} \|\psi_{\ell}^{\mathbf{a}}\nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}. \end{aligned} \quad (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 \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}$ is a face of an element $K \in \mathcal{T}_{\ell}^{\mathbf{a}}$, we have

$$ \|\lambda_{\ell}^{F}\|_{F} = \sup_{\substack{\mathbf{v}_{\ell} \in \text{RTN}_{\ell}^{\text{est}}(K) \\ \mathbf{v}_{\ell} \cdot \mathbf{n}_{K}|_{F} \neq 0 \\ (\mathbf{v}_{\ell}, \mathbf{r}_{\ell})_{K} = 0 \quad \forall \mathbf{r}_{\ell} \in [\mathbb{P}_{\ell}^{\text{est}}(-1)(K)]^{d}}} \frac{(\mathbf{v}_{\ell} \cdot \mathbf{n}_{K}, \lambda_{\ell}^{F})_{F}}{\|\mathbf{v}_{\ell} \cdot \mathbf{n}_{K}\|_{F}}. \quad (5.40) $$

Fix $\mathbf{v}_{\ell} \in \text{RTN}_{\ell}^{\text{est}}(K)$ with the constraints as in (5.40). Using this $\mathbf{v}_{\ell}$ as a test function in (5.34a) and the Cauchy–Schwarz inequality, we obtain

$$ \begin{aligned} (\mathbf{v}_{\ell} \cdot \mathbf{n}_{K}, \lambda_{\ell}^{F})_{F} &= (\nabla \cdot \mathbf{v}_{\ell}, \delta_{\ell}^{\mathbf{a}})_{K} - (\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}, \mathbf{v}_{\ell})_{K} \\ &\le \|\nabla \cdot \mathbf{v}_{\ell}\|_{K} \|\delta_{\ell}^{\mathbf{a}}\|_{K} + \|\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{K} \|\mathbf{v}_{\ell}\|_{K}. \end{aligned} \quad (5.41) $$

We now treat the two terms of (5.41) separately. For the first term, we start by employing the inverse inequality

$$ \|\nabla \cdot \mathbf{v}_{\ell}\|_{K} \lesssim h_{K}^{-1} \|\mathbf{v}_{\ell}\|_{K}. $$

Furthermore, by scaling arguments, under the constraints of (5.40), we have

$$ \|\mathbf{v}_{\ell}\|_{K} \lesssim h_{F}^{\frac{1}{2}} \|\mathbf{v}_{\ell} \cdot \mathbf{n}_{K}\|_{F}. \quad (5.42) $$

The first term of (5.41) can then be bounded as follows:

$$ \begin{aligned} \|\nabla \cdot \mathbf{v}_{\ell}\|_{K} \|\delta_{\ell}^{\mathbf{a}}\|_{K} & \lesssim h_{K}^{-1} \|\mathbf{v}_{\ell}\|_{K} \|\delta_{\ell}^{\mathbf{a}}\|_{K} \lesssim h_{F}^{-\frac{1}{2}} \|\mathbf{v}_{\ell} \cdot \mathbf{n}_{K}\|_{F} \|\delta_{\ell}^{\mathbf{a}}\|_{K} \\ & \lesssim h_{F}^{\frac{1}{2}} \|\mathbf{v}_{\ell} \cdot \mathbf{n}_{K}\|_{F} \|\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}, \end{aligned} \quad (5.43) $$

where we have also used the bound (5.37) and the mesh shape-regularity yielding $h_{\omega_\ell^\mathbf{a}} \approx h_K \approx h_F$. On the other hand, (5.42) leads to the following bound on the second term of (5.41)

$$ \|\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{K} \|\mathbf{v}_{\ell}\|_{K} \lesssim h_{F}^{\frac{1}{2}} \|\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{K} \|\mathbf{v}_{\ell} \cdot \mathbf{n}_{K}\|_{F}. \quad (5.44) $$

In order to bound the supremum in (5.40), we combine (5.43), (5.44), and (5.41), to get

$$ \|\lambda_{\ell}^{F}\|_{F} \lesssim h_{F}^{\frac{1}{2}} \|\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}. \quad (5.45) $$

Afterwards, $\|\psi_{\ell}^{\mathbf{a}}\|_{F} \le 1$ and Lemma 5.4 yield

$$ \|\psi_{\ell}^{\mathbf{a}}[\nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_{F}]\|_{F} \le \|[\nabla u_{\ell}^{\text{ex}}] \cdot \mathbf{n}_{F}\|_{F} \lesssim h_{F}^{-\frac{1}{2}} \|\nabla r^{\mathbf{a},\text{hp}}\|_{\omega_{F}}. \quad (5.46) $$

Using (5.45), (5.46), the shape-regularity yielding $h_F \approx h_{F'}$, $\forall F, F' \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}$, and the Cauchy–Schwarz inequality, we get

$$ \begin{aligned} \sum_{F \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}} \|\psi_{\ell}^{\mathbf{a}}[\nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_{F}]\|_{F} \|\lambda_{\ell}^{F}\|_{F} & \lesssim \sum_{F \in \mathcal{F}_{\ell}^{\mathbf{a},\text{int}}} \|\psi_{\ell}^{\mathbf{a}}[\nabla u_{\ell}^{\text{ex}} \cdot \mathbf{n}_{F}]\|_{F} h_{F}^{\frac{1}{2}} \|\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} \\ & \lesssim \|\nabla r^{\mathbf{a},\text{hp}}\|_{\omega_{\ell}^{\mathbf{a}}} \|\psi_{\ell}^{\mathbf{a}} \nabla u_{\ell}^{\text{ex}} + \sigma_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}. \end{aligned} \quad (5.47) $$

Finally, combining the bounds (5.39) and (5.47) with (5.36) proves the estimate (5.16) for $\mathbf{A} = \mathbf{I}$. $\square$
# 6 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_{\ell}^{\text{ex}}$ be the approximate solution of (4.1) on hp iteration $\ell \ge 0$. Let the pair $(\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+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_{\ell+1} = \mathbb{P}_{\mathbf{p}_{\ell+1}}(\mathcal{T}_{\ell+1}) \cap H_0^1(\Omega)$ be the finite element space to be used on iteration $(\ell+1)$ of the adaptive algorithm prescribed by Scheme 1. Let also $\underline{\eta}_{\mathcal{M}_{\ell}^{\text{ex}}}$ 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 $\eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell})$ from (4.6) is zero, in which case $u_{\ell}^{\text{ex}} = u$, and the adaptive loop terminates. Or the new numerical solution $u_{\ell+1}^{\text{ex}} \in V_{\ell+1}$ satisfies

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\| \le C_{\ell,\text{red}} \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\| \quad (6.1) $$

with

$$ 0 \le C_{\ell,\text{red}} := \sqrt{1 - \theta^2 \frac{\eta_{\mathcal{M}_{\ell}^{\text{ex}}}^2}{\eta^2(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})}} \le C_{\text{red}} := \sqrt{1 - \frac{\theta^2}{(d+1)^2 C_{\text{lb}}^2}} < 1, \quad (6.2) $$

so that $C_{\ell,\text{red}}$ is a fully computable bound on the error reduction factor and the generic constant $C_{\text{red}}$ only depends on the marking parameter $\theta$, the space dimension $d$, the mesh shape-regularity $\kappa_T$, and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient $\mathbf{A}$; here, $C_{\text{lb}}$ is the constant from (5.8) that can theoretically be chosen as (5.9).

*Proof.* Let us assume that $\eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta}) \ne 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 $(\ell+1)$, the Pythagorean relation reads

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\|^2 = \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\|^2 - \|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}})\|^2. \quad (6.3) $$

By the computable lower bound on the incremental error (4.17) and the marking criterion (4.7), we have

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}})\| \ge \|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}})\|_{\omega_{\ell}^{\text{ex}}} \ge \underline{\eta}_{\mathcal{M}_{\ell}^{\text{ex}}} = \frac{\eta_{\mathcal{M}_{\ell}^{\text{ex}}}}{\eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell})} \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}) \ge \theta \frac{\underline{\eta}_{\mathcal{M}_{\ell}^{\text{ex}}}}{\eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})} \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell}). \quad (6.4) $$

Hence, combining (6.3) with (6.4) and using the error estimate (4.6), we infer

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\|^2 \le \underbrace{\left(1 - \theta^2 \frac{\eta_{\mathcal{M}_{\ell}^{\text{ex}}}^2}{\eta^2(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta})}\right)}_{=C_{\ell,\text{red}}^2} \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\|^2. \quad (6.5) $$

*Step 2.* We show that the reduction factor $C_{\ell,\text{red}}$ is indeed bounded by the positive constant $C_{\text{red}}$ strictly smaller than one, again defined in (6.2). We decompose the elementwise estimators from (4.6) using the partition of unity $\sum_{\mathbf{a} \in V_K} \psi_{\mathbf{a}}^{\mathbf{a}}|_K = 1$, the triangle inequality, the Cauchy–Schwarz inequality, and the fact
that each simplex has (d + 1) vertices, as

$$
\begin{align}
\eta^2(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta}) &= \sum_{K \in \mathcal{M}_{\ell}^{\theta}} \left[ \|\mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}\|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \|f - \nabla \cdot \boldsymbol{\sigma}_{\ell}\|_{K} \right]^2 \nonumber \\
&= \sum_{K \in \mathcal{M}_{\ell}^{\theta}} \left[ \left\| \sum_{\mathbf{a} \in \mathcal{V}_{\ell,K}} (\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}) \right\|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \left\| \sum_{\mathbf{a} \in \mathcal{V}_{\ell,K}} (\psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{ex}} - \nabla \cdot \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}) \right\|_{K} \right]^2 \nonumber \\
&\le \sum_{K \in \mathcal{M}_{\ell}^{\theta}} \left[ \sum_{\mathbf{a} \in \mathcal{V}_{\ell,K}} \left( \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \|g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}})\|_{K} \right) \right]^2 \nonumber \\
&\le \sum_{K \in \mathcal{M}_{\ell}^{\theta}} (d+1) \sum_{\mathbf{a} \in \mathcal{V}_{\ell,K}} \left[ \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \|g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}})\|_{K} \right]^2 \nonumber \\
&\le (d+1) \sum_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \left[ \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \|g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}})\|_{K} \right]^2. \tag{6.6}
\end{align}
$$

We note that in (6.6), we have crucially employed the extended set of marked vertices $\mathcal{V}_{\ell}^{\sharp}$ of Section 4.2, as well as the divergence constraint from (4.4a) together with $\sum_{\mathbf{a} \in \mathcal{V}_{\ell,K}} \nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}^{\text{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 $\eta_{\mathcal{M}_{\ell}^{\sharp}}$ given by (4.17) yield

$$
\eta^2(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta}) \le (d+1) \sum_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} C_{\text{lb}}^2 \|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a},\text{hp}}\|_{\omega_{\ell}^{\mathbf{a}}}^2 \le (d+1)^2 C_{\text{lb}}^2 \eta_{\mathcal{M}_{\ell}^{\sharp}}^2. \quad (6.7)
$$

Finally, employing (6.7) to bound $C_{\ell,\text{red}}$ in (6.5) from above and defining the constant $C_{\text{red}}$ as in (6.2) finishes the proof. $\square$

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_{\ell \to \infty} \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}^{\text{ex}}) \| = 0.
$$

We finish this section by several remarks.

**Remark 6.3** (*p*-robust computable reduction factors $C_{\ell,\text{red}}$). We remark that the computable reduction factors $C_{\ell,\text{red}}$ from (6.2) are bounded by the constant $C_{\text{red}}$ that is independent of the polynomial degrees $p_\ell$. Moreover, in contrast to the available results, see, e.g., the reduction properties in [43, 55, 22, 19, 16, 17], the reduction factor $C_{\ell,\text{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 $\omega_{\ell}^{\mathbf{a}}$, $\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}$, 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 (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
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_\ell$. Thus, the strategy of Section 4.3.3 is not proven p-robust, and, moreover, theoretically requests $p_\ell$ 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 p-refinement 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_0^1(\omega_\ell^a)$ from (5.4) in place of $H_0^1(\omega_\ell^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_{\ell+1,K} = p_{\ell,K} + \lceil \lambda p_{\ell,K} \rceil$ for a parameter $\lambda$. The polynomial degree increase has thus theoretically be proportional to $p_{\ell,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.2 is, 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_{\ell+1,K} = p_{\ell,K} + 1$ is numerically sufficient to satisfy (5.8) and thus for a p-robust contraction factor $C_{\ell,red}$. Thus, the theoretical state-of-the-art strong p-refinement of the form $p_{\ell+1,K} = p_{\ell,K} + \lceil \lambda p_{\ell,K} \rceil$ 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 $\|\mathbf{A}^{\frac{1}{2}} \nabla r^{a,h}\|_{\omega_\ell^a}$, $\|\mathbf{A}^{\frac{1}{2}} \nabla r^{a,p}\|_{\omega_\ell^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 $h_K \|f - \nabla \cdot \boldsymbol{\sigma}_\ell\|_K / (\pi c_{A,K}^{1/2})$ can converge two orders of magnitude faster than the error and estimators, cf. [34, Remark 3.6] or [28]. Indeed, recall that $p_{\ell,K}^{\text{est}}$ is the maximum of the polynomial degrees $p_{\ell,K}$ over all elements sharing the vertex **a**. Then, from [28, Lemma 3.6],

$$ (f - \nabla \cdot \boldsymbol{\sigma}_\ell, v_h)_K = 0 \quad \forall v_h \in \mathbb{P}_{\min_{\ell \in \mathcal{V}_K}} p_{\ell,K}^{\text{est}}(K), \quad (6.8) $$

where the minimum goes over all vertices $\mathcal{V}_K$ of the element $K$. Thus, $f$ minus the divergence of $\boldsymbol{\sigma}_\ell$ is orthogonal to at least the polynomials of order $p_{\ell,K}$ in each mesh element $K \in \mathcal{T}_\ell$. Then, at least when the polynomial degrees are uniform around $K$,

$$ \frac{h_K}{\pi c_{A,K}^{1/2}} \|f - \nabla \cdot \boldsymbol{\sigma}_\ell\|_K \lesssim h_K \left( \frac{h_K}{p_{\ell,K}} \right)^{p_{\ell,K}+1} \|f\|_{H^{p_{\ell,K}+1}(K)} $$

for a regular-enough right-hand side $f$ (one additional division by $p_{\ell,K}$ is lost in (4.6) where only the division by the Poincaré constant $\pi$ is used for the purpose of obtaining a computable constant). As for $\eta_{osc,a}$ from (5.6), we similarly have, at least when $A$ is piecewise constant, $A \in \mathbb{P}_0 (T_\ell^a)^{d \times d}$,

$$ \frac{h_K}{\pi c_{A,K}^{1/2}} \|g_\ell^a - \Pi_\ell^a(g_\ell^a)\|_K = \frac{h_K}{\pi c_{A,K}^{1/2}} \|\psi_\ell^a f - \Pi_\ell^a(\psi_\ell^a f)\|_K \\ \lesssim h_K \left( \frac{h_K}{p_{\ell,K}^{\text{est}}} \right)^{p_{\ell,K}^{\text{est}}+1} \|\psi_\ell^a f\|_{H^{p_{\ell,K}^{\text{est}}+1}(K)} \lesssim \left( \frac{h_K}{p_{\ell,K}^{\text{est}}} \right)^{p_{\ell,K}^{\text{est}}+1} \|f\|_{H^{p_{\ell,K}^{\text{est}}+1}(K)} $$
which still converges one order of magnitude faster (in terms of $h_K/p_{\mathbf{a}}^{\text{est}}$) 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 \in H^{-1}(\Omega)$, the discussion is much more subtle, because in such a case, $\|f - \nabla \cdot \sigma_{\ell}\|_{K}$ or $\|g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}}(g_{\ell}^{\mathbf{a}})\|_{K}$ are not even defined. For this reason, $f \in L^{2}(\Omega)$ is typically assumed in a posteriori error analysis. For $f \in H^{-1}(\Omega)$ 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 \in H^{-1}(\Omega)$, 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^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\|$ 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.

# 7 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.

## 7.1 Adaptive sub-loop of ONE_SOLVER STEP and ESTIMATE

On the initial step $\ell = 0$, let the exact algebraic solver and hp-refinement criteria from Section 4 be used. Let henceforth $\ell \ge 1$ be the current iteration number. We employ the module ONE_SOLVER STEP exactly as in [26], with its output being an approximation $u_{\ell} := \sum_{n=1}^{N_{\ell}} (U_{\ell})_{n} \psi_{\ell}^{n} \in V_{\ell}$ of the unavailable exact Galerkin approximation $u_{\ell}^{\text{ex}}$ given by (4.1). We recall the notation in (4.2) and let the algebraic residual vector $R_{\ell}$ associated with $U_{\ell}$ be given by

$$R_{\ell} := F_{\ell} - A_{\ell}U_{\ell}. \qquad (7.1)$$

Moreover, we employ its functional representation $\mathbf{r}_{\ell} \in \mathbb{P}_{\mathbf{p}_{\ell}}(\mathcal{T}_{\ell})$, $\mathbf{r}_{\ell}|_{\partial\Omega} = 0$, satisfying

$$(\mathbf{r}_{\ell}, \psi_{\ell}^{n}) = (R_{\ell})_{n} \quad 1 \le n \le N_{\ell}. \qquad (7.2)$$

For details regarding the construction of such $\mathbf{r}_{\ell}$, 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 (7.1)

$$(\mathbf{r}_{\ell}, v_{\ell}) = (f, v_{\ell}) - (\mathbf{A}\nabla u_{\ell}, \nabla v_{\ell}) \quad v_{\ell} \in V_{\ell} \qquad (7.3)$$

that we solely use here, so that the forthcoming developments are independent of the choice of the basis $\psi_{\ell}^{n}$ in (4.2).

**Definition 7.1 (Discretization flux reconstruction $\sigma_{\ell,\text{dis}}$ by local minimization).** Let $u_{\ell}$ be the current inexact approximation. For each vertex $\mathbf{a} \in V_{\ell}$, let the local Raviart-Thomas space $V_{\ell}^{\mathbf{a}}$ be prescribed by (4.3) with, however, $p_{\ell}^{\text{est}}$ increased by one (this is related to p-robustness as discussed in [46, Remark 7.10]). Let the local discretization flux reconstruction $\sigma_{\ell,\text{dis}} \in V_{\ell}^{\mathbf{a}}$ be given by the following local minimization problem

$$\sigma_{\ell,\text{dis}}^{\mathbf{a}} := \arg \min_{\nabla \cdot \mathbf{v}_{\ell} = \Pi_{\ell}^{\mathbf{a}}(\psi_{\ell}^{\mathbf{a}} f - \nabla \psi_{\ell}^{\mathbf{a}} \cdot \mathbf{A} \nabla u_{\ell}) - \psi_{\ell}^{\mathbf{a}} \mathbf{r}_{\ell}} \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell} + \mathbf{A}^{-\frac{1}{2}} \mathbf{v}_{\ell}\|_{\omega_{\ell}^{\mathbf{a}}}. \qquad (7.4)$$
Then, the discretization flux reconstruction is defined as $\sigma_{\ell, \text{dis}} := \sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \sigma_{\ell, \text{dis}}^{\mathbf{a}}$, with each local contribution $\sigma_{\ell, \text{dis}}$ being extended by zero outside of $\omega_{\ell}^{\mathbf{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 $\sigma_{\ell, \text{alg}}$ (a piecewise in Raviart-Thomas polynomial) in $\mathbf{H}(\text{div}, \Omega)$, such that, thanks to this increase,

$$ \nabla \cdot \sigma_{\ell, \text{alg}} = \tau_{\ell}. \qquad (7.5) $$

Then, as established in [26], cf. also [33, 46], the following upper bounds on the total and algebraic error hold:

$$ \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell})\|}_{\text{total error}} \le \eta(u_{\ell}, \mathcal{T}_{\ell}) := \left\{ \sum_{K \in \mathcal{T}_{\ell}} \eta_{K}^{2}(u_{\ell}) \right\}^{\frac{1}{2}}, \quad \text{with} \quad \eta_{K}(u_{\ell}) := \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla u_{\ell} + \mathbf{A}^{-\frac{1}{2}}\sigma_{\ell, \text{dis}}\|_{K}}_{\eta_{\text{dis}, K}(u_{\ell})} + \underbrace{\|\mathbf{A}^{-\frac{1}{2}}\sigma_{\ell, \text{alg}}\|_{K}}_{\eta_{\text{alg}, K}(u_{\ell})} + \underbrace{\frac{h_{K}}{\pi c_{\mathbf{A}, K}^{1/2}} \|f - \nabla \cdot (\sigma_{\ell, \text{dis}} + \sigma_{\ell, \text{alg}})\|_{K}}_{\eta_{\text{osc}, K}(u_{\ell})} \quad (7.6a) $$

$$ \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell}^{\text{ex}} - u_{\ell})\|}_{\text{algebraic error}} \le \eta_{\text{alg}}(u_{\ell}, \mathcal{T}_{\ell}) := \left\{ \sum_{K \in \mathcal{T}_{\ell}} \eta_{\text{alg}, K}^{2}(u_{\ell}) \right\}^{\frac{1}{2}}. \quad (7.6b) $$

Next, for each vertex $\mathbf{a} \in \mathcal{V}_{\ell}$, let $V_{\ell}^{\mathbf{a}}$ be a suitable finite-dimensional subspace of $H_{*}^{1}(\omega_{\ell}^{\mathbf{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 $\rho_{\ell, \text{tot}} := \sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \psi_{\ell}^{\mathbf{a}} \rho_{\ell, \text{tot}} \in H_{0}^{1}(\Omega)$, where each vertex contribution $\rho_{\ell, \text{tot}}^{\mathbf{a}} \in V_{\ell}^{\mathbf{a}}$ solves the local primal finite element problem

$$ (\mathbf{A}\nabla\rho_{\ell, \text{tot}}^{\mathbf{a}}, \nabla v_{\ell})_{\omega_{\ell}^{\mathbf{a}}} = (f, \psi_{\ell}^{\mathbf{a}}v_{\ell})_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A}\nabla u_{\ell}, \nabla(\psi_{\ell}^{\mathbf{a}}v_{\ell}))_{\omega_{\ell}^{\mathbf{a}}} \quad \forall v_{\ell} \in V_{\ell}^{\mathbf{a}}. $$

Then, besides the upper bounds (7.6a) and (7.6b), the following lower bound on the total error is at our disposal, supposing $\rho_{\ell, \text{tot}} \neq 0$,

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell})\| \ge \frac{\sum_{\mathbf{a} \in \mathcal{V}_{\ell}} \|\mathbf{A}^{\frac{1}{2}}\nabla\rho_{\ell, \text{tot}}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}^{2}}{\|\mathbf{A}^{\frac{1}{2}}\nabla\rho_{\ell, \text{tot}}\|} =: \mu(u_{\ell}). \quad (7.7) $$

## 7.2 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_{\ell}$ is such that

$$ \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell}^{\text{ex}} - u_{\ell})\|}_{\text{algebraic error}} \le \tilde{\gamma}_{\ell} \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\|}_{\text{discretization error}}, \quad (7.8) $$

i.e. when the algebraic error is $\tilde{\gamma}_{\ell}$-times smaller than the discretization error. Recalling the Galerkin orthogonality of $u_{\ell}^{\text{ex}}$,

$$ \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell})\|}_{\text{total error}} = \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell}^{\text{ex}})\|}_{\text{discretization error}} + \underbrace{\|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell}^{\text{ex}} - u_{\ell})\|}_{\text{algebraic error}}, $$

requirement (7.8) is ensured via (7.6b) and (7.7) by the following criterion

$$ \underbrace{\eta_{\text{alg}}(u_{\ell}, \mathcal{T}_{\ell})}_{\text{algebraic error upper bound}} \le \tilde{\gamma}_{\ell} \underbrace{\sqrt{\mu^{2}(u_{\ell}) - \eta_{\text{alg}}^{2}(u_{\ell}, \mathcal{T}_{\ell})}}_{\text{discretization error lower bound}}, \quad (7.9) $$

under the condition $\eta_{\text{alg}}(u_{\ell}, \mathcal{T}_{\ell}) \le \mu(u_{\ell})$. As already used in [47, Section 6.3] and [26, Section 4.3], (7.9) is equivalent to criterion

$$ \eta_{\text{alg}}(u_{\ell}, \mathcal{T}_{\ell}) \le \gamma_{\ell} \mu(u_{\ell}) \quad (7.10) $$

with $\gamma_{\ell} = \tilde{\gamma}_{\ell}/(1+\tilde{\gamma}_{\ell}^{2})^{\frac{1}{2}} < 1$. The choice of the parameter $\gamma_{\ell}$ in (7.10) will follow from the theoretical analysis to be carried out in Section 9 below.
## 7.3 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_\ell$ 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 $\eta(u_\ell, \cdot)$ only its component corresponding to the discretization error, so the modified bulk chasing criterion to select a subset of marked vertices $V_\ell^\theta \subset V_\ell$ reads

$$ \eta_{\text{dis+osc}} \left( u_{\ell}, \bigcup_{a \in V_{\ell}^{\theta}} \mathcal{T}_{\ell}^{a} \right) \geq \theta \eta_{\text{dis+osc}}(u_{\ell}, \mathcal{T}_{\ell}), \quad (7.11) $$

where $\eta_{\text{dis+osc}}(u_\ell, S) := \{\sum_{K \in S} (\eta_{\text{dis},K}(u_\ell) + \eta_{\text{osc},K}(u_\ell))^2\}^{1/2}$, for any subset $S \subset \mathcal{T}_\ell$, and $\theta \in (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_\ell^\sharp$. In the module REFINE, the new pair $(\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+1})$ is determined on the basis of the two local problems (4.10), with an approximate solution $u_\ell$ in place of the exact Galerkin approximation $u_\ell^\text{ex}$ (recall that $u_\ell^\text{ex}$ is not at our disposal here), by one of the strategies of Sections 4.3.2 or 4.3.3.

## 7.4 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_\ell^\sharp$, let the residual lifting $r^{a, hp} \in V_\ell^{a, hp}$ be defined as the solution of the local problem (4.13), with $u_\ell^\text{ex}$ replaced by the inexact approximation $u_\ell$:

$$ (\mathbf{A}\nabla r^{a, hp}, \nabla v^{a, hp})_{\omega_\ell^a} = (f, v^{a, hp})_{\omega_\ell^a} - (\mathbf{A}\nabla u_\ell, \nabla v^{a, hp})_{\omega_\ell^a} \quad \forall v^{a, hp} \in V_\ell^{a, hp}, \quad (7.12) $$

with the local space $V_\ell^{a, hp}$ given by (4.12). Extending the residual liftings $r^{a, hp}$ by zero outside $\omega_\ell^a$, let $u_{\ell+1}^\text{ex} \in V_{\ell+1}$ be the (unavailable) exact Galerkin approximation on the next level. By [26, Lemma 5.1] used with the extended set $V_\ell^\sharp$ in place of $V_\ell^\theta$, cf. (4.17), the following lower bound holds true

$$ \| \mathbf{A}^{\frac{1}{2}} \nabla (u_{\ell+1}^{\text{ex}} - u_{\ell}) \|_{\omega_{\ell}^{\sharp}} \geq \underline{\eta}_{M_{\ell}^{\sharp}}^{*}, \quad \underline{\eta}_{M_{\ell}^{\sharp}}^{*} := \begin{cases} \frac{\sum_{a \in V_{\ell}^{\sharp}} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{a, hp} \|_{\omega_{\ell}^{a}}^{2}}{\| \mathbf{A}^{\frac{1}{2}} \nabla (\sum_{a \in V_{\ell}^{\sharp}} r^{a, hp}) \|_{\omega_{\ell}^{\sharp}}} & \text{if } \sum_{a \in V_{\ell}^{\sharp}} r^{a, hp} \neq 0, \\ 0 & \text{otherwise.} \end{cases} \quad (7.13) $$

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_{\ell+1}$ and $u_\ell$ 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 $\gamma_\ell$ 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).

## 8.1 Local stability

Recall that $\eta_{\text{osc},a}$ is defined by (5.6), with $g_\ell^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_{\ell} \in V_{\ell}$ satisfy the hat function orthogonality

$$ (f - \mathbf{r}_{\ell}, \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A}\nabla u_{\ell}, \nabla \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} = 0 \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}} $$

and let the local equilibrated flux $\sigma_{\ell,\text{dis}}^{\mathbf{a}}$ be constructed by (7.4). Let data oscillation do not dominate in that

$$ \eta_{\text{osc},\mathbf{a}} \le \frac{1}{(1+2C_{\text{st}})} \|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell} + \mathbf{A}^{-\frac{1}{2}} \sigma_{\ell,\text{dis}}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}} \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}. \quad (8.1) $$

Let the algebraic error flux reconstruction $\sigma_{\ell,\text{alg}} \in \mathbf{H}(\text{div}, \Omega)$ satisfy (7.5). Then, for a sufficiently fine patch space $V_{\ell}^{\mathbf{a},\text{hp}}$ from (4.12) and $\mathbf{r}^{\mathbf{a},\text{hp}}$ given by (7.12), there holds

$$ \left\{ \sum_{K \in \mathcal{T}_{\ell}^{\mathbf{a}}} \left[ \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{dis}}^{\mathbf{a}} \|_{K} + \frac{h_K}{\pi c_{\mathbf{A}, K}^{1/2}} \| g_{\ell}^{\mathbf{a}} - \Pi_{\ell}^{\mathbf{a}} (g_{\ell}^{\mathbf{a}}) \|_{K} \right]^2 \right\}^{\frac{1}{2}} \le C_{\text{lb}} ( \| \mathbf{A}^{\frac{1}{2}} \nabla \mathbf{r}^{\mathbf{a}, \text{hp}} \|_{\omega_{\ell}^{\mathbf{a}}} + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|_{\omega_{\ell}^{\mathbf{a}}} ) \quad (8.2) $$

for all $\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}$, where in particular $C_{\text{lb}}$ can be chosen as the maximum of (5.9) and 1.

*Proof.* The proof follows that of Proposition 5.1, upon replacing $u_{\ell}^{\text{ex}}$ by $u_{\ell}$ and $f$ by $f - \mathbf{r}_{\ell}$, and while working with the algebraic residual representation $\mathbf{r}_{\ell}$ 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_{\ell}^{\text{est}}$ by one in (7.4), $\Pi_{\ell}^{\mathbf{a}}(\psi_{\ell}^{\mathbf{a}}\mathbf{r}_{\ell}) = \psi_{\ell}^{\mathbf{a}}\mathbf{r}_{\ell}$. Then, on step 3, relying on (7.5), we can proceed as for [46, eq. (B.4)], so that, on step 4, we obtain

$$ \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{dis}}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} \le C_{\text{st}} C_{\text{cont,PF}} ( \| \mathbf{A}^{\frac{1}{2}} \nabla \mathbf{r}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|_{\omega_{\ell}^{\mathbf{a}}} ) + C_{\text{st}} \eta_{\text{osc},\mathbf{a}}, $$

from where the assertion follows by virtue of assumption (8.1). $\square$

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 \in \mathbb{P}_{\mathbf{p}_f}(\mathcal{T}_0)$ fulfill $\mathbf{p}_f \le \mathbf{p}_0 - 1$ and let the diffusion tensor $\mathbf{A}$ be piecewise constant with respect to $\mathcal{T}_0$, $\mathbf{A} \in \mathbb{P}_0(\mathcal{T}_0)^{d \times d}$. Let the approximate solution $u_{\ell} \in V_{\ell}$ satisfy the hat function orthogonality

$$ (f - \mathbf{r}_{\ell}, \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} - (\mathbf{A}\nabla u_{\ell}, \nabla \psi_{\ell}^{\mathbf{a}})_{\omega_{\ell}^{\mathbf{a}}} = 0 \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\text{int}}. $$

Let the local discretization flux reconstruction $\sigma_{\ell,\text{dis}}^{\mathbf{a}}$ be constructed by (7.4) and let the algebraic error flux reconstruction $\sigma_{\ell,\text{alg}}$, a piecewise Raviart-Thomas polynomial in $\mathbf{H}(\text{div}, \Omega)$, satisfy (7.5). Finally, let the patch space $V_{\ell}^{\mathbf{a},\text{hp}}$ from (4.12) involve the interior nodes (4.15) or the polynomial degree increase (4.16) and let $\mathbf{r}^{\mathbf{a},\text{hp}}$ be given by (7.12). Then there exists a constant $C_{\text{lb}} \ge 1$ only depending on the space dimension $d$, the mesh shape-regularity $\kappa_{\mathcal{T}}$, the ratio of the largest and the smallest eigenvalue of the diffusion coefficient $\mathbf{A}$, and the maximal polynomial degree $p_{\text{max}}$ such that

$$ \| \psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{dis}}^{\mathbf{a}} \|_{\omega_{\ell}^{\mathbf{a}}} \le C_{\text{lb}} ( \| \mathbf{A}^{\frac{1}{2}} \nabla \mathbf{r}^{\mathbf{a}, \text{hp}} \|_{\omega_{\ell}^{\mathbf{a}}} + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|_{\omega_{\ell}^{\mathbf{a}}} ) \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}. \quad (8.3) $$

*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

$$ \| \psi_{\ell}^{\mathbf{a}} (f + \Delta u_{\ell} - \mathbf{r}_{\ell}) \|_{K} \le \| \psi_{\ell}^{\mathbf{a}} \|_{\infty, K} \| f + \Delta u_{\ell} - \mathbf{r}_{\ell} \|_{K} \le \| f + \Delta u_{\ell} \|_{K} + \| \mathbf{r}_{\ell} \|_{K} \quad \forall K \in \mathcal{T}_{\ell}^{\mathbf{a}}, \quad (8.4) $$

and to treat the term $\|\mathbf{r}_{\ell}\|_{K}$, we use property (7.5) of the algebraic error flux reconstruction $\sigma_{\ell,\text{alg}}$ and an inverse inequality yielding

$$ \|\mathbf{r}_{\ell}\|_{K} \le \|\nabla \cdot \boldsymbol{\sigma}_{\ell, \text{alg}}\|_{K} \lesssim h_{K}^{-1} \|\boldsymbol{\sigma}_{\ell, \text{alg}}\|_{K}. \quad (8.5) $$

$\square$
Let

$$ C_{\text{lb}}^{\text{a}} := \max \left\{ \left\{ \sum_{K \in \mathcal{T}_{\ell}^{\text{a}}} \left[ \| \psi_{\ell}^{\text{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{dis}} \|_{K} + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \| g_{\ell}^{\text{a}} - \Pi_{\ell}^{\text{a}}(g_{\ell}^{\text{a}}) \|_{K} \right]^{2} \right\}^{\frac{1}{2}} \right\}, \quad 1 \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}, \quad (8.6) $$

$$ C_{\text{lb}}^{\ell} := \max_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} C_{\text{lb}}^{\text{a}}. $$

From (8.2) or (8.3), we see that $C_{\text{lb}}^{\ell} \le C_{\text{lb}}$.

## 8.2 Global stability

We show here how the local stability (8.2) or (8.3) implies a corresponding global stability, as in the passage from (5.8) to (6.7) in the exact solver setting.

**Lemma 8.3** (Discrete global stability of the flux equilibration, inexact solver). Let $\ell \ge 1$, let $\eta_{\text{dis+osc}}(u_\ell, \mathcal{M}_\ell^\theta)$ be given by (7.11), let $\eta_{\mathcal{M}_\ell^\sharp}^*$ 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

$$ \xi^2 := \frac{1}{2(d+1)^2 C_{\text{lb}}^2} \le \frac{1}{2(d+1)^2 (C_{\text{lb}}^{\ell})^2} \le \frac{\left(\eta_{\mathcal{M}_\ell^\sharp}^*\right)^2 + \|\mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}}\|^2}{\eta_{\text{dis+osc}}^2(u_\ell, \mathcal{M}_\ell^\theta)} =: \xi_\ell^2. \quad (8.7) $$

**Remark 8.4** (Boundedness of $\xi_\ell$ from above). Observe that unless $\eta(u_\ell, \mathcal{T}_\ell) = 0$, in which case $u_\ell = u$ and the adaptive loop terminates, (7.11) together with (9.7) below implies $\eta_{\text{dis+osc}}(u_\ell, \mathcal{M}_\ell^\theta) \ne 0$ and thus $\xi_\ell < \infty$; actually the first two inequalities in (9.9) below give $\theta^2 \xi_\ell^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 \in \mathcal{T}_\ell$,

$$ \nabla \cdot \boldsymbol{\sigma}_{\ell, \text{dis}} = \sum_{\mathbf{a} \in \mathcal{V}_{\ell, K}} \nabla \cdot \boldsymbol{\sigma}_{\ell, \text{dis}}^{\text{a}} = \sum_{\mathbf{a} \in \mathcal{V}_{\ell, K}} \left\{ \Pi_{\ell}^{\text{a}} (\psi_{\ell}^{\text{a}} f - \nabla \psi_{\ell}^{\text{a}} \cdot \mathbf{A} \nabla u_{\ell}) - \psi_{\ell}^{\text{a}} \mathbf{r}_{\ell} \right\} = \sum_{\mathbf{a} \in \mathcal{V}_{\ell, K}} \Pi_{\ell}^{\text{a}} (g_{\ell}^{\text{a}}) - \nabla \cdot \boldsymbol{\sigma}_{\ell, \text{alg}}. $$

Thus, employing the notation from (7.6a) and that below (7.11), we obtain

$$ \eta_{\text{dis+osc}}^2(u_\ell, \mathcal{M}_\ell^\theta) = \sum_{K \in \mathcal{M}_\ell^\theta} \left[ \| \mathbf{A}^{\frac{1}{2}} \nabla u_\ell + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{dis}} \|_K + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \| f - \nabla \cdot (\boldsymbol{\sigma}_{\ell, \text{dis}} + \boldsymbol{\sigma}_{\ell, \text{alg}}) \|_K \right]^2 \\ \le (d+1) \sum_{\mathbf{a} \in \mathcal{V}_\ell^\sharp} \sum_{K \in \mathcal{T}_\ell^\text{a}} \left[ \| \psi_\ell^\text{a} \mathbf{A}^{\frac{1}{2}} \nabla u_\ell + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{dis}} \|_K + \frac{h_K}{\pi c_{\mathbf{A},K}^{1/2}} \| g_\ell^\text{a} - \Pi_\ell^\text{a}(g_\ell^\text{a}) \|_K \right]^2, \quad (8.8) $$

as in (6.6). Recall now that $C_{\text{lb}}^{\ell}$ from (8.6) are bounded by $C_{\text{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 \le 2a^2+2b^2$ together with the localization of the lower bound $\eta_{\mathcal{M}_\ell^\sharp}^*$ in the spirit of (4.18) leads to

$$ \eta_{\text{dis+osc}}^2(u_\ell, \mathcal{M}_\ell^\theta) \le 2(d+1)(C_{\text{lb}}^\ell)^2 \left( \sum_{\mathbf{a} \in \mathcal{V}_\ell^\sharp} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}, hp} \|_{\omega_\ell^\mathbf{a}}^2 + \sum_{\mathbf{a} \in \mathcal{V}_\ell^\sharp} \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|_{\omega_\ell^\mathbf{a}}^2 \right) \\ \le 2(d+1)^2 (C_{\text{lb}}^\ell)^2 \left( \frac{\sum_{\mathbf{a} \in \mathcal{V}_\ell^\sharp} \| \mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a}, hp} \|_{\omega_\ell^\mathbf{a}}^2}{d+1} + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|_{\omega_\ell^\mathbf{a}}^2 \right) \quad (8.9) \\ \le 2(d+1)^2 (C_{\text{lb}}^\ell)^2 \left( \left( \eta_{\mathcal{M}_\ell^\sharp}^* \right)^2 + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|_{\omega_\ell^\mathbf{a}}^2 \right), $$

which yields the second inequality in the assertion (8.7) and concludes the proof. $\square$
9 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 $\theta \in (0, 1]$ is the Dörfler marking parameter from (7.11) and the notation $\xi_\ell$
and $\xi$ from (8.7).

**Theorem 9.1** (Guaranteed contraction, “adaptive residual lifting” strategy of Section 4.3.2, inexact solver).
Let $u$ be the weak solution of (1.1) and let the current inexact approximation $u_\ell$, $\ell \ge 0$, be obtained by the
adaptive sub-loop of Section 7.1 employing the stopping criterion (7.10) with the parameter $\gamma_\ell$ satisfying

$$
0 < \gamma_{\ell}^{2} < \min_{\ell'=\ell-1, \ell} \left\{ \frac{\theta^{2} \xi_{\ell'}^{2}}{4(1 + 3\theta^{2} \xi_{\ell'}^{2})} \right\}, \quad \ell \geq 1. \qquad (9.1)
$$

Owing to the exact initial solve, set $\gamma_\ell = 0$ and $\sigma_{\ell,\text{alg}} = 0$ for $\ell = 0$. Let the pair $(\mathcal{T}_{\ell+1}, \mathbf{p}_{\ell+1})$ be obtained by the module REFINE of Section 7.3 with the “adaptive residual lifting” strategy of Section 4.3.2 and let $\eta_{\mathcal{M}_{\ell}^{\sharp}}^*$ be the computable discrete lower bound defined by (7.13). Moreover, let $u_{\ell+1} \in V_{\ell+1}$ be the inexact finite element approximation on iteration $\ell+1$ of the adaptive loop prescribed by Scheme 2, satisfying the stopping criterion (7.10) and conditions (9.1) with the iteration counter $\ell+1$ in place of $\ell$. Let finally the data oscillation do not dominate in that (8.1) holds. Then, two options arise. Either the total a posteriori error estimate $\eta(u_\ell, \mathcal{T}_\ell)$ from (7.6a) is zero, in which case $u_\ell = u$, and the adaptive loop terminates. Or $u_{\ell+1} \in V_{\ell+1}$ satisfies

$$
\| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell+1}) \| \le C_{\ell, \text{red}} \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}) \| \quad (9.2)
$$

with

$$
0 \le C_{\ell,\text{red}} := \frac{\sqrt{1 - \frac{\left(\frac{\eta_{\mathcal{M}_{\ell}^{\sharp}}^{*}}{\eta^{2}(u_{\ell}, \mathcal{T}_{\ell})}\right)^{2} + \|\mathbf{A}^{-\frac{1}{2}}\sigma_{\ell,\text{alg}}\|^{2}}}{\sqrt{1 - \gamma_{\ell+1}^{2}}} + \gamma_{\ell}^{2}}{C_{\text{red}}^{\#} := \sqrt{1 - \frac{\theta^{2}}{8(d+1)^{2}C_{\text{lb}}^{2}}}} < 1, \quad (9.3)
$$

so that $C_{\ell,\text{red}}$ is a fully computable bound on the error reduction factor and the generic constant $C_{\text{red}}^{\#}$ only depends on the marking parameter $\theta$, the space dimension $d$, the mesh shape-regularity $\kappa_{\mathcal{T}}$, and the ratio of the largest and the smallest eigenvalue of the diffusion coefficient $\mathbf{A}$; here, $C_{\text{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

$$
\| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell+1}) \|^2 = \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell+1}) \|^2 + \| \mathbf{A}^{\frac{1}{2}} \nabla (u_{\ell+1}^{\text{ex}} - u_{\ell+1}) \|^2. \quad (9.4)
$$

We proceed in five steps.

Step 1 (bound on $\|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1}^{\text{ex}})\|$). Let $\ell \ge 1$ and consider only the non trivial case when $\eta(u_\ell, \mathcal{T}_\ell) \ne 0$. Note first that the Galerkin orthogonality of $u_{\ell+1}^{\text{ex}}$ gives

$$
\| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell+1}^{\text{ex}}) \|^2 = \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}) \|^2 - \| \mathbf{A}^{\frac{1}{2}} \nabla (u_{\ell+1}^{\text{ex}} - u_{\ell}) \|^2. \quad (9.5)
$$

Employing the discrete lower bound (7.13), adding and subtracting $\|\mathbf{A}^{-\frac{1}{2}}\boldsymbol{\sigma}_{\ell,\text{alg}}\|^2$, using the adaptive stopping criterion (7.10), and relying on the total error bounds (7.7) and (7.6a), we obtain

$$
\begin{align}
\| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell+1}^{\text{ex}}) \|^2 &\le \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}) \|^2 - \frac{\left( \eta_{\mathcal{M}_{\ell}^{\sharp}}^{*} \right)^2 + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|^2}{\eta^2(u_{\ell}, \mathcal{T}_{\ell})} \eta^2(u_{\ell}, \mathcal{T}_{\ell}) + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|^2 \nonumber \\
&\le \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}) \|^2 - \frac{\left( \eta_{\mathcal{M}_{\ell}^{\sharp}}^{*} \right)^2 + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|^2}{\eta^2(u_{\ell}, \mathcal{T}_{\ell})} \eta^2(u_{\ell}, \mathcal{T}_{\ell}) + \gamma_{\ell}^2 \mu^2(u_{\ell}) \tag{9.6} \\
&\le \underbrace{\left( 1 - \frac{\left( \eta_{\mathcal{M}_{\ell}^{\sharp}}^{*} \right)^2 + \| \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell, \text{alg}} \|^2}{\eta^2(u_{\ell}, \mathcal{T}_{\ell})} + \gamma_{\ell}^2 \right)}_{=: (C_{\ell, \text{red}}^{*})^2} \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}) \|^2. \nonumber
\end{align}
$$
As we suppose an exact solve on the initial mesh, the claim (9.6) for $\ell = 0$ immediately holds with $\gamma_\ell = 0$ and $\sigma_{\ell,\text{alg}} = 0$.

Step 2 (bound on $\|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell+1})\|$). The algebraic error estimate $\eta_{\text{alg}}(u_{\ell+1}, T_{\ell+1})$ from (7.6b), the stopping criterion (7.10), and the total error lower bound (7.7), all on step $\ell + 1$, lead to

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell+1})\| \le \eta_{\text{alg}}(u_{\ell+1}, T_{\ell+1}) \le \gamma_{\ell+1} \mu(u_{\ell+1}) \le \gamma_{\ell+1} \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1})\|. $$

Step 3 (proof of the bound (9.2)). The error reduction property (9.2) follows from (9.4) and steps 1 and 2, yielding

$$ \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1})\|^2 \le (C_{\ell,\text{red}}^*)^2 \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_\ell)\|^2 + \gamma_{\ell+1}^2 \|\mathbf{A}^{\frac{1}{2}}\nabla(u - u_{\ell+1})\|^2, $$

since $\gamma_{\ell+1} < 1$ from (9.1).

Step 4 (intermediate bound on $C_{\ell,\text{red}}^*$ from (9.6)). The triangle inequality together with the stopping criterion (7.10) and the total error lower bound (7.7) and upper bound (7.6a), yield

$$ \eta(u_\ell, T_\ell) \le \eta_{\text{alg}}(u_\ell, T_\ell) + \eta_{\text{dis+osc}}(u_\ell, T_\ell) \le \gamma_\ell \mu(u_\ell) + \eta_{\text{dis+osc}}(u_\ell, T_\ell) \le \gamma_\ell \eta(u_\ell, T_\ell) + \eta_{\text{dis+osc}}(u_\ell, T_\ell). $$

Thus, since $\gamma_\ell < 1$ from (9.1), we can bound the total error estimate $\eta(u_\ell, T_\ell)$ only by means of its component corresponding to the discretization error,

$$ \eta(u_\ell, T_\ell) \le \frac{\eta_{\text{dis+osc}}(u_\ell, T_\ell)}{(1 - \gamma_\ell)}. \qquad (9.7) $$

Employing (9.7) and the marking criterion (7.11) to bound $(C_{\ell,\text{red}}^*)^2$ from (9.6) leads to

$$ (C_{\ell,\text{red}}^*)^2 \le 1 - (1 - \gamma_\ell)^2 \frac{\left(\eta_{M_\ell^\text{alg}}^*\right)^2 + \|\mathbf{A}^{-\frac{1}{2}}\boldsymbol{\sigma}_{\ell,\text{alg}}\|^2}{\eta_{\text{dis+osc}}^2(u_\ell, T_\ell)} + \gamma_\ell^2 \le 1 - (1 - \gamma_\ell)^2 \theta^2 \frac{\left(\eta_{M_\ell^\text{alg}}^*\right)^2 + \|\mathbf{A}^{-\frac{1}{2}}\boldsymbol{\sigma}_{\ell,\text{alg}}\|^2}{\eta_{\text{dis+osc}}^2(u_\ell, M_\ell^\theta)} + \gamma_\ell^2. \qquad (9.8) $$

Now, using the notation from (8.7) and developing by the Young inequality $2\gamma_\ell \le 4\gamma_\ell^2 + 1/4$, we see

$$ (C_{\ell,\text{red}}^*)^2 \le 1 - (1 - \gamma_\ell)^2 \theta^2 \xi_\ell^2 + \gamma_\ell^2 \le 1 - \frac{3}{4} \theta^2 \xi_\ell^2 + 3\gamma_\ell^2 \theta^2 \xi_\ell^2 + \gamma_\ell^2. $$

Thus, using the second condition in (9.1) and then the bound (8.7) from Lemma 8.3, we obtain

$$ 0 \le C_{\ell,\text{red}}^* < \sqrt{1 - \frac{1}{2}\theta^2\xi_\ell^2} \le C_{\text{red}}^* := \sqrt{1 - \frac{\theta^2}{4(d+1)^2 C_{\text{lb}}^2}} < 1. \qquad (9.9) $$

Step 5 (bound (9.3) on $C_{\ell,\text{red}}$). The first condition in (9.1) gives

$$ 0 < \gamma_{\ell+1}^2 \le \frac{\theta^2 \xi_\ell^2}{4(1 + 3\theta^2 \xi_\ell^2)} \le \frac{\theta^2 \xi_\ell^2}{4 - \theta^2 \xi_\ell^2}. $$

Recalling that (9.9) gives $\theta^2 \xi_\ell^2 < 2 < 4$, we then have

$$ \gamma_{\ell+1}^2 \le \frac{\theta^2 \xi_\ell^2}{4 - \theta^2 \xi_\ell^2} \iff \frac{1}{1 - \gamma_{\ell+1}^2} \le \frac{1 - \frac{\theta^2 \xi_\ell^2}{4}}{1 - \frac{\theta^2 \xi_\ell^2}{2}}. $$

Thus, from (9.9), we infer

$$ \frac{(C_{\ell,\text{red}}^*)^2}{1 - \gamma_{\ell+1}^2} \le \left(1 - \frac{\theta^2 \xi_\ell^2}{2}\right) \frac{1 - \frac{\theta^2 \xi_\ell^2}{4}}{1 - \frac{\theta^2 \xi_\ell^2}{2}} = 1 - \frac{\theta^2 \xi_\ell^2}{4} \le 1 - \frac{\theta^2}{8(d+1)^2 C_{\text{lb}}^2} = (C_{\text{red}}^*)^2 < 1, $$

where we have also employed the bound (8.7) from Lemma 8.3. $\square$
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_{\ell \to \infty} \| \mathbf{A}^{\frac{1}{2}} \nabla (u - u_{\ell}) \| = 0. $$

**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 $\gamma_{\ell} = \gamma_{\ell+1} = 0$ and $\sigma_{\ell,\text{alg}} = 0$, they give the same contraction factor $C_{\ell,\text{red}}$ when the Dörfler criterion (4.7) is employed with equality, i.e., $\eta(u_{\ell}^{\text{ex}}, \mathcal{M}_{\ell}^{\theta}) = \theta \eta(u_{\ell}^{\text{ex}}, \mathcal{T}_{\ell})$.

**Remark 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.3 when the polynomial degrees are bounded by $p_{\text{max}}$, with the non p-robust constant $C_{\text{lb}}$ from (8.3) in place of the p-robust constant $C_{\text{lb}}$ from (8.2).

**Remark 9.5** (Stopping parameter $\gamma_{\ell}$). Our analysis shows that choosing the stopping parameter $\gamma_{\ell}$ following (9.1) (practically, since both $\theta$ and $\xi_{\ell}$ 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 $\gamma_{\ell}$, leading in particular to the reasonable choice of $\gamma_{\ell} \approx 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

$$ 0 < \frac{\theta^2 \xi^2}{4(1 + 3\theta^2 \xi^2)} \le \frac{\theta^2 \xi_{\ell}^2}{4(1 + 3\theta^2 \xi_{\ell}^2)} \quad \forall \ell \ge 1. \qquad (9.10) $$

This would give a valid (theoretical, since $\xi$ from (8.7) is generally unknown) choice for $\gamma_{\ell}$, but we rather stress that (9.10) shows that picking, e.g., $\gamma_{\ell}$ as a fixed fraction of the min in (9.1) is uniformly bounded from below in all the iteration counter $\ell$, the mesh sizes, and the polynomial degrees.

**Remark 9.6** (Contraction factors $C_{\text{red}}$, $C_{\text{red}}^{*}$, and $C_{\text{red}}^{\#}$). We remark that the contraction factors $C_{\text{red}}$ from (6.2), $C_{\text{red}}^{*}$ from (9.9), and $C_{\text{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.

# 10 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 $\Omega = (-1, 1) \times (-1, 1) \setminus [0, 1] \times [-1, 0]$ with $f = 0$, $\mathbf{A} = \mathbf{I}$, and the weak solution, in polar coordinates,

$$ u(r, \varphi) = r^{\frac{2}{3}} \sin \left( \frac{2\varphi}{3} \right). $$

We start on a coarse criss-cross grid $\mathcal{T}_0$ with $\max_{K \in \mathcal{T}_0} 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 $\|\nabla(u - u_{\ell}^{\text{ex}})\|/\|\nabla 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 $\|\nabla(u - u_{\ell}^{\text{ex}})\| \le C_1 \exp(-C_2 \text{DoF}_{\ell}^{\frac{1}{3}})$ 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 $\mathcal{T}_{25}$
<figure><img src="image_8.png" /><figcaption>**Figure 4:** Relative energy error $||\nabla(u - u_{\ell}^{\text{ex}})||/||\nabla u||$ as a function of DoF$\frac{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 $\mathcal{T}_{25}$ 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.</figcaption></figure>

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_{\text{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_{\ell,\text{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

$$I_{\text{red}}^{\text{eff}} := \frac{C_{\ell,\text{red}}}{\frac{\|\nabla(u-u_{\ell+1}^{\text{ex}})\|}{\|\nabla(u-u_{\ell}^{\text{ex}})\|}} \quad (10.1)$$

We also verify the sharpness of the underlying discrete lower bound $\eta_{\mathcal{M}_{\ell}^{\sharp}}$ given by (4.17) in terms of the effectivity index defined as

$$I_{\text{LB}}^{\text{eff}} := \frac{\|\nabla(u_{\ell+1}^{\text{ex}} - u_{\ell}^{\text{ex}})\|_{\omega_{\ell}^{\sharp}}}{\eta_{\mathcal{M}_{\ell}^{\sharp}}} \quad (10.2)$$

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 $\mathbf{A} = \mathbf{I}$ here, so that no data oscillations arise from $f$ and $\mathbf{A}$), let us define

$$C_{\text{lb}}^{\mathbf{a}} := \frac{\|\psi_{\ell}^{\mathbf{a}} \mathbf{A}^{\frac{1}{2}} \nabla u_{\ell}^{\text{ex}} + \mathbf{A}^{-\frac{1}{2}} \boldsymbol{\sigma}_{\ell}^{\mathbf{a}}\|_{\omega_{\ell}^{\mathbf{a}}}}{\|\mathbf{A}^{\frac{1}{2}} \nabla r^{\mathbf{a},hp}\|_{\omega_{\ell}^{\mathbf{a}}}} \quad \forall \mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}, \quad \bar{C}_{\text{lb}}^{\ell} := \max_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} C_{\text{lb}}^{\mathbf{a}}, \quad \text{and} \quad \underline{C}_{\text{lb}}^{\ell} := \min_{\mathbf{a} \in \mathcal{V}_{\ell}^{\sharp}} C_{\text{lb}}^{\mathbf{a}} \quad (10.3)$$

after one newest-vertex bisection step or polynomial increase to $p+1$ in each marked patch $\omega_{\ell}^{\mathbf{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
<figure><img src="image_9.png" /><figcaption>**Figure 5:** Effectivity index for the error reduction factor estimate $C_{\ell,\text{red}}$ of Theorem 6.1 given by (10.1) (left); effectivity index for the discrete lower bound $\eta_{\underline{M}_{\ell}^{\text{ad}}}$ 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.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6:** Maximal and minimal values over all marked vertices of the local stability constants $C_{\text{lb}}^{\text{a}}$ given by $\bar{C}_{\text{lb}}^{\ell}$ and $\underline{C}_{\text{lb}}^{\ell}$ from (10.3). The “adaptive residual lifting” hp strategy of Section 4.3.2 corresponding to Figure 4 with maximal polynomial degree 8 (left), maximal polynomial degree 15 (right). Exact algebraic solver.</figcaption></figure>

both $\bar{C}_{\text{lb}}^{\ell}$ and $\underline{C}_{\text{lb}}^{\ell}$ in Figure 6, left (recall from Figure 4 that the maximal polynomial degree produced by the strategy on the last mesh $\mathcal{T}_{25}$ equals 8). We see that $\bar{C}_{\text{lb}}^{\ell} \le 2.2 < 10$, which explains why no iteration on h- or p-refinements appears above for our choice $C_{\text{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 $\bar{C}_{\text{lb}}^{\ell}$ and $\underline{C}_{\text{lb}}^{\ell}$ essentially do not change with respect to the left part of Figure 6, with $\bar{C}_{\text{lb}}^{\ell} \le 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_{\text{lb}} = 10$.
To support this numerical observation further, we plot in Figure 7, left, $\bar{C}_{lb}^\ell$ and $\underline{C}_{lb}^\ell$ 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_{\ell,\text{red}}$ and of $\eta_{M_\ell^2}$, and the values of $\bar{C}_{lb}^\ell$ as well as of $\underline{C}_{lb}^\ell$. We thus only focus on the stopping parameter $\gamma_\ell$ 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 $\gamma_\ell$ 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, $\gamma_\ell$ 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).

<figure><img src="image_11.png" /><figcaption>**Figure 7:** Maximal and minimal values over all marked vertices of the local stability constants $C_{lb}^a$ given by $\bar{C}_{lb}^\ell$ and $\underline{C}_{lb}^\ell$ 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 $\gamma_\ell$, the “interior node, $p+d-1$” hp strategy of Section 7 compared to the parameter $\gamma_\ell$ as in [26, Theorem 5.4]. Inexact algebraic solver (**right**).</figcaption></figure>

# 11 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 $\ell$ 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 $\ell$ 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 $\gamma_\ell$ 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 $\gamma_\ell$ 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.

References

[1] M. AINSWORTH AND J. T. ODEN, *A posteriori error estimation in finite element analysis*, Pure and Applied Mathematics (New York), Wiley-Interscience [John Wiley & Sons], New York, 2000.

[2] M. ARIOLI, E. H. GEORGOULIS, AND D. LOGHIN, *Stopping criteria for adaptive finite element solvers*, SIAM J. Sci. Comput., 35 (2013), pp. A1537-A1559.

[3] I. BABUŠKA AND A. MILLER, *A feedback finite element method with a posteriori error estimation. I. The finite element method and some basic properties of the a posteriori error estimator*, Comput. Methods Appl. Mech. Engrg., 61 (1987), pp. 1-40.

[4] I. BABUŠKA AND B. GUO, *The h-p version of finite element method, part 1: The basic approximation results*, Comp. Mech., (1986), pp. 21-41.

[5] ——, *The h-p version of finite element method, part 2: General results and application*, Comp. Mech., (1986), pp. 203-220.

[6] R. E. BANK, A. PARSANIA, AND S. SAUTER, *Saturation estimates for hp-finite element methods*, Comput. Vis. Sci., 16 (2013), pp. 195-217.

[7] R. BECKER, C. JOHNSON, AND R. RANNACHER, *Adaptive error control for multigrid finite element methods*, Computing, 55 (1995), pp. 271-288.

[8] R. BECKER AND S. MAO, *Convergence and quasi-optimal complexity of a simple adaptive finite element method*, M2AN Math. Model. Numer. Anal., 43 (2009), pp. 1203-1219.

[9] P. BINEV, *Tree approximation for hp-adaptivity*, SIAM J. Numer. Anal., 56 (2018), pp. 3346-3357.

[10] P. BINEV, W. DAHMEN, AND R. DeVORE, *Adaptive finite element methods with convergence rates*, Numer. Math., 97 (2004), pp. 219-268.

[11] J. BLECHTA, J. MÁLEK, AND M. VOHRALÍK, *Localization of the $W^{-1,q}$ norm for local a posteriori efficiency*, IMA J. Numer. Anal., 40 (2020), pp. 914-950.

[12] D. BRAESS, V. PILLWEIN, AND J. SCHÖBERL, *Equilibrated residual error estimates are p-robust*, Comput. Methods Appl. Mech. Engrg., 198 (2009), pp. 1189-1197.

[13] D. BRAESS AND J. SCHÖBERL, *Equilibrated residual error estimator for edge elements*, Math. Comp., 77 (2008), pp. 651-672.

[14] F. BREZZI AND M. FORTIN, *Mixed and hybrid finite element methods*, vol. 15 of Springer Series in Computational Mathematics, Springer-Verlag, New York, 1991.

[15] M. BÜRG AND W. DÖRFLER, *Convergence of an adaptive hp finite element strategy in higher space-dimensions*, Appl. Numer. Math., 61 (2011), pp. 1132-1146.

[16] C. CANUTO, R. H. NOCHETTO, R. P. STEVENSON, AND M. VERANI, *Convergence and optimality of hp-AFEM*, Numer. Math., 135 (2017), pp. 1073-1119.
[17] ——, *On p-robust saturation for hp-AFEM*, Comput. Math. Appl., 73 (2017), pp. 2004–2022.

[18] ——, *A saturation property for the spectral-Galerkin approximation of a Dirichlet problem in a square*, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 987–1003.

[19] C. CARSTENSEN, M. FEISCHL, M. PAGE, AND D. PRAETORIUS, *Axioms of adaptivity*, Comput. Math. Appl., 67 (2014), pp. 1195–1253.

[20] C. CARSTENSEN AND S. A. FUNKEN, *Fully reliable localized error control in the FEM*, SIAM J. Sci. Comput., 21 (1999/00), pp. 1465–1484.

[21] J. M. CASCÓN, C. KREUZER, R. H. NOCHETTO, AND K. G. SIEBERT, *Quasi-optimal convergence rate for an adaptive finite element method*, SIAM J. Numer. Anal., 46 (2008), pp. 2524–2550.

[22] J. M. CASCÓN AND R. H. NOCHETTO, *Quasioptimal cardinality of AFEM driven by nonresidual estimators*, IMA J. Numer. Anal., 32 (2012), pp. 1–29.

[23] P. CIARLET, JR. AND M. VOHRALÍK, *Localization of global norms and robust a posteriori error control for transmission problems with sign-changing coefficients*, ESAIM Math. Model. Numer. Anal., 52 (2018), pp. 2037–2064.

[24] P. DANIEL, *Adaptive hp-finite elements with guaranteed error contraction and inexact multilevel solvers*, Ph.D. thesis, Sorbonne University, March 2019. https://hal.inria.fr/tel-02104982.

[25] P. DANIEL, A. ERN, I. SMEARS, AND M. VOHRALÍK, *An adaptive hp-refinement strategy with computable guaranteed bound on the error reduction factor*, Comput. Math. Appl., 76 (2018), pp. 967–983.

[26] P. DANIEL, A. ERN, AND M. VOHRALÍK, *An adaptive hp-refinement strategy with inexact solvers and computable guaranteed bound on the error reduction factor*, Comput. Methods Appl. Mech. Engrg., 359 (2020), p. 112607.

[27] P. DESTUYNDER AND B. MÉTIVET, *Explicit error bounds in a conforming finite element method*, Math. Comp., 68 (1999), pp. 1379–1396.

[28] V. DOLEJŠÍ, A. ERN, AND M. VOHRALÍK, *hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems*, SIAM J. Sci. Comput., 38 (2016), pp. A3220–A3246.

[29] W. DÖRFLER, *A convergent adaptive algorithm for Poisson's equation*, SIAM J. Numer. Anal., 33 (1996), pp. 1106–1124.

[30] W. DÖRFLER AND V. HEUVELINE, *Convergence of an adaptive hp finite element strategy in one space dimension*, Appl. Numer. Math., 57 (2007), pp. 1108–1124.

[31] C. ERATH, G. GANTNER, AND D. PRAETORIUS, *Optimal convergence behavior of adaptive FEM driven by simple $(h-h/2)$-type error estimators*, Comput. Math. Appl., 79 (2020), pp. 623–642.

[32] A. ERN, I. SMEARS, AND M. VOHRALÍK, *Discrete p-robust $H$-(div)-liftings and a posteriori estimates for elliptic problems with $H^{-1}$ source terms*, Calcolo, 54 (2017), pp. 1009–1025.

[33] A. ERN AND M. VOHRALÍK, *Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs*, SIAM J. Sci. Comput., 35 (2013), pp. A1761–A1791.

[34] ——, *Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations*, SIAM J. Numer. Anal., 53 (2015), pp. 1058–1081.

[35] A. ERN AND M. VOHRALÍK, *Stable broken $H^1$ and $H$-(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions*, Math. Comp., 89 (2020), pp. 551–594.

[36] G. GANTNER, A. HABERL, D. PRAETORIUS, AND B. STIFTNER, *Rate optimal adaptive FEM with inexact solver for nonlinear operators*, IMA J. Numer. Anal., 38 (2018), pp. 1797–1831.
[37] M. HOLST, R. SZYPOWSKI, AND Y. ZHU, *Adaptive finite element methods with inexact solvers for the nonlinear Poisson-Boltzmann equation*, in Domain decomposition methods in science and engineering XX, vol. 91 of Lect. Notes Comput. Sci. Eng., Springer, Heidelberg, 2013, pp. 167–174.

[38] P. JIRÁNEK, Z. STRAKOŠ, AND M. VOHRALÍK, *A posteriori error estimates including algebraic error and stopping criteria for iterative solvers*, SIAM J. Sci. Comput., 32 (2010), pp. 1567–1590.

[39] C. KREUZER AND K. G. SIEBERT, *Decay rates of adaptive finite elements with Dörfler marking*, Numer. Math., 117 (2011), pp. 679–716.

[40] C. KREUZER AND A. VEESEL, *Oscillation in a posteriori error estimation*, Numer. Math., 148 (2021), pp. 43–78.

[41] W. F. MITCHELL, *Adaptive refinement for arbitrary finite-element spaces with hierarchical bases*, J. Comput. Appl. Math., 36 (1991), pp. 65–78.

[42] W. F. MITCHELL AND M. A. MCCLAIN, *A comparison of hp-adaptive strategies for elliptic partial differential equations (long version)*, NISTIR 7824, National Institute of Standards and Technology, (2011).

[43] P. MORIN, R. H. NOCHETTO, AND K. G. SIEBERT, *Convergence of adaptive finite element methods*, SIAM Rev., 44 (2002), pp. 631–658. Revised reprint of “Data oscillation and convergence of adaptive FEM” [SIAM J. Numer. Anal. **38** (2000), no. 2, 466–488; MR1770058 (2001g:65157)].

[44] ——, *Local problems on stars: a posteriori error estimators, convergence, and performance*, Math. Comp., 72 (2003), pp. 1067–1097.

[45] P. MORIN, K. G. SIEBERT, AND A. VEESEL, *A basic convergence result for conforming adaptive finite elements*, Math. Models Methods Appl. Sci., 18 (2008), pp. 707–737.

[46] J. PAPEŽ, U. RÜDE, M. VOHRALÍK, AND B. WOHLMUTH, *Sharp algebraic and total a posteriori error bounds for h and p finite elements via a multilevel approach. Recovering mass balance in any situation*, Comput. Methods Appl. Mech. Engrg., 371 (2020), p. 113243.

[47] J. PAPEŽ, Z. STRAKOŠ, AND M. VOHRALÍK, *Estimating and localizing the algebraic and total numerical errors using flux reconstructions*, Numer. Math., 138 (2018), pp. 681–721.

[48] C.-M. PFEILER AND D. PRAETORIUS, *Dörfler marking with minimal cardinality is a linear complexity problem*, Math. Comp., 89 (2020), pp. 2735–2752.

[49] A. QUARTERONI AND A. VALLI, *Numerical approximation of partial differential equations*, vol. 23 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1994.

[50] V. REY, C. REY, AND P. GOSSELET, *A strict error bound with separated contributions of the discretization and of the iterative solver in non-overlapping domain decomposition methods*, Comput. Methods Appl. Mech. Engrg., 270 (2014), pp. 293–303.

[51] J. E. ROBERTS AND J.-M. THOMAS, *Mixed and hybrid methods*, in Handbook of Numerical Analysis, Vol. II, North-Holland, Amsterdam, 1991, pp. 523–639.

[52] D. SCHÖTZAU AND C. SCHWAB, *Exponential convergence of hp-FEM for elliptic problems in polyhedra: mixed boundary conditions and anisotropic polynomial degrees*, Found. Comput. Math., 18 (2018), pp. 595–660.

[53] E. G. SEWELL, *Automatic generation of triangulations for piecewise polynomial approximation*, ProQuest LLC, Ann Arbor, MI, 1972. Thesis (Ph.D.)–Purdue University.

[54] R. STEVENSON, *An optimal adaptive finite element method*, SIAM J. Numer. Anal., 42 (2005), pp. 2188–2217.
[55] ——, *Optimality of a standard adaptive finite element method*, Found. Comput. Math., 7 (2007), pp. 245–269.

[56] R. VERFÜRTH, *A posteriori error estimation techniques for finite element methods*, Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, 2013.

[57] M. VOHRALÍK, *Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods*, Math. Comp., 79 (2010), pp. 2001–2032.