ADAPTIVE FEM WITH QUASI-OPTIMAL OVERALL COST FOR NONSYMMETRIC LINEAR ELLIPTIC PDES

. We consider a general nonsymmetric second-order linear elliptic PDE in the framework of the Lax–Milgram lemma. We formulate and analyze an adaptive finite element algorithm with arbitrary polynomial degree that steers the adaptive mesh refinement and the inexact iterative solution of the arising linear systems. More precisely, the iterative solver employs, as an outer loop, the so-called Zarantonello iteration to sym-metrize the system and, as an inner loop, a uniformly contractive algebraic solver, e.g., an optimally preconditioned conjugate gradient method or an optimal geometric multi-grid algorithm. We prove that the proposed inexact adaptive iteratively symmetrized finite element method (AISFEM) leads to full linear convergence and, for sufficiently small adaptivity parameters, to optimal convergence rates with respect to the overall computational cost, i.e., the total computational time. Numerical experiments underline the theory.


Introduction
The mathematical understanding of optimal adaptivity for finite element methods (AFEMs) has reached a high level of maturity; see, e.g., [BDD04; Ste07; CKNS08; KS11; CN12; FFP14; CFPP14] for some contributions to linear PDEs.While the focus is usually on optimal convergence rates with respect to the degrees of freedom [BDD04; CKNS08; KS11; CN12; FFP14; CFPP14], the cumulative nature of adaptivity should rather ask for optimal convergence rates with respect to the overall computational cost, i.e., the overall elapsed computational time.This, usually called optimal complexity, has been thoroughly analyzed for adaptive wavelet methods [CDD01;CDD03] and it has also been addressed in the seminal work [Ste07] on AFEM for the Poisson model problem.Recent works [GHPS21; HPW21; HPSV21] considered optimal complexity for energy minimization problems and, in particular, for symmetric linear elliptic PDEs.In contrast to this, optimal complexity for nonsymmetric linear elliptic PDEs remained an open question due to the lack of a contractive algebraic solver that is compatible with the variational structure of the PDE.Closing this gap is the topic of the present work.While the canonical candidate for solving the nonsymmetric discrete systems would be GMRES, we take a different path that is motivated by up-to-date proofs of the Lax-Milgram lemma and closely related to the Richardson iteration used in the context of optimal adaptive wavelet methods.Some comments on the challenges presented by GMRES and related future work are given below.
As a model problem, we consider the nonsymmetric second-order linear elliptic PDE Find u ⋆ ∈ X := H 1 0 (Ω) such that b(u ⋆ , v) = F (v) for all v ∈ X . (1.2) To ensure the existence and uniqueness of u ⋆ ∈ H 1 0 (Ω), we assume that the bilinear form b(•, •) is continuous and elliptic on H 1 0 (Ω) so that the Lax-Milgram lemma applies.To discretize (1.2), we employ a conforming finite element method based on a conforming simplicial triangulation T ℓ of Ω and a fixed polynomial degree m ∈ N.With the finite element formulation reads: Existence and uniqueness of u ⋆ ℓ follow again from the Lax-Milgram lemma.Note that (1.4) leads to a nonsymmetric, yet positive definite linear system of equations.To derive an optimal nonsymmetric algebraic solver, we follow the constructive proof of the Lax-Milgram lemma and reduce the discrete formulations (1.4) to symmetric problems by employing the so-called Zarantonello symmetrization (sometimes referred to as Banach-Picard fixedpoint iteration).To this end, we define the bilinear form associated with the principal part of the PDE by a(u, v) := ⟨A∇u , ∇v⟩ Ω for all u, v ∈ X . (1.5) Note that a(•, •) is continuous and elliptic on X and consult Section 2 for details.For a given damping parameter δ > 0, define the Zarantonello mapping Φ ℓ (δ; •) : (1.6) see [Zar60] or [Zei90,Section 25.4].The Riesz-Fischer theorem (and also the Lax-Milgram lemma) proves existence and uniqueness of Φ ℓ (δ; u ℓ ) ∈ X ℓ , i.e., the Zarantonello operator is well-defined.In particular, u ⋆ ℓ = Φ(δ; u ⋆ ℓ ) is the only fixpoint of Φ(δ; •) for any δ > 0.Moreover, choosing δ suitably small will lead to a contractive method to approximate u ⋆ ℓ in the spirit of the Banach fixpoint theorem with respect to the a(•, •)-induced energy norm ⦀v⦀ := a(v, v) 1/2 .At this point, it thus remains to treat a symmetric, positive definite (SPD) linear system of equations corresponding to (1.6), that can be solved iteratively in practice for instance by the use of either a conjugate gradient (CG) method with an optimal preconditioner, see e.g., [CNX12], or an optimal geometric multigrid (MG) solver, see e.g., [WZ17;IMPS22].
The proposed adaptive strategy of this work, hereafter referred to as AISFEM, begins with the initial guess u 0,0 0 := u 0,j 0 := u 0,⋆ 0 := 0 ∈ X 0 associated to a coarse mesh T 0 .Finite element approximations u k,j ℓ ∈ X ℓ are successively computed, where ℓ ∈ N 0 is the meshrefinement index of the ℓ-th adaptively refined mesh.More precisely, u k,j ℓ is obtained after j algebraic solver steps in the k-th step of the Zarantonello symmetrization approximating the unique u k,⋆ ℓ := Φ ℓ (δ; u k−1,j ℓ ) ∈ X ℓ , where u k−1,j ℓ ∈ X ℓ denotes the final approximation of u k−1,⋆ ℓ when the algebraic solver is adaptively terminated.In particular, our analysis provides stopping criteria for the algebraic solver as well as the (perturbed) Zarantonello symmetrization.We give a schematic view of our approach in Figure 1; see Algorithm A in Section 3 below for the formal statement.Overall, the adaptive strategy thus leads to a triple index set ℓ is used by the AISFEM Algorithm A}, ( equipped with the natural lexicographic order |•, •, •|.This enables us to present the main contributions of this work: First, in the spirit of [GHPS21; HPSV21], we prove that the quasi-error which is the sum of the overall error plus the algebraic solver error plus the residual error estimator, is linearly convergent with respect to the order of is computed earlier than u k,j ℓ within the (sequential) adaptive loop and |ℓ, k, j| − |ℓ ′ , k ′ , j ′ | ∈ N 0 is the overall number of discretization, symmetrization, and algebraic solver steps in between.In explicit terms, Theorem 4.1 proves the existence of constants C lin > 0 and 0 < q lin < 1 as well as an index ℓ 0 ∈ N 0 such that, for all (1.9) The threshold level ℓ 0 ∈ N 0 arises from the lack of Galerkin orthogonality with respect to the a(•, •)-induced energy norm leading to a more involved analysis.Second, as shown in Corollary 4.2, this implies that, for any s > 0, there holds the equivalence (1.10) The interpretation of (1.10) is that the AISFEM algorithm leads to algebraic convergence rate s > 0 with respect to the degrees of freedom (finite left-hand side) if and only if it leads to algebraic convergence rate s with respect to the overall computational cost (finite right-hand side), i.e., with respect to the computational time.Third, extending available results from the literature [CN12; FFP14; BHP17], Theorem 4.3 proves that, for sufficiently small adaptivity parameters, the proposed algorithm has optimal complexity (which follows from optimal rates with respect to the degrees of freedom and (1.10)).Finally, we admit that the proposed strategy hinges crucially on the appropriate (sufficiently small) choice of the Zarantonello parameter δ > 0 in (1.6) as well as on the parameter λ alg > 0 in the stopping criterion for the algebraic solver in Algorithm A(i.b.II) below.If these parameters are chosen too large, the proposed method may fail to converge.Besides this restriction, linear convergence (1.9) is guaranteed for any choice of the other adaptivity parameters λ sym , θ, C mark (see Algorithm A below).
Outline.The remainder of the work is organized as follows.Section 2 focuses on the setting and underlying assumptions.In Section 3, we present the AISFEM algorithm in full detail and highlight some of its properties.The main results of this work are presented in Section 4, the proofs of which are given in Section 5. Numerical experiments in Section 6 underline the theoretical results, before the short Section 7 concludes our results and outlines future work.Throughout, A ≲ B denotes A ≤ c B with a generic constant c > 0 that is independent of the discretization, but may depend on all problem parameters.Moreover, A ≃ B abbreviates A ≲ B ≲ A.

Preliminaries
In this section, we state all prerequisites to formulate the AISFEM algorithm (Algorithm A in Section 3 below).In particular, we collect the contraction properties of the Zarantonello symmetrization, the algebraic solver, the mesh-refinement strategy, and the required properties of the a posteriori error estimator.
Remark 2.1.The contractive Zarantonello symmetrization and hence the results of this work hold in an abstract framework beyond that of the introduction in Section 1.More precisely, the analysis allows for an abstract separable Hilbert space X over K ∈ {R, C} with norm ‖•‖ X and a weak formulation (2.1), where a(•, •) is a Hermitian and continuous sesquilinear form on X and K : X → X ′ is a compact linear operator such that b(•, •) is elliptic and continuous on X .Provided that a contractive algebraic solver is used (see Section 2.5), the analysis thus also applies to other boundary conditions (e.g., mixed Dirichlet-Neumann-Robin instead of homogeneous Dirichlet boundary conditions used in the introduction).

Mesh refinement.
From now on, let T 0 be a given conforming triangulation of Ω ⊂ R d with d ≥ 1 which is admissible in the sense of [Ste08] for d ≥ 3.For mesh refinement, we employ newest vertex bisection (NVB); see [AFF + 15] for d = 1, [Ste08] for d ≥ 2 and [KPP13] for d = 2 with non-admissible T 0 .For each triangulation T H and marked elements M H ⊆ T H , let T h := refine(T H , M H ) be the coarsest conforming triangulation where all T ∈ M H have been refined, i.e., M H ⊆ T H \T h .We write T h ∈ T(T H ) if T h results from T H by finitely many steps of refinement and, for N ∈ N 0 , we write To abbreviate notation, let T := T(T 0 ).Throughout, each triangulation T H ∈ T is associated with a finitedimensional finite element space X H ⊂ X , see (1.3), and refinement Within the setting of AFEM, we will work with a hierarchy {T ℓ } ℓ∈N 0 generated by NVB refinements from the initial mesh T 0 .

2.3.
A posteriori error estimator and axioms of adaptivity.For T H ∈ T, let be the local contributions of some computable error estimator.We define To abbreviate notation, let η H (v H ) := η H (T H ; v H ). Furthermore, we suppose that η H satisfies the following axioms of adaptivity from [CFPP14] with constants C stab , C rel , C drel > 0 and 0 < q red < 1 only depending on the dimension d, the polynomial degree m, and shape regularity of T 0 : (A1) stability: For all T H ∈ T and T h ∈ T(T H ), all v h ∈ X h and all v H ∈ X H , and every U H ⊆ T H ∩ T h , it holds that (A2) reduction: For all T H ∈ T and T h ∈ T(T H ), and all v H ∈ X H , it holds that (A3) reliability: For all T H ∈ T, the exact solutions u ⋆ ∈ X of (1.2) and u ⋆ H ∈ X H of (1.4) satisfy that (A4) discrete reliability: For all T H ∈ T and T h ∈ T(T H ), the corresponding exact discrete solutions satisfy that ).We note that these axioms (A1)-(A4) are satisfied for the standard residual error estimators; see Section 6 below for the model problem (1.1) from the introduction.
2.5.Contractive algebraic solver.We assume that we have at hand an iterative algebraic solver with iteration step Ψ H : X ′ × X H → X H .This means, given a linear and continuous functional G ∈ X ′ and an approximation w H ∈ X H of the unique solution the algebraic solver returns an improved Ψ H (G; w H ) ∈ X H in the sense that there exists a constant 0 < q alg < 1, which is independent of G and X H , such that (2.9) To simplify notation when the right-hand side G is complicated or lengthy (as for the Zarantonello iteration (1.6)), we shall write Ψ H (w ⋆ H ; •) instead of Ψ H (G; •), even though w ⋆ H is unknown and will never be computed.
In the framework of AFEM, possible examples for such contractive solvers include optimally preconditioned conjugate gradient methods or optimal geometric multigrid methods, see, e.g., [CNX12] or [WZ17], respectively, for approaches focused on lowest-order discretizations and [IMPS22] for an optimal multigrid method which is also robust with respect to the polynomial degree.

Completely adaptive algorithm
In the following, we formulate an inexact adaptive iteratively symmetrized finite element method (AISFEM) in the spirit of [HPSV21].For ease of presentation, we make the following conventions: Algorithm A defines certain terminal indices ℓ, k[ℓ], j[ℓ, k], indicated by underlining.We shall omit the arguments of k and j if these are clear from the context, e.g., we simply write A similar convention will be used for triple indices, e.g., (ℓ, k, j) = (ℓ, k, j[ℓ, k]), etc.
Remark 3.1.To give an interpretation of the stopping criteria in Step (i.b.II) and Step (i.d) of Algorithm A, we note the following: Since the algebraic solver is contractive (2.9), the term ⦀u k,j ℓ − u k,j−1 ℓ ⦀ provides a posteriori error control of the algebraic error ⦀u k,⋆ ℓ − u k,j ℓ ⦀, i.e., Moreover, for sufficiently small λ alg > 0 and ongoing Zarantonello iterations, also the perturbed Zarantonello symmetrization is a contraction; see Lemma 5.1 below.With the same reasoning as for the algebraic solver, the term ⦀u ).With this understanding and the interpretation that the error estimator η ℓ (u k,j ℓ ) controls the discretization error ⦀u ⋆ − u ⋆ ℓ ⦀ (which is indeed true for u k,j ℓ = u k,j ℓ ), the heuristics behind the stopping criteria is as follows: We stop the algebraic solver in Algorithm A(i.b.II) provided that the algebraic error ⦀u k,⋆ ℓ −u k,j ℓ ⦀ is of the level of the discretization error plus the symmetrization error.Moreover, we stop the (perturbed) Zarantonello symmetrization in Algorithm A(i.d) provided that the symmetrization error ℓ ⦀ is of the level of the discretization error.Up to the factors λ alg and λ sym , this ensures that all three error sources of ⦀u ⋆ − u k,j ℓ ⦀ are equibalanced.For the analysis of Algorithm A, we recall that the set Q from (1.7) is given by ℓ is used in Algorithm A}.Together with this set, we define Note that these definitions are consistent with that of Algorithm A, but also cover the cases that the ℓ-loop, the k-loop, or the j-loop in the algorithm do not terminate, respectively.We note that formally where the latter case is excluded by Lemma 3.2.
On Q, we define a total order by Our first observation is that the algebraic solver in the innermost loop of Algorithm A always terminates.
Proof.Let (ℓ, k, 0) ∈ Q.We argue by contradiction and assume that the stopping criterion in Algorithm A(i.b.II) always fails and hence j[ℓ, k] = ∞.By assumption (2.9), the algebraic solver is contractive and hence convergent with limit u k,⋆ ℓ Moreover, by failure of the stopping criterion in Algorithm A(i.b.II), we thus obtain that ℓ by uniqueness of the fixpoint.In particular, the initial guess is already the exact solution of the linear Zarantonello system and hence the algebraic solver guarantees that u k,j ℓ = u k,⋆ ℓ for all j ∈ N 0 .Consequently, the stopping criterion in Algorithm A(i.b.II) will be satisfied for j = 1.This contradicts our assumption, and hence we conclude that j[ℓ, k] < ∞. □ Remark 3.3.For the mathematical tractability, we formulated Algorithm A in a way that #Q = ∞.Any practical implementation will aim to provide a sufficiently accurate approximation u k,j ℓ in finite time.More precisely, Algorithm A will then be terminated after Algorithm A(i.b.II) if where τ > 0 is a user-specified tolerance.For τ = 0, finite termination yields that u by uniqueness of the fixpoint of the contractive solver and the contractive Zarantonello symmetrization, respectively.Finally, the first summand in Remark 3.4.Up to the algebraic stopping criterion in Algorithm A(i.b.II), the AISFEM algorithm coincides with the adaptive algorithm from [HPSV21], where the (perturbed) Zarantonello iteration is employed for an adaptive iteratively linearized finite element method for the solution of an energy minimization problem with strongly monotone nonlinearity in the corresponding Euler-Lagrange equations.However, the present analysis is much more refined than that of [HPSV21]: (i) To guarantee full linear convergence, [HPSV21, Theorem 4] requires θ sufficiently small, λ sym sufficiently small with respect to θ, and λ alg sufficiently small with respect to λ sym .In contrast, the present analysis proves full linear convergence for arbitrary 0 < θ ≤ 1 and 0 < λ sym ≤ 1, and only requires λ alg to be sufficiently small to preserve the contraction of the perturbed Zarantonello iteration (see Lemma 5.1 below in comparison to [HPSV21, Lemma 6]).
(ii) Despite the linear model problem, our analytical setting is more involved: the compact perturbation in (2.1) prevents the use of energy arguments that guarantee a Pythagorean-type identity in terms of the energy error (see, e.g., [HPSV21;HPW21]).Instead, we first need to exploit a priori convergence of Algorithm A (see Lemma 5.3) to deduce a quasi-Pythagorean estimate in Lemma 5.4, which then allows proving linear convergence (Theorem 4.1).As a consequence (and beyond the results of [HPSV21]), this finally yields that, for arbitrary θ and λ sym , the convergence rates with respect to the number of the degrees of freedom and with respect to the overall computational work coincide (Corollary 4.2).
The following proposition provides a computable upper bound for the energy error ⦀u ⋆ − u k,j ℓ ⦀.Since Algorithm A follows the structure of [HPSV21, Algorithm 1], the proof can be obtained analogously to [HPSV21, Proposition 2] and is thus omitted here.

Main results
In the following, we formulate the main results of the present work.We refer to Section 5 for the proofs and Section 6 for numerical experiments, which underline these theoretical results.First, recall from (2.7) that a sufficiently small parameter δ > 0 ensures contraction of the Zarantonello mapping and hence with 0 < q sym < 1.The following theorem states full linear convergence of the quasi-error.
Theorem 4.1 (full linear convergence of AISFEM).Suppose that δ > 0 is sufficiently small and that the estimator satisfies (A1)-(A3).Choose λ ⋆ alg > 0 depending only on q alg from (2.9) and q sym from (4.1) such that 0 < q sym := q sym + 2 Then, for arbitrary 0 < θ ≤ 1 and 0 < λ sym ≤ 1, there exists 0 < λ ′ alg ≤ λ ⋆ alg such that Algorithm A, for all 0 < λ alg ≤ λ ′ alg , guarantees full linear convergence: There exist constants C lin > 0 and 0 < q lin < 1 as well as an index ℓ 0 ∈ N 0 with ℓ 0 ≤ ℓ such that the quasi-error The constants C lin and q lin as well as the index ℓ 0 depend only on C stab , C rel , q red , q sym , q alg , θ, λ sym , λ alg , and While the proof of Theorem 4.1 is postponed to Section 5.5, we shall immediately prove the following important consequence of Theorem 4.1: Algorithm A guarantees that rates with respect to the number of degrees of freedom coincide with rates with respect to the overall computational cost.
Corollary 4.2.Let s > 0. Under the assumptions of Theorem 4.1, the output of Algorithm A guarantees that This yields the equivalence Proof.The lower bound in (4.5) is obvious.To prove the upper bound, without loss of generality, we may assume that M (s) < ∞.By definition of M (s), it follows that For |ℓ, k, j| ≥ |ℓ ′ , k ′ , j ′ | and ℓ ′ ≥ ℓ 0 , full linear convergence (4.4) can be rewritten as The geometric series yields that Rearranging this estimate, we see that Taking the supremum over all (ℓ, k, j) ∈ Q with ℓ ≥ ℓ 0 , we prove the second estimate in (4.5).Moreover, i.e., the sets over which we compute the suprema in (4.5)-(4.6)differ only by finitely many index triples.This and (4.5) thus prove the equivalence in (4.6).□ To present our second main result on quasi-optimal computational cost, we first introduce the notion of approximation classes.For T ∈ T and s > 0, define ‖u ⋆ ‖ As(T ) := sup with u ⋆ opt and η opt denoting the exact discrete solution and the estimator on the optimal triangulation T opt ∈ T N (T ), respectively.When (4.9) is finite, this means that a decrease of the error plus estimator with rate s is possible along optimal meshes obtained by refining T .
Theorem 4.3 (optimal computational complexity).Suppose that δ > 0 is sufficiently small and that the estimator satisfies where Choose 0 < λ sym < λ ⋆ sym sufficiently small such that Then, for any 0 < λ alg ≤ λ ′ alg with λ ′ alg > 0 from Theorem 4.1, Algorithm A guarantees, for all s > 0, that sup where ℓ 0 ∈ N is the index from Theorem 4.1.The constant c opt > 0 depends only on C Céa = L/α, C stab , C rel , s, and the use of NVB refinement; the constant C opt > 0 depends only on C stab , C drel , C mark , C Céa = L/α, C ′ rel , C lin , q lin , #T ℓ 0 , q red , λ sym , q sym , θ, s, and the use of NVB refinement.In particular, this proves the equivalence which yields optimal complexity of Algorithm A.
The proof is postponed to Section 5.6.

Contraction of perturbed Zarantonello symmetrization.
Recall that for δ < 2 δ ⋆ = 2 α/L 2 , the Zarantonello mapping is a contraction (2.7).However, Algorithm A does not compute u k,⋆ ℓ ) exactly, but relies on an approximation u k,j ℓ ≈ u k,⋆ ℓ .The next lemma states that, for a sufficiently small stopping parameter λ alg > 0 in Algorithm A, the Zarantonello symmetrization remains a contraction under this perturbation (up to the final iteration).Its proof essentially follows along the lines of [HPSV21, Lemma 6].However, the present work considers a stopping criterion of the algebraic solver in Algorithm A(i.b.II) which allows to choose λ alg independently of λ sym .
Lemma 5.1.Let λ ⋆ alg > 0 and 0 < q sym < 1 as in Theorem 4.1.Then, for all stopping parameters 0 < λ alg ≤ λ ⋆ alg and λ sym > 0, it holds that (5.1) . By using the triangle inequality and the contraction (4.1) of the unperturbed Zarantonello iteration, we obtain that 2) It remains to treat the algebraic error term and to show that it is sufficiently contractive.We use the contraction (2.9) of the algebraic solver, i.e., the met algebraic stopping criterion in Algorithm A(i.b.II), and the not met stopping criterion in Algorithm A(i.d) to obtain that Combining the last estimate with (5.2) and rearranging the terms lead us to This concludes the proof of (5.1).Now suppose that k = k[ℓ].By the met algebraic stopping criterion in Algorithm A(i.b.II) followed by the met stopping criterion of the Zarantonello iteration in Algorithm A(i.d), we obtain that Together with the contraction (5.3) of the algebraic solver, this yields that (5.4) By contraction (4.1) of the unperturbed Zarantonello iteration, we obtain that This concludes also the proof of (5.1 + ).□ An important consequence of the contraction (5.1) of the perturbed Zarantonello iteration is that k[ℓ] = ∞ implies that the exact solution is already discrete u ⋆ = u ⋆ ℓ ∈ X ℓ .
Lemma 5.2.Suppose that the estimator satisfies stability (A1) and reliability (A3), and that the perturbed Zarantonello iteration is contractive (5.1).Then, ℓ < ∞ implies that Proof.Since j[ℓ, k] < ∞ by virtue of Lemma 3.2, it follows for ℓ < ∞ that k[ℓ] = ∞ and hence by the not met stopping criterion in Algorithm A(i.d) that ⦀ for all k ∈ N. Since the perturbed Zarantonello iteration is convergent (see Lemma 5.1) with limit u ⋆ ℓ (and thus (u k,j ℓ ) k∈N 0 is a Cauchy sequence), we infer that This proves η ℓ (u ⋆ ℓ ) = 0, whence with reliability (A3), we conclude u ⋆ ℓ = u ⋆ .□ 5.2.A priori convergence.For general second-order linear elliptic PDEs, an a priori convergence result is required to ensure that there holds a quasi-Pythagorean estimate; see Lemma 5.4 below.
Lemma 5.3 (a priori convergence).With ℓ ∈ N 0 ∪{∞} from (3.1), define the discrete limit space X ∞ := closure ℓ ℓ=0 X ℓ .Then, there exists and it holds that (5.6) In particular, this implies Moreover, with C Céa = L/α from (2.5), there holds the Céa-type estimate (5.7) (5.8) ∞ , the Céa lemma (2.5) holds with u ⋆ being replaced by u ⋆ ∞ , and the definition of X ∞ proves that Reliability (5.8) follows from the triangle inequality, nestedness of spaces X ℓ ⊆ X ∞ , and the Céa lemma (2.5), since ).This concludes the proof.□ 5.3.Quasi-Pythagorean estimate.While symmetric PDEs satisfy a Pythagorean identity in the energy norm (with ε = 0 and ℓ 0 = 0 in (5.9) below), the situation is more involved for nonsymmetric PDEs.The following result generalizes [BHP17, Lemma 18] by considering general v ℓ ∈ X ℓ and by additionally proving the lower bound in (5.9).Moreover, it is given here in terms of the a priori limit u ⋆ ∞ .Although the proof follows essentially that of [BHP17], we include it for the sake of completeness.

Auxiliary contraction estimates. The following lemma extends [GHPS21,
Lemma 10] to the present setting with a quasi-Pythagorean estimate.

Using this with the not met stopping criterion in Algo
(5.21) In this case, we are thus led to (5.20) Up to the choice of the parameters ε and ν, this proves (5.17a) for any 0 < λ alg ≤ λ ⋆ alg .
This concludes the proof with q 2 lin := where, compared with (4.3), the quasi-error ∆ k,j ℓ has been redefined.Later, we shall conclude that indeed u ⋆ ∞ = u ⋆ so that both definitions coincide.
Step 1.In the first step, we prove that (5.25) Together with reliability (5.8) and stability (A1), the definition of ∆ k,j ℓ shows that The contraction of the (unperturbed) Zarantonello iteration (4.1) proves that Furthermore, the contraction of the algebraic solver (5.3) proves that Combining the last three estimates with the not met stopping criterion of the algebraic solver in Algorithm A(i.b.II) for 1 ≤ j < j[ℓ, k], we conclude that Finally, the triangle inequality and the contraction (5.3) imply (5.25).
Step 3. In this step, we prove that (5.28) for all (ℓ, k, j) ∈ Q, where the hidden constants depend only on ν.Together with (5.26) from Step 2, it thus only remains to prove To this end, let (ℓ, k, j) ∈ Q with k ≥ 1.From contraction (4.1) of the unperturbed Zarantonello symmetrization and nested iteration u k,0 ℓ = u k−1,j ℓ , we get that ).
The Céa lemma (5.7) proves that Combining the last two estimates, we arrive at This concludes the proof of (5.28).
Step 5.For (ℓ, k, j) ∈ Q with ℓ ≥ ℓ 0 , the preceding steps show that With linear convergence (5.17) of Λ k ℓ from Lemma 5.5 and the geometric series, we thus see that (5.17) According to basic calculus (see, e.g., [CFPP14, Lemma 4.9]), this is equivalent to linear convergence with respect to the lexicographic order on Q, i.e., for all (ℓ, k, j), (ℓ where the constants C lin > 0 and 0 < q lin < 1 depend only on C sum .This also yields that and hence u ⋆ ∞ = u ⋆ .In particular, the definitions of ∆ k,j ℓ from (4.3) and (5.24) coincide.Overall, we thus conclude the proof of linear convergence (4.4).□ 5.6.Proof of Theorem 4.3.The proof of Theorem 4.3 requires the following auxiliary lemma stating that the error estimator η ℓ (u k,j ℓ ) of the inexact but available final iterate of Algorithm A is equivalent to the error estimator η ℓ (u ⋆ ℓ ) of the (unknown) exact solution u ⋆ ℓ .While the statement is similar to [HPSV21, Lemma 7], the present proof provides a minor clarification of the involved constant.

Rearranging the terms and noting that
We recall from [BHP17, Lemma 22] that for all T H ∈ T and all T h ∈ T(T H ), it holds that (5.48) Overall, we have thus shown that This concludes the proof of the upper bound in (5.33b) and hence that of (4.11).
Optimality of the iteratively symmetrized solver.Optimality of AISFEM is possible when the inherent symmetrization and algebraic procedures are treated efficiently.In Figure 3, we present the time required for our iteratively symmetrized solver compared to the Matlab built-in direct solver (backslash) of the linear system related to (1.4).We note that the displayed timings are comparing the direct solve time itself with the remaining time (including the setup of the Zarantonello system, computation of the error estimator, and mesh refinement).Hence, the presented numbers favor the built-in direct solver over the Matlab-implemented multigrid code.Nevertheless, the combination of the Zarantonello symmetrization with the optimal local multigrid solver from [IMPS22] appears to be of comparable speed to the built-in direct solver with the observation that as the dimension of the linear system increases, the backslash performance begins to degrade.Moreover, Figure 4 shows that the iteration numbers of the solver remain uniformly bounded in the levels for various choices of the parameters λ sym and θ.Note that when λ sym decreases, a higher accuracy of the Zarantonello symmetrization is required.Therefore, the iteration numbers are expected to increase with smaller λ sym as seen in Figure 4 (left).Moreover, the iteration numbers are also expected to increase as θ becomes larger.This is due to the aggressive refinement leading to hierarchies of low numbers of levels but with considerable increase in the dimension of the linear systems.This may lead to the conclusion that θ should be chosen very small in order to have less iterations per level, but studying the cumulative solver steps in Figure 4 (right) shows that this is not the best strategy.
Figure 3. Optimality of the combined iterative solver for the diffusionconvection-reaction problem on the L-shaped domain from Section 6.1.Cumulative time for the direct solve and AISFEM over the computational costs.
Parameter study of AISFEM.We now investigate which parameters yield the best contraction in the iteratively symmetrized steps A(ii)-(iii).Since the parameters depend on the contraction factors q alg from (2.9) and q sym from (4.1), we study a setting where the exact discrete solution u ⋆ ℓ to (1.4) and the exact Zarantonello solution u k,⋆ ℓ to (1.6) are computed.Then, we compute q alg (ℓ, k, j) for (ℓ, k, j) ∈ Q and define the level-wise contraction factors q alg (ℓ) as the maximum over all q alg (ℓ, k, j) for fixed ℓ ∈ N 0 and analogously for q sym .From now on, we fix the polynomial degree m = 2 and the parameters .Uniform contraction of the iterative solver for the diffusionconvection-reaction problem on the L-shaped domain from Section 6.1.Experimental contraction factors q alg , q sym and q sym for various choices of the symmetrization stopping parameter λ sym with fixed θ = 0.5 (left) and different marking parameter θ with fixed λ sym = 0.1 (right).⋆ alg for various choices of the symmetrization stopping parameter λ sym with fixed θ = 0.5 (left) and different marking parameter θ with fixed λ sym = 0.1 (right), where we emphasize the double scaling of the y-axis for λ ⋆ alg resp.q sym in both figures.
Figure 7 shows that even for a strong convection combined with a strong singularity at the origin, the adaptive algorithm recovers the optimal convergence rates −m/2 for several polynomial degrees m both with respect to the cumulative costs and computational time.
In Figure 4 we see that the number of solver steps per level ℓ behaves similarly to the diffusion-convection-reaction problem on the L-shape from Section 6.1 with an increase due to the stronger singularity.Furthermore, Figure 8 displays upper bounds on λ alg ≤ where F ∈ X ′ is a linear and continuous functional, a(•, •) is a symmetric, continuous, and elliptic bilinear form on X , and K : X → X ′ is a compact operator such that the bilinear form b(•, •) is still elliptic on X .Let ⦀ • ⦀ denote the a(•, •)-induced energy norm.For the discrete formulation b(u ⋆ ℓ , v ℓ ) = F (v ℓ ) for all v ℓ ∈ X ℓ , (7.2) we require an (abstract) inexact iterative solver with iteration map given by Φ ℓ (F ; •) : X ℓ → X ℓ that contracts the error in the energy norm, i.e., where the contraction constant 0 < q sym < 1 is independent of u 0 ℓ ∈ X ℓ .Under such assumptions and with the usual residual a posteriori error estimator η ℓ (•) (satisfying the abstract assumptions (A1)-(A4)) on nested conforming discrete spaces X ℓ ⊆ X ℓ+1 ⊂ X , the present work proves that the analysis from [GHPS21] can be generalized from symmetric PDEs (with K = 0) to the general formulation (7.1): Restricting Algorithm A to the outer ℓ-loop (for mesh refinement) and the inner k-loop (for the solver associated to Φ ℓ ), we obtain a simplified index set ℓ is computed by the simplified algorithm} (7.4) together with the canonical step counter |ℓ, k| ∈ N 0 on Q defined analogously to (3.2).Then, Lemma 5.2 (lucky non-termination of the solver), Lemma 5.3 (a priori convergence), Lemma 5.4 (quasi-Pythagorean estimate), and Lemma 5.5 (contraction of weighted discretization and solver error) hold verbatim (and the proof of Lemma 5.4 indeed relies on the compactness of K) if we replace u k,j ℓ in the given proofs by u k ℓ in the current solver setting.Therefore, we obtain full linear convergence in the spirit of Theorem 4.1: For arbitrary adaptivity parameters 0 < θ ≤ 1 and λ sym > 0, there exist constants C lin > 0 and 0 < q lin < 1 as well as an index ℓ 0 ∈ N 0 such that for all (ℓ ′ , k ′ ), (ℓ, k) ∈ Q with |ℓ ′ , k ′ | ≤ |ℓ, k| and ℓ ′ ≥ ℓ 0 , (7.5) where ∆ k ℓ := ⦀u ⋆ − u k ℓ ⦀ + η ℓ (u k ℓ ) denotes the corresponding quasi-error.In particular, also Corollary 4.2 holds verbatim with Q replaced by Q and ∆ k,j ℓ replaced by ∆ k ℓ , i.e., convergence rates with respect to the number of degrees of freedom coincide with rates with respect to the overall computational cost.Finally, it is easy to check that also Theorem 4.3 holds verbatim and proves that, for sufficiently small adaptivity parameters 0 < θ ≪ 1 and 0 < λ sym ≪ 1 in the sense of (4.10), it holds that which yields optimal complexity of the simplified algorithm.
In the current analysis, the combined Zarantonello symmetrization with a contractive SPD algebraic solver is used as one solver module to guarantee (7.3) for u k ℓ := u k,j ℓ (see Lemma 5.1, where contraction, however, only holds for 1 ≤ k < k[ℓ]), leading to all results being formulated over the triple index set Q ⊂ N 3 0 (see Section 3-4).We note that another choice for solving the arising nonsymmetric FEM systems would be preconditioned GMRES (see, e.g., [SS86;Saa03]), where an optimal preconditioner for the symmetric part would be employed.Then, it is well-known from the field-ofvalue analysis (see, e.g., [Sta97]) that the algebraic solver would satisfy a generalized contraction for the algebraic residual (in a discrete vector norm).However, the link between the algebraic residual and the functional setting appears to be open.Moreover, the a posteriori error control of the algebraic error for such a GMRES solver is still to be developed.
While these questions are left for future work, we already note some results that can be achieved along the arguments of [GHPS21]: If the solver Φ ℓ (F ; •) provides iterates (u k ℓ ) k∈N 0 satisfying only the generalized contraction together with the a posteriori error control where C sym , C ′ sym > 0 and 0 < q sym < 1 are given constants independently of u 0 ℓ ∈ X ℓ , then full linear convergence (7.5) can be proved for all 0 < θ ≤ 1 under the additional assumption that λ sym has to be sufficiently small.However, the proof of full linear convergence (7.5) for arbitrary 0 < θ ≤ 1 and arbitrary λ sym > 0 is open, while optimal complexity (7.6) for sufficiently small 0 < θ < 1 and λ sym in the sense of (4.10) remains valid (even with the same proof).

Figure 1 .
Figure 1.Schematic view of the AISFEM algorithm components.
.43) with the constant C mark ≥ 1 from Algorithm A. Recall the definition θ mark (4.10)

10 3
Figure5.Uniform contraction of the iterative solver for the diffusionconvection-reaction problem on the L-shaped domain from Section 6.1.Experimental contraction factors q alg , q sym and q sym for various choices of the symmetrization stopping parameter λ sym with fixed θ = 0.5 (left) and different marking parameter θ with fixed λ sym = 0.1 (right).