---
abstract: 'We deal with the numerical solution of linear partial differential equations (PDEs) with focus on the goal-oriented error estimates including algebraic errors arising by an inaccurate solution of the corresponding algebraic systems. The goal-oriented error estimates require the solution of the primal as well as dual algebraic systems. We solve both systems simultaneously using the bi-conjugate gradient method which allows to control the algebraic errors of both systems. We develop a stopping criterion which is cheap to evaluate and guarantees that the estimation of the algebraic error is smaller than the estimation of the discretization error. Using this criterion and an adaptive mesh refinement technique, we obtain an efficient and robust method for the numerical solution of PDEs, which is demonstrated by several numerical experiments.'
author:
- 'V[í]{}t Dolej[š]{}[í]{}'
- 'Petr Tich[ý]{}'
bibliography:
- 'ref.bib'
date: 'Received: date / Accepted: date'
title: 'On efficient numerical solution of linear algebraic systems arising in goal-oriented error estimates[^1] '
---

[example.eps]{} gsave newpath 20 20 moveto 20 220 lineto 220 220 lineto 220 20 lineto closepath 2 setlinewidth gsave .4 setgray fill grestore stroke grestore

Introduction
============

Let $\Om \subset \IR^d$ be a polygonal domain with the boundary $\gom=\partial\Om$. We consider an abstract partial differential equation in the form $$\begin{aligned}
  \label{eq:PP}
  \krl \u = f,\end{aligned}$$ where $\u:\Om \to \IR$ is the unknown solution, $\krl$ is a linear differential operator and $f$ is the right-hand side. The problem has to be accompanied by suitable boundary conditions. Further, let $\Vh$ be a finite-dimensional space of functions, ($\dim\Vh=\Nh< \infty$), where an approximation of $\u$ is sought. We note that the space $\Vh$ is not given a priori but it is generated as a sequence of spaces by a suitable mesh adaptive method. By $\uh\in\Vh$ we denote the discrete solution which approximates $u$.

In many practical application, we are not interested in the solution of itself but rather in a [*quantity of interest*]{}, which is the value of a certain, a priori known, solution-dependent [*target functional*]{} $\J(\u)$. For example, the functional $\J$ is a mean value of the solution over a subset of the computational domain or its boundary. Therefore, we need not to estimate the error $\norm{\u-\uh}{}$ in an energy norm (cf. [@AinsworthOden; @verfurth-book2]) but the error $\J(u)-\J(\uh)$. We assume that $J$ is a linear functional.

The necessity to estimate the error of the target functional gives rise to the [*goal-oriented error estimates*]{}, cf. the pioneering works summarized in [@RannacherBook; @BeckerRannacher01; @GileSuli02]. This approach was further developed for many problems, let us mention [@OdenPrudhomme_CAMWA01; @Korotov_JCAM06; @SolDemko04] dealing with linear elliptic problems, [@KuzminMoller_JCAM10] dealing with steady linear hyperbolic problems. Moreover, extensions to nonlinear problem were presented, e.g., in [@Rey2Gosselet_CMAME14; @Rey2Gosselet_CMAME15] for elasticity problems and, e.g., in [@LoseilleDervieuxAlauzet_JCP10; @Richter_IJNMF10; @HH06:SIPG2; @GeHaHo09; @Hartman08; @BalanWoopenMay16] for computational fluid dynamics, for a survey see [@FidkowskiDarmofal_AIAA11].

The goal-oriented error estimates require, except the solution of the original ([*primal*]{}) problem , also to solve the [*dual*]{} (or adjoint) problem $$\begin{aligned}
  \label{eq:DP}
  \krl^* \z = \J,\end{aligned}$$ where $\krl^*$ is the [*dual operator*]{} to $\krl$, $\J$ is the target functional and $\z:\Om\to\IR$ is the dual solution. The problem has to be accompanied by suitable boundary conditions as well. It is necessary to approximate the solution of by $\zhp\in\Vhp$ where $\Vhp$ is a richer space than $\Vh$. One possibility is to discretize directly on $\Vhp$ (e.g., [@Korotov_JCAM06]) but then we have to solve two different algebraic problems. Another way is to discretize on $\Vh$ and define $\zhp:= \krr(\zh)$, where $\zh\in\Vh$ is a numerical approximation of and $\krr:\Vh\to\Vhp$ is a suitable higher order reconstruction (e.g., [@Richter_IJNMF10; @CarpioPrietoBermejo_SISC13]).

The advantage of the latter approach is that the discretizations of the primal and dual problems and using the same space $\Vh$ are equivalent to two mutually transposed linear algebraic systems which can be beneficial in practical solutions. Namely, we obtain $$\begin{aligned}
  \label{i0}
  \mA\bx = \bb\qquad\mbox{and}\qquad \mA^{\T}\by=\bc,\end{aligned}$$ where $\mA\in\IR^{\Nh\times\Nh}$ is the matrix arising from the discretization of $\krl$ from , $\mA^{\T}$ is the transpose matrix of $\mA$, $\bb\in\IR^{\Nh}$ and $\bc\in\IR^{\Nh}$ represent the discretization of the right-hand sides of and , respectively, and $\bx\in\IR^{\Nh}$ and $\by\in\IR^{\Nh}$ are the vectors corresponding to the discrete solutions $\uh\in\Vh$ and $\zh\in\Vh$ of and , respectively.

In many situations, it is advantageous to solve systems iteratively since (i) approximate solutions satisfying the prescribed tolerance are sufficient as an output of the computation; (ii) very good initial approximations of $\bx$ and $\by$ are typically available from previous level of mesh adaptation. Let $\uhk\in\Vh$ and $\zhk\in\Vh$ denote functions corresponding to the $k$-th iterations $\bxk$ and $\byk$, respectively. Then the computed approximations $\uhk$ and $\zhk$ are influenced also by the algebraic error arising from the inexact solution of systems . Let us note that even if a direct solver for is used then these systems are not solved exactly (only on the level of machine accuracy) and the computed approximations suffers from algebraic errors too, see [@ArioliLiiesenMiedlerStrakos13].

Having only approximate solution of , the Galerkin orthogonality of the error is violated and many standard a posteriori error estimation techniques can not be used. Therefore, the algebraic error has to be taken into account in a posteriori error analysis. Additionally, the algebraic error estimation is important for the setting of a suitable stopping criterion for iterative solvers and therefore for the optimization of the computational costs. The algebraic error estimates in the framework of energy norms were treated, e.g., in [@Arioli_NM04; @Picasso_CNME09; @VohralikStrakos10] for conforming finite element method in combination with the conjugate gradient method. In the framework of goal-oriented error estimates, a possible violation of Galerkin orthogonality was taken into account, e.g., in [@KuzminKorotov_MCS10]. The influence of the algebraic error in goal-oriented error estimates was considered in [@Meidner2009Goal] for iterative algebraic solvers (including multigrid methods). The computational error was expressed and estimated as a sum of the [discretization]{} and [algebraic]{} errors. The primal and dual algebraic problems were solved alternatively and after some number of iterations (or one multigrid cycle), the algebraic and discretization components of the errors were estimated. Then the iterative process was stopped when the algebraic error estimate was sufficiently smaller than the discretization error estimate. We also mention the recent paper [@MallikVohralikYousef_JCAM20] using a completely deferent approach for elliptic problems. This technique is based on $H(\mathrm{div})$-conforming flux reconstructions and $H^1$-conforming potential reconstructions and yields to a guaranteed upper error bound.

In [@DolRoskovec_AM17], we employed an approach similar to [@Meidner2009Goal] and studied the convergence of the error estimator with respect to algebraic iterations of the primal and dual problems. We observed a delay in the convergence for insufficiently accurate resolution of the primal or the dual problem. We arrived to the conclusion that it is difficult to control the efficiency of the computational process when the primal and dual problems are solved alternatively.

Therefore, in this paper, we develop a technique to solve the primal and dual problems simultaneously using the [*bi-conjugate gradient*]{} (BiCG) method. At each BiCG iteration, the approximations of $\uh$ and $\zh$ are available and the estimate of the algebraic error can be computed and the accuracy and efficiency can be controlled. Motivated by [@StTi2011], we proposed a stopping criterion for the BiCG solver which is very cheap for evaluation and in contrary to [@Meidner2009Goal; @DolRoskovec_AM17], it does not require the evaluation of the discretization error estimator. The use of the BiCG solver for the primal and dual problems with the proposed criterion in a combination the mesh adaptation leads to an efficient numerical method for the solution of .

Since the presented approach can be applied to with a general linear operator $\krl$ discretized by any numerical scheme based on a variational formulation, we express the problem considered and its numerical approximation only in the abstract form. However, in order to demonstrate the applicability of this technique, we present the numerical solution of purely elliptic and convection-diffusion problems by the $hp$-adaptive discontinuous Galerkin method (DGM) on possibly anisotropic meshes. The content of the rest of the paper is the following. In Section \[sec:GO\], we shortly summarize the goal-oriented error estimate technique including two variants of the expression of the algebraic errors. In Section \[sec:algeb\], we introduce the algebraic representation of the discrete problems and discuss several ways for the evaluation of the quantity of interest and their corresponding error estimates. These possibilities are numerically tested in Section \[sec:numerF\], where we solve the Laplace problem on fixed meshes. Moreover, in Section \[sec:adapt\], we describe the mesh adaptation process and proposed new stopping criteria for the BiCG solver. Their computational performance is demonstrated in Section \[sec:numerA\] for the Laplace and convection-dominated problems. We conclude with several remarks and discuss a possible extension of this technique for nonlinear problems in Section \[sec:sum\].

Framework for the goal-oriented error estimates {#sec:GO}
===============================================

We briefly recall a general framework for the goal-oriented error estimates. More details can be found, e.g., in [@BeckerRannacher01; @GileSuli02].

Primal problem
--------------

Let the [*weak formulation*]{} of the primal problem be given by $$\begin{aligned}
  \label{eq:PPw}
  a( \u, \vp) = \ell(\vp)\qquad \forall \vp\in \V,\end{aligned}$$ where $\u\in \V$ is a weak solution, $a(\cdot,\cdot):\V\times \V\to\IR$ is a bilinear form, $\ell(\cdot):\V\to\IR$ is a linear form and $V$ is a Hilbert space. We assume that is well-posed, i.e., it admits a unique weak solution.

For the numerical approximation of , let $\Vh,\ h>0$ be a finite element space of functions defined on $\Om$, typically piecewise polynomial functions related to the partition of $\Om$ onto a set of finite elements $\Th$. Moreover, let $\Wh$ be a functional space such that $\V\subset \Wh$ and $\Vh\subset\Wh$. For conforming finite element methods ($\Vh\subset\V$), we simply put $\Wh:=\V$. However, for nonconforming methods (the case $\Vh\not\subset\V$), the choice of $\Wh$ is more delicate. E.g., for [*discontinuous Galerkin method*]{}, we employ the so-called [*broken Sobolev space*]{} and put $$\begin{aligned}
  \label{HsTh}
  \Wh:= H^s(\Th) = \{ v \in L^{2}(\Om); \, \restr{v}{K} \in H^{s}(K) \,
  \forall K \in \Th\},\end{aligned}$$ where $\Th$ is a mesh partition of $\Om$ and $s>0$ is a suitable Sobolev index, e.g., for the second-order operator $\krl$, we have $s=2$.

Let $$\begin{aligned}
  \label{ah}
  \ah:\Wh\times\Wh\to\IR\qquad \mbox{and} \qquad  \ell_h:\Wh\to\IR\end{aligned}$$ be a bilinear and a linear forms corresponding to the discretization of the left-hand and right-hand sides of , respectively, by a particular numerical method.

We say that $\uh\in\Vh$ is the [*discrete solution*]{} of the [*primal problem*]{} if $$\begin{aligned}
  \label{eq:PPh}
    \ah(\uh,\vp_h) = \ell_h(\vp_h)\qquad \forall \vp_h\in\Vh.\end{aligned}$$ We assume that the numerical scheme is [*consistent*]{}, i.e., $$\begin{aligned}
  \label{eq:cons}
  \ah(\u,\vp) = \ell_h(\vp)\qquad \forall \vp\in\Wh\end{aligned}$$ where $\u\in\V$ is the weak solution of . This implies the [*Galerkin orthogonality of the error*]{} of the primal problem $$\begin{aligned}
  \label{eq:GO}
  \ah(\uh - u, \vp_h) = 0 \qquad\forall \vp_h\in\Vh.\end{aligned}$$ Finally, we define the [*residual of the primal problem*]{} by $$\begin{aligned}
  \label{res}
  \res{\uh}{\vp}:= \ell_h(\vp) - \ah(\uh, \vp) = \ah(\u-\uh,\vp) ,\quad \vp\in\Wh,\end{aligned}$$ where the last equality follows from the consistency .

Quantity of interest and the dual problem
-----------------------------------------

As mentioned in the introduction, we are interested in a [*sufficiently accurate*]{} approximation of the [*quantity of interest*]{} $J(u)\in\IR$, where $$\begin{aligned}
  \label{Ju}
  \J:\Wh\to \IR\end{aligned}$$ is a linear functional. Typically, it is defined as a weighted mean value of $u$ over the computational domain $\Om$ or its boundary $\gom$.

In order to estimate the error $\J(\u) - J(\uh)$, we consider the [*adjoint*]{} (or [*dual*]{}) [*problem*]{} and its discretization. We say that $\zh$ is the [*discrete solution*]{} of the [*dual problem*]{} if $$\begin{aligned}
  \label{eq:DPh}
  \ah(\psi_h, \zh) = \J(\psi_h) \qquad \forall \psi_h\in\Vh,\end{aligned}$$ where $\ah$ and $\J$ are given by and , respectively. Moreover, we assume that the numerical scheme is [*adjoint consistent*]{}, i.e., $$\begin{aligned}
  \label{eq:consD}
  \ah(\psi, \z) = \J(\psi)\qquad \forall \psi\in\Wh,\end{aligned}$$ where $\z$ is the weak solution of the dual problem . This implies the [*Galerkin orthogonality of the error*]{} of the dual problem $$\begin{aligned}
  \label{eq:GOD}
  \ah(\psi_h, \zh - \z) = 0 \qquad\forall \psi_h\in\Vh.\end{aligned}$$ Finally, we define the [*residual of the dual problem*]{} by $$\begin{aligned}
  \label{resD}
  \resD{\zh}{\psi}:= \J(\psi) - \ah(\psi,\zh) = \ah(\psi, \z-\zh),\quad \psi\in\Wh,\end{aligned}$$ where the last equality follows from the adjoint consistency .

The problems and represent two linear algebraic systems whose efficient solution by an iterative solver is developed in §\[sec:algeb\]. Due to iterative and rounding errors, the “exact” discrete solutions $\uh$ and $\zh$ are not available, we have only their approximation $\uhk\in\Vh$ and $\zhk\in\Vh$, $k=1,2,\dots$. However, they do not fulfil the Galerkin orthogonalities and .

Abstract goal-oriented error estimates {#sec:GOEE}
--------------------------------------

Obviously, for adjoint consistent discretization we have (due to and ) the equivalence between the quantity of interest $\J(\u)$ and the right-hand side of evaluated for the dual solution $\z,$ i.e. $$\begin{aligned}
  \label{eq:eqi}
 \ell_h(\z) = \ah(\u, \z) = \J(\u).\end{aligned}$$ Similarly, from and , we obtain the discrete variant of as $$\begin{aligned}
  \label{eq:eqi_h}
 \ell_h(\zh) = \ah(\uh, \zh) = \J(\uh),\end{aligned}$$ where $\uh$ and $\zh$ are the discrete solutions of the primal and dual problems, respectively. Consequently, we have an [*error equivalence*]{} $$\begin{aligned}
  \label{eq:eqi_err}
  \J(\u -\uh) =  \J(\u) - \J(\uh) = \ell_h(\z) - \ell_h(\zh) = \ell_h(\z -\zh).\end{aligned}$$ Therefore, the difference $\ell_h(\z -\zh)$ can be used as an error estimate of the quantity of interest as well.

First, we present the primal and dual error identities for the error of the quantity of interest for the algebraically exact discrete solution $\uh$ and $\zh$ of and , respectively. Using the adjoint consistency , the Galerkin orthogonality and relation , we get the [*primal error identity*]{} $$\begin{aligned}
 \label{eq:JP}
 \J(\u - \uh) &=  \ah(\u - \uh, \z)=  \ah(\u - \uh, \z - v_h)  \\
   &= \res{\uh}{\z -  v_h}  \quad \forall v_h \in \Vh, \notag\end{aligned}$$ where $\res{\zh}{\cdot}$ denotes the residual of the primal problem given by .

Similarly, exploiting in addition the Galerkin orthogonality and relation , we get the [*dual error identity*]{} $$\begin{aligned}
\label{eq:JD}
\J(\u - \uh) &=  \ah(\u - \uh, \z - \zh)
=   \ah(\u - w_h, \z - \zh ) \\
&=   \resD{\zh}{\u -  w_h}  \quad \forall w_h \in \Vh, \notag\end{aligned}$$ where $\resD{\zh}{\cdot}$ denotes the residual of the dual problem given by .

Abstract goal-oriented error estimates including algebraic errors {#sec:GOEE_alg}
-----------------------------------------------------------------

As mentioned above, the discrete solutions $\uh$ and $\zh$ fulfilling and , respectively, are non-available but we have only $\uhk\in\Vh$ and $\zhk\in\Vh$, $k=1,2,\dots$ denoting their approximations given by an algebraic iterative solver. Hence, the error representations and are useless and we have to estimate the error $J(\u-\uhk)$.

In the same way as in [@Meidner2009Goal; @DolRoskovec_AM17], we employ the adjoint consistency , identity $\z = \z-v_h + v_h,\ v_h\in\Vh$ and relation , we get the [*primal error identity*]{} including [*algebraic errors*]{} $$\begin{aligned}
 \label{eq:JP_alg0}
 \J(\u - \uhk) &=  \ah(\u - \uhk, \z)=  \ah(\u - \uhk, \z - v_h) + \ah(\u - \uhk, v_h)
 \\
 &   = \res{\uhk}{\z -  v_h} + \res{\uhk}{v_h} \quad \forall v_h \in \Vh,\notag\end{aligned}$$ where $\res{\uhk}{\cdot}$ is given by . Inserting $v_h:=\zhk$ in , we obtain $$\begin{aligned}
 \label{eq:JP_alg}
 \J(\u - \uhk) = \res{\uhk}{\z -  \zhk} + \res{\uhk}{\zhk}=: \errSk + \errAk,\end{aligned}$$ where the quantity $\errSk=\res{\uhk}{\z -  \zhk}$ represents the discretization error of the primal problem since it coincides with for $\uhk=\uh$. Further, the quantity $\errAk=\res{\uhk}{\zhk}$ represents the algebraic error of the primal problem since it vanishes for $\uhk=\uh$. Let us note that and hold without the Galerkin orthogonality .

In order to derive the analogue of including algebraic error, we take into account the error equivalence . The quantity $\ell_h(\z - \zhk)$ exhibits the analogue of the error $ \J(\u - \uhk)$ since both quantities are equal for $\uhk=\uh$ and $\zhk=\zh$. However, if $\uhk$ and $\zhk$ are arbitrary iterations then, in general, $$\begin{aligned}
  \label{eq:equi_alg}
  \J(\u - \uhk) \not= \ell_h(\z - \zhk).\end{aligned}$$ Nevertheless, we show in Section \[sec:algeb\] (cf. ) that if $\uhk$ and $\zhk$ are obtained by the bi-conjugate gradient (BiCG) method with vanishing initial approximations (in the exact arithmetic) then $$\begin{aligned}
  \label{eq:equi_algB}
  \J(\u - \uhk) =\ell_h(\z - \zhk)\qquad \mbox{ for } k=1,2,\dots.\end{aligned}$$

Now, using the consistency , identity $\z = \z-w_h + w_h,\ w_h\in\Vh$ and relation , we get the [*dual error identity*]{} including [*algebraic errors*]{} $$\begin{aligned}
  \label{eq:JD_alg0}
  \ell_h(\z - \zhk) &=  \ah(\u , \z-\zhk)=  \ah(\u - w_h, \z - \zhk) + \ah( w_h, \z - \zhk)
  \\
  &   = \resD{\zhk}{\u -  w_h} + \resD{\zhk}{w_h} 
  \quad \forall w_h \in \Vh,\notag\end{aligned}$$ where $\resD{\zhk}{\cdot}$ is given by . Putting $w_h:=\uhk$ in , we obtain $$\begin{aligned}
  \label{eq:JD_alg}
  \ell_h(\z - \zhk) = \resD{\zhk}{\u -  \uhk} + \resD{\zhk}{\uhk} =: \errDSk + \errDAk,\end{aligned}$$ where, similarly as the quantities $\errSk$ and $\errAk$ in , the quantities $\errDSk$ and $\errDAk$ represent the discretization and algebraic error of the dual problem. Let us note that and hold without the Galerkin orthogonality .

Computable goal-oriented error estimates {#sec:GOEE_com}
----------------------------------------

Whereas the algebraic errors $\errAk$ and $\errDAk$ from and , respectively, are computable quantities, the discretization errors $\errSk$ and $\errDSk$ require the knowledge of the exact dual and primal and solutions $\z$ and $\u$, respectively. In practical computations, they should be approximated by a higher-order reconstruction denoted here by $$\begin{aligned}
  \label{krr}
  \huh=\krr(\uhk),\qquad \hzh=\krr(\zhk),\end{aligned}$$ where $\uhk\in\Vh$ and $\zhk\in\Vh$ denote approximations of $\uh$ and $\zh$, $\krr:\Vh\to\Vhp$ denotes a reconstruction operator and $\Vhp$ is a “richer” finite dimensional space, for examples, see the papers cited in Introduction. Replacing $\u$ and $\z$ in and by $\huh$ and $\hzh$, we obtain the computable approximations $\etaSk$ and $\etaDSk$ of the discretization errors $\errSk$ and $\errDSk$ by $$\begin{aligned}
  \label{etas}
  &\errSk = \res{\uhk}{\z -  \zhk} \approx \res{\uhk}{\hzh -  \zhk} =: \etaSk, \\
  &\errDSk = \resD{\zhk}{\u -  \uhk} \approx \res{\zhk}{\huh -  \uhk} =: \etaDSk,\qquad
  k=1,2,\dots.
  \notag\end{aligned}$$ The quantities $\etaSk$ and $\etaDSk$ are called the estimates of the primal and dual discretization errors, respectively.

The approximations together with and lead to the computable estimate of the total error by $$\begin{aligned}
  \label{EEP}
  \J(\u - \uhk) & \approx   \etaSk + \etaAk, \qquad 
       \ell_h(\z - \zhk)  \approx \etaDSk + \etaDAk, \quad k=1,2,\dots,\end{aligned}$$ where, in order to have a consistent notation, we put $$\begin{aligned}
  \label{etasA}
  \etaAk := \errAk = \res{\uhk}{\zhk}, \quad
  \etaDAk := \errDAk = \resD{\zhk}{\uhk}.\end{aligned}$$ The terms $\etaAk$ and $\etaDAk$ are equal to the [*algebraic errors*]{} of the primal and dual problems and , respectively, and they are independent on the used higher-order reconstruction $\krr$ from .

\[rem:guaranteed\] The approximations do not give a guaranteed error estimate since, as usual, we neglected the terms $\res{\uhk}{\z -  \hzh} $ and $\res{\zhk}{\u -  \huh} $. Hence, the guaranteed error estimate requires an estimation of these terms. We refer to [@NochettoVeeserVerani_IMAJNA09; @Ainsworth2012Guaranteed] where this problem was treated for symmetric elliptic problems discretized by the conforming finite element method. This is also a subject of our further research.

An alternative representation of the algebraic error {#sec:alter}
----------------------------------------------------

In , we expressed the discretization and algebraic parts of the computational error using the approach from [@Meidner2009Goal], namely $$\begin{aligned}
 \label{JP1}
 \J(\u - \uhk) = \res{\uhk}{\z -  \zhk} + \res{\uhk}{\zhk} = \errSk +\errAk.\end{aligned}$$ On the other hand, it is possible to decompose the total error $\J(\u - \uhk)$ in the alternative way as $$\begin{aligned}
  \label{JP2}
  \J(\u - \uhk) = \J(\u - \uh) + \J(\uh - \uhk)=: \estSk + \estAk.\end{aligned}$$ The first term in is the discretization error given by and it satisfies $$\begin{aligned}
  \label{JP3}
   \estSk = \J(\u - \uh) =  \res{\uh}{\z -  v_h}, \quad  v_h \in \Vh\end{aligned}$$ and it is independent of $k$. Obviously, although $\errSk$ in and $\estSk$ in both represent the discretization error, they differ. Moreover, $\estSk$ cannot be evaluated even if we replace $\z$ by $\hzh$ since $\uh$ is unavailable.

The second term $\estAk$ in represents the algebraic error. It cannot be expressed in a residual form but in Section \[sec:algeb\], we present a technique which is able to estimate this term in a cheap way.

Hence, we have two representations of the algebraic error, the first one from given by $\etaAk=\res{\uhk}{\zhk}$ and the second one from given by $\estAk:=\J(\uh - \uhk)$. Both representations are “algebraically consistent” which means that if $\uhk\to \uh$ and $\zhk\to \zh$ for $k\to \infty$ then $$\begin{aligned}
  \label{JP4}
  \etaAk \to 0 \quad \mbox{ for } k\to \infty
  \qquad \mbox{ and }\qquad 
  \estAk \to 0 \quad \mbox{ for } k\to \infty.\end{aligned}$$ However, the speed of convergence for $\etaAk$ and $\estAk$ can differ substantially, see numerical experiments in Section \[sec:numerA\].

Moreover, for $\uhk\to \uh$ and $\zhk\to \zh$, both discretization errors $\errSk$ and $\estSk$ are closer and closer since (cf. and ) $$\begin{aligned}
  \label{JP4a}
  | \errSk - \estSk| = | \estAk - \errAk| \leq | \estAk| + |\errAk| \to 0.\end{aligned}$$

Finally, we present the dual analogue of – in the form $$\begin{aligned}
  \label{JP5}
  \ell_h(\z - \zhk) &= \ell_h(\z - \zh) + \ell_h(\zh - \zhk) 
  =: \estDSk + \estDAk,\end{aligned}$$ where the term $\estDAk$ can be estimated in the same manner as $\estAk$ in , cf. Section \[sec:algeb\].

Solution of primal and dual discretized problems {#sec:algeb}
================================================

In this section, we introduce the algebraic representation of the goal-oriented error estimates from the previous section, present the BiCG method allowing a simultaneous solution of the primal and dual problems and discuss several possibilities algebraic errors estimates and stopping criteria for the iterative solver.

Algebraic representation {#sec:alg1}
------------------------

Let $\{\vp_i:\Om\to\IR,\ i=1,\dots, \Nh \}$ be a basis of the finite-dimensional space $\Vh$. We define the matrix $\mA\in\IR^{\Nh\times\Nh}$ by $$\begin{aligned}
  \label{mA}
  \mA = \{\mA_{i,j}\}_{i,j=1}^{\Nh},\quad \mA_{i,j} := \ah(\vp_i,\vp_j),\ i,j=1,\dots, \Nh,\end{aligned}$$ where $\ah$ is the bilinear form defined by . Then the primal and dual discrete problems and are equivalent to the solution of two linear algebraic systems $$\begin{aligned}
  \label{AA1}
  \mA\bx = \bb\qquad\mbox{and}\qquad \mA^{\T}\by=\bc,\end{aligned}$$ respectively, where $\bx\in\IR^\Nh$ and $\by\in\IR^\Nh$ are the algebraic representation of the primal and dual solutions given by $$\begin{aligned}
  \label{AA2}
  & \bx = (\x_1,\dots, \x_{\Nh})^{\T} \ \leftrightarrow \uh =\sum\nolimits_{i=1}^\Nh \x_i \vp_i, \\
  \quad 
  \mbox{and} \quad &   \by = (\y_1,\dots, \y_{\Nh})^{\T} \ \leftrightarrow
  \zh =\sum\nolimits_{i=1}^\Nh \y_i \vp_i,
  \notag\end{aligned}$$ respectively, and $\bb\in\IR^\Nh$ and $\bc\in\IR^\Nh$ are the algebraic representation of the right-hand sides of the primal and dual problems given by $$\begin{aligned}
  \label{AA3}
  & \bb = (\b_1,\dots, \b_{\Nh})^{\T} \ \leftrightarrow \b_i = \ell_h(\vp_i),\ i=1,\dots,\Nh,
  \\ 
   \mbox{and} \quad
   &  \bc = (\c_1,\dots, \c_{\Nh})^{\T} \ \leftrightarrow \c_i =\J(\vp_i),\quad i=1,\dots,\Nh,
   \notag\end{aligned}$$ respectively. Using –, we obtain the equivalencies $$\begin{aligned}
  \label{AA4}
  \J(\uh)  &   = \sum_{i=1}^\Nh \x_i \J(\vp_i)
  = \bcT \bx = \bcT \mA^{-1} \bb 
   = \by^{\T} \bb 
    = \sum_{i=1}^{\Nh} \y_i \ell_h(\vp_i)
    = \ell_h(\zh), $$ which exhibits the algebraic analogue of .

Similarly, as in let $\bxk$ and $\byk$ be the algebraic analogues of the approximations $\uhk$ and $\zhk$, respectively. The corresponding residuals of are given by $$\begin{aligned}
  \label{AA5}
  \brk := \bb - \mA \bxk =\mA (\bx - \bxk) \quad \mbox{ and } \quad
  \bsk := \bc - \mA^{\T} \byk = \mA^{\T} (\by -  \byk).\end{aligned}$$ By a standard manipulation, we derive the following correspondence between the discretization and its algebraic representation which will be used in next paragraphs. Similarly as in , we have $$\begin{aligned}
  \label{aa1}
  \J(\uhk)  & = \bcT \bxk, \qquad 
  \ell_h(\zhk)   = \bykT \bb.\end{aligned}$$ From and , we obtain $$\begin{aligned}
  \label{aa3}
  \ah(\uhk, \zhk ) & = \bykT \mA \bxk.\end{aligned}$$ Additionally, employing and , we derive $$\begin{aligned}
  \label{aa4}
  \ah(\uh - \uhk, \zh - \zhk ) & = (\by - \byk)^{\T} \mA (\bx - \bxk) = 
  \bskT \mA^{-1} \brk.\end{aligned}$$ Finally, using , , and –, we obtain $$\begin{aligned}
  \label{aa5}
  \etaAk & = \ell(\zhk) - \ah(\uhk ,\zhk)
  = \bykT \bb - \bykT \mA \bxk
  = \bykT \brk, \\
    \etaDAk & = \J(\uhk) - \ah(\uhk ,\zhk)
  = \bcT \bxk - \bykT \mA \bxk
  = \bskT \bxk. \notag\end{aligned}$$

Approximating the quantity of interest using iterates $\bxk$ and $\byk$ {#sec:eqi}
-----------------------------------------------------------------------

Let the approximations $\bxk$ and $\byk$, and the corresponding residual vectors $\brk$ and $\bsk$, computed by some iterative method for solving linear systems , be given. We introduce several possibilities, how to approximate the quantity of interest $\J(\uh) = \bcT \mA^{-1} \bb = \ell_h(\zh)$ using these vectors. For each possibility, we express the quantity of interest as a sum of two terms, where the former one is a computable value approximating the quantity of interest and the latter one is an incomputable term which represents the error of the approximation. Note that the identities derived below rely only on the relations $$\label{eqn:residuals}
  \brk = \bb - \mA \bxk \quad \mbox{ and } \quad
  \bsk = \bc - \mA^{\T} \byk.$$

-   Using , a simple manipulation gives $$\begin{aligned}
        \label{BB1}
        \bcT \mA^{-1} \bb = \bcT\bx = \bcT\bxk +  \bcT( \bx - \bxk)
        =  \underbrace{\bcT\bxk}_{\mathrm{approximation}}
        + \underbrace{\bcT \mA^{-1} \brk}_{\mathrm{error}}.
      \end{aligned}$$ The first term on the right-hand side of is computable and then it can be used for the approximation $\bcT \mA^{-1} \bb \approx \bcT\bxk$ and the second term $\bcT \mA^{-1} \brk$ represents the corresponding error. Similarly, for the dual form, we have $$\begin{aligned}
        \label{BB2}
        \bcT \mA^{-1} \bb = \by^{\T}\bb = \bykT \bb +  ( \by - \byk)^{\T}\bb
        = \underbrace{\bykT \bb}_{\mathrm{approximation}} +
        \underbrace{ \bskT\mA^{-1} \bb}_{\mathrm{error}},
      \end{aligned}$$ thus $\bykT \bb$ is a computable approximation of $\bcT \mA^{-1} \bb$ and $\bskT\mA^{-1} \bb$ the corresponding error.

-   We carry out more sophisticated manipulation, take into account , , and get $$\begin{aligned}
        \label{BB3}
        \bcT \mA^{-1} \bb & = \bcT\bxk +  \bcT \mA^{-1}\mA( \bx - \bxk) 
         =  \bcT\bxk + \by^{\T} \mA(\bx-\bxk)   \\
        & =   \bcT\bxk + \bykT (\bb- \mA\bxk) + (\by-\byk)^{\T}\mA (\bx- \bxk)
        \notag \\
        & = \underbrace{\bcT\bxk + \bykT \brk}_{\mathrm{approximation}}
        + \underbrace{\bskT\mA^{-1} \brk}_{\mathrm{error}}. \notag
       \end{aligned}$$ Using the algebraic-discretization equivalence relations , , and , the identity can be written in the equivalent form $$\begin{aligned}
         \label{BB3a}
         \J(\uh) = \underbrace{\J(\uhk)+ \etaAk}_{\mathrm{approximation}}
         + \underbrace{\ah(\uh - \uhk, \zh - \zhk )}_{\mathrm{error}}.
       \end{aligned}$$ Similarly, we can derive the dual relation $$\begin{aligned}
         \label{BB4}
         \bcT \mA^{-1} \bb
         = \underbrace{\bykT\bb + \bskT \bxk}_{\mathrm{approximation}}
         + \underbrace{\bskT\mA^{-1} \brk}_{\mathrm{error}},
       \end{aligned}$$ which is equivalent to $$\begin{aligned}
         \label{BB4a}
         \ell_h(\zh) =\underbrace{ \ell_h(\zhk) + \etaDAk}_{\mathrm{approximation}}
         + \underbrace{\ah(\uh - \uhk, \zh - \zhk )}_{\mathrm{error}}.
       \end{aligned}$$

As mentioned at the beginning of this section, both evaluations and can be used for any iterative solver generating approximations $\bxk$ and $\byk$, and residuals $\brk$ and $\bsk$, $k=1,2,\dots$. If the norms of the residual vectors $\brk$ and $\bsk$ tend to be small, one can expect that $$\begin{aligned}
  \label{BB5}
  | \bskT\mA^{-1} \brk | \ll |\bcT \mA^{-1} \brk| \qquad \mbox{and}\qquad
  |\bskT\mA^{-1} \brk| \ll |  \bskT\mA^{-1} \bb|.\end{aligned}$$ In other words, one can expect that the approximation is more accurate than since the corresponding error terms tend to be smaller.

Approximating the quantity of interest using BiCG iterates {#sec:bicg}
----------------------------------------------------------

As mentioned in Introduction, our aim is to solve systems by an iterative method, which allows to solve the primal and dual problems simultaneously. Since we intend to solve large and sparse systems, we need to pick a method with low memory requirements. Then, a natural choice is to use the [*preconditioned bi-conjugate gradient*]{} (BiCG) method (Algorithm \[alg:BiCG\]) introduced in [@Fl1976]; see also [@B:BaBeCh1994; @StTi2011]. Let $\mP$ be a suitable preconditioner (its choice depends on the particular discretization scheme). BiCG is a short-term recurrence Krylov subspace method which generates (if no breakdown occurs) approximations $\bxk \in \bx_0 + \mP^{-1}\mathcal{K}_k (\mA \mP^{-1},\br_0)$ and $\byk
\in \by_0 + \mP^{-\T}\mathcal{K}_k (\mA^{\T}\mP^{-\T},\bs_0)$ such that $$\label{eqn:bicg}
    \brk \perp \mP^{-\T}\mathcal{K}_k (\mA^{\T}\mP^{-\T},\bs_0),\qquad \bsk \perp \mP^{-1}\mathcal{K}_k (\mA\mP^{-1},\br_0),$$ where $\mathcal{K}_k (\mA,\br_0)$ denotes the $k$th Krylov subspace generated by $\mA$ and $\br_0$. The determining conditions imply that $$\begin{aligned}
    \label{orthog}
    (\byk-\by_0)^{\T} \brk = 0, \qquad 
    \bskT (\bxk-\bx_0) = 0.
  \end{aligned}$$ Note that during the BiCG finite precision computations, the bi-orthogonality conditions are usually not satisfied, and this fact cannot be ignored in our considerations. However, not all properties of BiCG vectors are lost in finite precision arithmetic. Because of the choice of coefficients $\alpha_{k}$ and $\beta_{k+1}$ (cf. Algorithm \[alg:BiCG\]), the bi-orthogonality of two consecutive vectors (local bi-orthogonality) is usually well preserved. This fact can be exploited to derive a more efficient way of approximating the quantities of interrest; for more details, see [@StTi2011].

Algorithm \[alg:BiCG\] shows the BiCG algorithm which generates sequences of primal and dual approximations $\{\bxk\}$ and $\{\byk\}$, respectively. Note that $\tilde{\br}_{k+1}=\mP^{-1} \br_{k+1}$ and $\tilde{\bs}_{k+1}=\mP^{-\T} \bs_{k+1}$ is equivalent to solving the system $$\begin{aligned}
  \label{mP}
  \mP\:\tilde{\br}_{k+1}=\br_{k+1}\quad\mbox{and}\quad \mP^{\T}\:\tilde{\bs}_{k+1}=\bs_{k+1},\end{aligned}$$ respectively. At lines 2 and 14 of Algorithm \[alg:BiCG\] we compute an additional sequence $\{\bxi_k^B\}$. The meaning of quantities $\bxi_k^B$ is explained in the text below. The algorithm has to be furnished by a suitable stopping criterion, which is discussed in Section \[sec:estim\].

**input** $\mA$, $\bx_{0}$, $\by_{0}$, $\mP$

$\bxi_{0}^{B}=0$

$\br_{0}=\bb-\mA\bx_{0},$ $\bs_{0}=\bc-\mA^{\T}\by_{0},$

$\bp_{0}=\mP^{-1} \br_{0}$

$\bq_{0}=\mP^{-\T} \bs_{0}$

$\tilde{\br}_{0}=\bp_{0}$

$\alpha_{k}=\frac{\bs_{k}^{\T}\tilde{\br}_{k}}{\bq_{k}^{\T}\mA\bp_{k}}$

$\bx_{k+1}=\bx_{k}+\alpha_{k}\bp_{k}$, $\by_{k+1}=\by_{k}+\alpha_k\bq_{k}$

$\br_{k+1}=\br_{k}-\alpha_{k}\mA\bp_{k}$, $\bs_{k+1}=\bs_{k}-\alpha_k\mA^{\T}\bq_{k}$

$\tilde{\br}_{k+1}=\mP^{-1} \br_{k+1}$, $\tilde{\bs}_{k+1}=\mP^{-\T} \bs_{k+1}$

$\beta_{k+1}=\frac{\bs_{k+1}^{\T}\tilde{\br}_{k+1}}{\bs_{k}^{\T}\tilde{\br}_{k}}$

$\bp_{k+1}=\tilde{\br}_{k+1}+\beta_{k+1}\bp_{k}$, $\bq_{k+1}=\tilde{\bs}_{k+1}+{\beta}_{k+1}\bq_{k}$

$\bxi_{k+1}^{B}=\bxi_{k}^{B}+\alpha_{k}\,\bs_{k}^{\T}\tilde{\br}_{k}$

To approximate $\bcT \mA^{-1} \bb$ in BiCG, we can use technique from [@StTi2011], which is based on the local bi-orthogonality of the BiCG vectors.

First, we present some manipulations. Let $\bx_0\in\IR^\Nh$ and $\by_0\in\IR^{\Nh}$ be the initial approximations of $\bx$ and $\by$, respectively. Using , , with $k=0$, and denoting $$\begin{aligned}
  \label{xis}
  \bxiP:= \bcT \bx_0 +   \by_0^{\T} \br_0, \qquad
  \bxiD:= \by_0^{\T} \bb +   \bx_0^{\T} \bs_0, \qquad 
  \bxiB:=  \bs_0 ^{\T} \mA^{-1} \br_0,\end{aligned}$$ we obtain the identities $$\begin{aligned}
  \label{BA4}
    \bcT\mA^{-1}\bb = \bxiP + \bxiB = \bxiD + \bxiB\qquad \Rightarrow \quad \bxiP = \bxiD.\end{aligned}$$ Based on technique from [@StTi2011] for estimating $\bxiB$ using $\bxiB_k$, we are now ready to introduce the next variant of approximating the quantity of interest.

-   In [@StTi2011 (3.13)] , it has been shown that $$\begin{aligned}
        \label{xi5}
        \bcT\mA^{-1}\bb  =  \underbrace{\bxiP + \bxi_{k}^{B}}_{\mathrm{approximation}}
        + \underbrace{\bskT \mA^{-1} \brk}_{\mathrm{error}},
        \quad k=0,1,\dots,
      \end{aligned}$$ with $$\begin{aligned}
        \label{xi3}
         \bxi_{k}^{B}:=\sum_{n=0}^{k-1} \alpha_{n}\,\bs_{n}^{\T}\tilde{\br}_{n},
      \end{aligned}$$ where scalars $\alpha_n$ and vectors $\bs_{n}$, $\tilde{\br}_{n}$, $n=0,1,\dots$ are defined by Algorithm \[alg:BiCG\]. The authors of [@StTi2011] derive this formula using the assumptions and the local bi-orthogonality conditions only. Therefore, if the consecutive vectors are almost bi-orthogonal and if the recursively computed residuals approximately agree with the true residuals during finite precision computations, then the identity holds (up to some small inaccuracy), and can be used for approximating the quantity of interest. The identity is mathematically equivalent to the identities and , and represents the third possibility of approximating the quantity of interest. The advantage of using is that we do not have to compute additional scalar products. The quantity $\bxi_{k}^{B}$ can be computed in BiCG almost for free, since the scalar products $\bs_{n}^{\T}\tilde{\br}_{n}$ are used in BiCG to compute the coefficients $\alpha_n$ and $\beta_{n+1}$.

Note that in exact arithmetic, the orthogonality relations hold. Then, comparing the identities , , and , we obtain $${\bcT\bxk + \bykT \brk}= 
{\bcT\bxk + (\byk-\by_0)^{\T} \brk+ \by_0^{\T} \brk } =  
{\bcT\bxk + \by_0^{\T} \brk }$$ so that $$\begin{aligned}
  \label{xi6}
  \underbrace{\bxiP + \bxi_{k}^{B}}_{\PPc}
  = \underbrace{\bcT\bxk + \bykT \brk}_{\PPb} = 
  \underbrace{\bcT\bxk}_{\PPa} + \by_0^{\T} \brk = \J(\uhk) + \by_0^{\T} \brk.\end{aligned}$$ In particular, if $\by_0=0$, then all the evaluations – are identical for the BiCG method in the exact arithmetic. Similarly, from , , and using the orthogonality , we obtain the dual counterpart relation $$\begin{aligned}
  \label{xi7}
  \underbrace{\bxiD + \bxi_{k}^{B}}_{\PPc}
  = \underbrace{\bykT\bb + \bskT \bxk}_{\PPb}
  = \underbrace{\bykT\bb}_{\PPa} + \bskT \bx_0 = \ell_h(\zhk) + \bskT \bx_0.\end{aligned}$$

In finite precision arithmetic, the first equalities in and still hold up to some small inaccuracy. However, if the orthogonality conditions are not (approximately) satisfied during finite precision computations, then the second equalities in and do not (approximately) hold; for more details and examples see [@StTi2011]. Note that in our experiments in Section \[sec:numerF\], the orthogonality conditions are well preserved and, therefore, the evaluations – provide almost the same results.

Finally, let us mention that if $\bx_0=0$ and $\by_0=0$, the relations and imply $$\begin{aligned}
  \label{xi8}
  \J(\uhk) =\ell_h(\zhk),\end{aligned}$$ which together with imply . We recall that the last relation is valid only if $\uhk$ and $\zhk$ are obtained by the BiCG method in exact arithmetic.

Estimation of the algebraic error {#sec:estim}
---------------------------------

In Sections \[sec:eqi\] and \[sec:bicg\], we presented three possibilities of the evaluations of the quantity of interest – . In this paragraph, we discuss the estimates of the errors of these evaluations, i.e., estimation of quantities $\estAk$ and $\estDAk$ introduced in and in , respectively. We use the standard approach when $\nu>0$ additional algebraic solver steps are performed and the difference between $k$ and $k+\nu$ iterates is used for the estimate of the error.

-   To estimate the error $\bcT \mA^{-1} \brk$ of evaluation , we subtract the identity in iterations $k$ and $k+\nu$, and obtain $$\begin{aligned}
      \label{EE3}
        \bcT \mA^{-1} \brk
        =  \underbrace{\bcT(\bxknu-\bxk)}_{\mathrm{estimate}}
        + \,\bcT \mA^{-1} \brknu.
      \end{aligned}$$

    Similarly, the dual identity gives $$\begin{aligned}
        \bskT \mA^{-1} \bb
        =  \underbrace{(\byknu-\byk)^\T \bb}_{\mathrm{estimate}}
        + \,\bsknuT \mA^{-1} \bb.
      \end{aligned}$$

-   Subtracting the identities in iterations $k$ and $k+\nu$ we can express the error $\bskT \mA^{-1} \brk$ of evaluation as $$\begin{aligned}
         \label{EE6}
        \bskT \mA^{-1} \brk
        =  \underbrace{\bcT(\bxknu-\bxk)+\byknuT\brknu-\bykT\brk}_{\mathrm{estimate}}
        + \,\bsknuT \mA^{-1} \brknu.
      \end{aligned}$$ Analogously, using we obtain $$\begin{aligned}
        \bskT \mA^{-1} \brk
        =  \underbrace{(\byknu-\byk)^\T \bb+\bsknuT\bxknu-\bskT\bxk}_{\mathrm{estimate}}
        + \,{\bsknuT \mA^{-1} \brknu}.
      \end{aligned}$$

-   Finally, considering the identity in iterations $k$ and $k+\nu$ we can express the error $\bskT \mA^{-1} \brk$ of evaluation as $$\begin{aligned}
        \label{EE16}
         \bskT \mA^{-1} \brk
       = \underbrace{\bxiBknu - \bxiBk }_{\mathrm{estimate}}
        +\, \bsknuT \mA^{-1} \brknu.
      \end{aligned}$$

Similarly as in Section \[sec:eqi\], starting with $\bx_0=0$ and $\by_0=0$, all errors of evaluations – as well as their estimates – are identical for the BiCG method in exact arithmetic. However, in finite precision arithmetic when the orthogonality is violated, the error of the evaluation and its estimate can substantially differ from the errors of the evaluations – and their estimates –.

Numerical experiments on fixed meshes {#sec:numerF}
=====================================

In this section we present the first collections of numerical experiments where approximation space $\Vh$ is fixed. The aim is to demonstrate the accuracy of the approximation of the quantity of interest $\J(\uh) = \bcT \mA^{-1} \bb = \ell_h(\zh)$ from Sections \[sec:eqi\] and \[sec:bicg\], and the estimates of the errors of these these approximations from Section \[sec:estim\].

We consider a second order elliptic problem which is discretized by the symmetric interior penalty Galerkin (SIPG) method using a piecewise polynomial but discontinuous approximation. SIPG method guarantees the the primal as well as dual consistencies and , respectively. For the definitions of the forms $\ah$ and $\ell_h$, we refer to [@ESCO-18], the detailed analysis can be found, e.g., in [@DGM-book]. All numerical examples presented in this paper were carried out using the in-house code ADGFEM [@ADGFEM] written in gfortran in double precision with processor i7-2620M CPU 2.70GHz (Ubuntu 16.04).

Elliptic problem on a “cross” domain {#sec:cross}
------------------------------------

We consider the example from [@Ainsworth2012Guaranteed Example 2] $$\begin{aligned}
  \label{eq:CR}
 - \Delta u  = 1 \quad \text{ in } \Om,  \qquad 
  u  = 0  \quad \text{ on } \partial\Om,\end{aligned}$$ where the “cross” domain $\Om = (-2,2) \times (-1,1) \cup (-1,1) \times (-2,2)$ and $\Delta$ denotes the Laplace operator. The target functional is defined as the mean value of the solution over the square $\Om_J = [1.2,1.4]\times [0.2,0.4],$ i.e. $J(u) = \frac{1}{|\Om_J|}\intX{\Om}{}{ j_\Om(x) u(x) }{x},$ where $j_\Om$ is the characteristic function of the square $\Om_J$, see Figure \[fig:CR\_meshes\], left. The exact value of $J(u)$ is unknown but we use the reference value $0.407617863684$, which was computed in [@Ainsworth2012Guaranteed] on an adaptively refined mesh with more than 15 million triangles.

The presence of interior obtuse angles of $\Om$ gives the singularities of the weak solution of . We carried out the computations on two triangular meshes, the first one is (quasi-)uniform having 3742 triangles and the second one is adaptively refined in the vicinity of interior angles and it has 4000 triangles, see Figure \[fig:CR\_meshes\]. For both meshes we used the SIPG method with $P_2$ and $P_4$ polynomial approximations

![Cross domain, computational domain $\Om$ an initial mesh and the domain of interest $\Om_J$ (small red square) (left), the finer uniform mesh (center) and the adaptively refined mesh (right).[]{data-label="fig:CR_meshes"}](Figs/CR_U0mesh.pdf "fig:"){height="33.00000%"} ![Cross domain, computational domain $\Om$ an initial mesh and the domain of interest $\Om_J$ (small red square) (left), the finer uniform mesh (center) and the adaptively refined mesh (right).[]{data-label="fig:CR_meshes"}](Figs/CR_U1mesh.pdf "fig:"){height="33.00000%"} ![Cross domain, computational domain $\Om$ an initial mesh and the domain of interest $\Om_J$ (small red square) (left), the finer uniform mesh (center) and the adaptively refined mesh (right).[]{data-label="fig:CR_meshes"}](Figs/CR_A2mesh.pdf "fig:"){height="33.00000%"}

Approximation of the quantity of interest {#sec:crossR}
-----------------------------------------

For each of the four corresponding discrete problems (uniform/adapted mesh and $P_2$/$P_4$ approximations) we carried out the solution of the corresponding algebraic systems by the BiCG method from Algorithm \[alg:BiCG\]. Since the method has tendency to stagnate after some number of iterations, we restarted the computations once after 400 BiCG iterations. Table \[tab:CR\] shows the limit values of $\J(\uh) = \bcT\mA^{-1}\bb$.

  -------------- --------- --------------- --------------------
   mesh uniform   $P_2  $   $\DoF=$ 22452   0.4071783143507507
   mesh uniform   $P_4  $   $\DoF=$ 56130   0.4075262691478035
   mesh adapted   $P_2  $   $\DoF=$ 24000   0.4076152998044911
   mesh adapted   $P_4  $   $\DoF=$ 60000   0.4076169203362077
  -------------- --------- --------------- --------------------

  : Cross domain, computation on fixed meshes, the limit values of the quanitity of interest $\J(\uh) = \bcT\mA^{-1}\bb$ after three restarts on two meshes (uniform/adapted) and $P_2$ and $P_4$ approximations.[]{data-label="tab:CR"}

Figures \[fig:BiCG\_convU1P2\] – \[fig:BiCG\_convA2P4\] show the results obtained within the first 400 BiCG iterations[^2], namely

-   the convergence of the errors three types of the approximation of the quantity of interest $\J(\uh) = \bcT\mA^{-1}\bb$, namely

    -   : $ | \bcT\bxk - \J(\uh) | $ following from ,

    -   : $ | \bcT\bxk + \bykT \brk -  \J(\uh) | $ following from ,

    -   : $ | \bxiP + \bxi_{k}^{B} -  \J(\uh) | $ following from ,

    where $\J(\uh)$ is the limit value from Table \[tab:CR\] and $|\cdot|$ is the absolute value,

-   the convergence of the values of three different estimations of the error of the approximation of the quanitity of interest, namely

    -   : $ | \bcT\delta \bxk |$ following from ,

    -   : $ | \bcT \delta \bxk +\delta(\bykT\brk) |$ following from ,

    -   : $ | \delta \bxiBk|$ following from ,

    where symbol $\delta$ denotes the difference between $k$-th and $(k+\nu)$-th iterations, i.e, $\delta \bxk= \bxknu - \bxk$, $\delta(\bykT\brk)= \byknuT\brknu-\bykT\brk$, $\delta \bxiBk = \bxiBknu - \bxiBk$, in the experiments, we put uniquely $\nu=10$,

-   the quantity

    -   $ \frac{| (\byk-\by_0)^\T \brk  |}{\|\byk-\by_0 \|\, \| \brk\|} $ measuring the lost of the orthogonality, cf. .

We observe the following.

1.  The errors of all approximations of the quantity of interest (solid lines) decrease until they reach some level of accuracy. After reaching this level they stagnate. This is in agreement with the results in [@Gr1997] where its is showen that the residual norms $\| \brk \|$ of the Krylov subspace methods like BiCG, reach (if they converge) the level of accuracy close to $\varepsilon \|A\|  \Theta $ where $\varepsilon$ is the machine precision and $\Theta = \max\{\|\bxk\|\}$ (unpreconditioned case).

2.  All approximation of the quantity of interest – (solid lines) have almost identical convergence as long as the lost of orthogonality is small (dotted-dashed line). When the lost of orthogonality starts to play more important role the values of the approximation of the quantity of interest slightly differ.

3.  Similarly, all the estimations of the error of approximations – (dashed lines) have almost identical convergence as long as the lost of orthogonality is small. If the evaluations of the quantity of interests start to stagnate the estimates – underestimate the error. This process is the slowest for but it is presented too.

From these observations, we conclude that any error estimate – can be used for the stopping criterion. They underestimate the error only when the computational process stagnates and then it make no sense to proceed with next iterative steps. Based on the argumentation in Section \[sec:estim\], criteria and are less sensitive to the lost of orthogonality. Finally, criterion is cheaper to evaluate than but, in our case, it does not play any essential role in the whole computational process.

![Convergence of the algebraic error and its estimates for the BiCG solver, case uniform mesh, $P_2$ approximation, the errors of the approximation of the quantity of interest , , (solid lines), their approximation , , (dashed lines) and the lost of the orthogonality (dotted-dashed line).[]{data-label="fig:BiCG_convU1P2"}](Figs/BiCG_convU1P2_P.pdf){width="100.00000%"}

![Convergence of the algebraic error and its estimates for the BiCG solver, case uniform mesh, $P_4$ approximation, the errors of the approximation of the quantity of interest , , (solid lines), their approximation , , (dashed lines) and the lost of the orthogonality (dotted-dashed line).[]{data-label="fig:BiCG_convU1P4"}](Figs/BiCG_convU1P4_P.pdf){width="100.00000%"}

![Convergence of the algebraic error and its estimates for the BiCG solver, case adapted mesh, $P_2$ approximation, the errors of the approximation of the quantity of interest , , (solid lines), their approximation , , (dashed lines) and the lost of the orthogonality (dotted-dashed line).[]{data-label="fig:BiCG_convA2P2"}](Figs/BiCG_convA2P2_P.pdf){width="100.00000%"}

![Convergence of the algebraic error and its estimates for the BiCG solver, case adapted mesh, $P_4$ approximation, the errors of the approximation of the quantity of interest , , (solid lines), their approximation , , (dashed lines) and the lost of the orthogonality (dotted-dashed line).[]{data-label="fig:BiCG_convA2P4"}](Figs/BiCG_convA2P4_P.pdf){width="100.00000%"}

Adaptive mesh refinement and algebraic stopping criteria {#sec:adapt}
========================================================

The goal of the numerical solution of is to obtain a numerical approximation $\tuh$ such that (cf. ) $$\begin{aligned}
  \label{goal}
  |\J(\u) - \J(\tuh)| \approx   \etaSk + \etaAk \leq \TOL,\end{aligned}$$ where $\TOL>0$ is the given tolerance. The adaptive mesh refinement allows to reduce the computation costs necessary to achieve .

Mesh adaptation algorithm
-------------------------

The idea of the mesh adaptive algorithm is to start on an initial coarse mesh $\Thz$ (dimension of the corresponding space $\Vhz$ is small). Then for $m=0,1,\dots$, we discretize and solve both primal and dual problem on $\Vhm$ and estimate the error of the quantity of interest. If the estimate does not fulfil then we adapt the mesh and create a new one $\ThmP$ and proceed with the computation. Algorithm \[alg:AMAhp\] shows the abstract form of the adaptive mesh refinement algorithm. It means that we obtain the approximations of the primal and dual solutions $\uhmk\in\Vhm$ and $\zhmk\in\Vhm$, where subscript $m$-th corresponds to the level of mesh adaptation and subscript $k$-th corresponds to the algebraic iteration.

[ let $\Thz$ be the initial (coarse) mesh let $\Vhz$ be the corresponding finite element space set the algebraic problems corresponding to the discretization of and , respectively, on $\Vhm$, i.e., $$\begin{aligned}
  \label{AAm}
  \mA_m \bx_{m} = \bb_{m}\qquad\mbox{and}\qquad \mA_{m}^{\T}\by_m=\bc_m,\end{aligned}$$ where $\bb_{m},\bc_m \in\R^{N_m}$, $\mA_m\in\R^{N_m\times N_m}$, $\dim\Vhm=N_m$, apply an iterative algebraic solver evaluate $\bx_{m,k}$ and $\by_{m,k}$ approximations of EXIT define $\uhmk,\zhmk\in\Vhm$ corresponding to $\bx_{m,k}$ and $\by_{m,k}$ (the outputs of Algorithm \[alg:BiCG\]) define higher-order reconstructions $\huhm=\krr(\uhmk)$ and $\hzhm=\krr(\zhmk)$ using employing –, evaluate the goal oriented error estimates $$\begin{aligned}
    \label{etasm}
    \etamk := \etaSmk + \etaAmk, \qquad \etaDmk := \etaDSmk + \etaDAmk
  \end{aligned}$$ where $$\begin{aligned}
    \label{etasm1}
    &\etaSmk := \res{\uhmk}{\hzhm -  \zhmk},  &&\etaAmk := \res{\uhmk}{\zhmk}, \\
    &\etaDSmk := \resD{\zhmk}{\huhm -  \uhmk},  &&\etaDAmk := \resD{\zhmk}{\uhmk}, \notag
  \end{aligned}$$ , based on a localization of $\etamk$, adapt the mesh, i.e., modify the size and/or the shape of elements and possibly also the polynomial approximation degrees, the new mesh is $\ThmP$ and the corresponding space $\VhmP$, ]{}

Step 18 of Algorithm \[alg:AMAhp\] (mesh adaptation) depends on the used discretization method and the refinement technique. In order to demonstrate the robustness of the presented stopping criteria, we use the anisotropic $hp$-mesh adaptation approach from [@ESCO-18], which generate anisotropic meshes (consisting of possibly thin and long triangular elements) and varying polynomial approximation degree.

The crucial aspect is the algebraic stopping criterion in step 8 of Algorithm \[alg:AMAhp\] which is discussed in the next section.

Standard stopping criteria for the solution of algebraic systems {#sec:stop1}
----------------------------------------------------------------

We focus on the algebraic stopping criterion in Step 8 of Algorithm \[alg:AMAhp\]. Obviously, too strong criterion leads to many algebraic iterations without a gain of accuracy. On the other hand, the weak criterion leads to an under-solving of which affect the mesh adaptation process. Typically too many mesh elements are generated.

Often the residual stopping criteria for are used, i.e., $$\begin{aligned}
  \label{RES}
  \| \br_{m,k} \| := \| \bb_m - \mA_m \bx_{m,k} \| \leq \TOLA, \qquad
  \| \bs_{m,k} \| := \| \bc_m - \mA_m^{\T} \by_{m,k} \| \leq \TOLA, \end{aligned}$$ or their preconditioned variant $$\begin{aligned}
  \label{RESP}
  \| \mP^{-1} (\bb_m - \mA_m \bx_{m,k}) \| \leq \TOLA, \qquad
  \| \mP^{-\T} (\bc_m - \mA_m^{\T} \by_{m,k}) \| \leq \TOLA, \end{aligned}$$ where $\mP$ represents a suitable preconditioner, cf. . This criterion is easy to evaluate but the choice of the suitable tolerance $\TOLA>0$ is difficult since this criterion has no relation to the discretization error.

This drawback was eliminated in [@Meidner2009Goal] by the following stopping criterion $$\begin{aligned}
  \label{aDWR}
  \etaAmk \leq \cA \etaSmk,\qquad
  \etaDAmk \leq \cA \etaDSmk,\end{aligned}$$ where the primal/dual estimates of the algebraic/discretization errors are defined by and $\cA\in(0,1)$ is a suitable constant. This means that Steps 12-14 of Algorithm \[alg:AMAhp\] are moved inside the inner loop (after Step 7). We call this criterion as [*algebraic goal-oriented*]{} stopping criterion.

The conditions allow to control the size of the algebraic error. However, a strong drawback of are to computational costs. Whereas the evaluation of $ \etaAmk $ and $ \etaDAmk$ is cheap, the computation of $ \etaSmk $ and $ \etaDSmk$ is much more expensive namely due to the necessity to perform the higher-order reconstructions $\huhm$ and $\hzhm$ (Step 13). The computational costs can be reduced by testing only after some number of iterations, e.g., after 1-3 restarts of iterative solver. In [@Meidner2009Goal], the first condition of was tested after the performing of one cycle of multigrid method.

New stopping criteria for the solution of algebraic systems {#sec:stop2}
-----------------------------------------------------------

In Section \[sec:alter\], we introduced two alternative formula for the decomposition of the computational errors into the discretization and algebraic parts, namely using and , we have $$\begin{aligned}
{3}
    \label{ST1}
    \J(\u - \uhmk) &  = \res{\uhmk}{\z -  \zhmk} + \res{\uhmk}{\zhmk} &&=\errSmk +\errAmk,
    \\
    \J(\u - \uhmk) &  =  \J(\u - \uhm) + \J(\uhm - \uhmk)&&= \estSmk +\estAmk. \notag
  \end{aligned}$$ Similarly, and imply the dual counterpart $$\begin{aligned}
{3}
    \label{ST2}
    \ell_h(\z - \zhmk) & = \resD{\zhmk}{\u -  \uhmk} + \resD{\zhmk}{\uhmk} &&= \errDSmk + \errDAmk, \\
    \ell_h(\z - \zhmk) &= \ell_h(\z - \zhm) + \ell_h(\zhm - \zhmk)   &&= \estDSmk + \estDAmk.
    \notag
  \end{aligned}$$ The discussion presented therein shows that both quantities $\errAmk\not=\estAmk$ corresponds to the algebraic error of the primal problem and similarly, $\errDAmk\not=\estDAmk$ corresponds to the algebraic error of the dual one.

In Section \[sec:estim\], we presented several techniques estimating the quantities $\estAmk$ and $\estDAmk$ by techniques – . On the other hand, quantities $\errAmk$ and $\errDAmk$ are available during the BiCG iterative method since $$\begin{aligned}
  \label{ST3}
  \errAmk &= \etaAmk = \res{\uhmk}{\zhmk} = \bykT \brk, \\
  \errDAmk &= \etaDAmk = \resD{\zhmk}{\uhmk} = \bskT \bxk, \notag\end{aligned}$$ cf. and .

Figure \[fig:adaptBiCG\] shows a typical dependence of these quantities w.r.t. the number of BiCG iterations during mesh adaptation process by Algorithm \[alg:AMAhp\] for $m=0,\dots, 5$. In the top figure, we plot the total error $J(\u-\uhmk)$, the exact algebraic error $\estAmk=J(\uhm-\uhmk)$, the algebraic error $\errAmk=\bskT \bxk$ and the estimates and of $\estAmk$ from Section \[sec:estim\]. The bottom figure shows the dual counterparts. Estimate gives the same graphs as so we do not show it.

![Typical example of the adaptive computations for mesh adaptation steps $m=0,\dots,5$, dependence of primal (top) and dual (bottom) quantities w.r.t. the number of BiCG steps.[]{data-label="fig:adaptBiCG"}](Figs/BiCG_adapt.pdf){width="\textwidth"}

![Typical example of the adaptive computations for mesh adaptation steps $m=0,\dots,5$, dependence of primal (top) and dual (bottom) quantities w.r.t. the number of BiCG steps.[]{data-label="fig:adaptBiCG"}](Figs/BiCG_adaptD.pdf){width="\textwidth"}

We observe that estimate approximates $\estAmk$ similarly as for the moderate values of accuracy but much better on the level close to the machine accuracy. This is in agreement with theoretical considerations in Section \[sec:estim\], but this effect is not observed on the initial mesh ($m=0$) when we do not have a good initial approximation. Further, both algebraic error representations converge but with different speed. Whereas on the initial mesh we have $\eta_{\mathrm{A},0,k} \ll  \zeta_{\mathrm{A},0,k}$ for $k\le 100$, starting from $m=1$ we observe $\eta_{\mathrm{A},m,k} \gg  \zeta_{\mathrm{A},m,k}$. Similar behaviour is observed for the dual quantities.

Based on this observation, we consider both pairs of quantities $(\etaAmk,  \estAmk)$ and $(\etaDAmk,  \estDAmk)$ for the definition of the stopping criterion, namely, we evaluate the quantities $$\begin{aligned}
  \label{estMK1}
  \ESTAmk := & | \delta \bxiBk| + | \etaAmk| = |\bxiBknu - \bxiBk| + | \by_{k,m}^{\T} \br_{k,m}|, \\
  \ESTDAmk := & | \delta \bxiBk| + | \etaDAmk| = |\bxiBknu - \bxiBk| + | \bs_{k,m}^{\T} \bx_{k,m}|, \notag\end{aligned}$$ corresponding to the $k$-iteration of the BiCG Algorithm \[alg:BiCG\] on $m$-level of mesh adaptation. Let us recall that quantities $\ESTAmk$ and $\ESTDAmk$ are available during the whole BiCG iterative process with negligible computational costs.

Now we define a new algebraic stopping criterion (called hereafter $\sigma$-stopping criterion) for Algorithm \[alg:AMAhp\] (Step 8) by $$\begin{aligned}
  \label{xiB}
  \ESTAmk  \leq \cA \TOL\qquad \mbox{and} \qquad
  \ESTDAmk \leq \cA \TOL, \end{aligned}$$ where $ \TOL>0$ is the (global) tolerance from and $\cA\in(0,1)$ is a suitable constant. In contrary to the stopping criterion , the proposed one does not take into account the (estimate of) the discretization error and consequently the strong advantage of is that its evaluation is very fast.

We can summarized the idea of the original stopping criterion as follows: “[*algebraic system is solved as long the algebraic error is $\cA$-times smaller than the discretization error*]{}”. On the other hand, the idea of the new stopping criterion is the following: “[*algebraic system is solved as long the algebraic error is $\cA$-times smaller than the prescribed tolerance for the discretization error*]{}”.

The stopping criterion does not look to much efficient since for starting levels of mesh adaptation, when $\etaSmk \gg\TOL$ and $\etaDSmk \gg\TOL$, the iterative solver is stopped when the algebraic error is much lower than the discretization one. It is true but on the other hand performing of additional several teens of BiCG iterations is typically faster than the evaluation of $\etaSmk$ and $\etaDSmk$. Moreover, approximate solutions computed using BiCG on starting levels serve as good initial approximations for BiCG on higher levels, which improve then the solving process significantly.

Numerical experiments on adaptively refined meshes {#sec:numerA}
==================================================

We demonstrate the computational performance of the stopping criteria from Sections \[sec:stop1\] and \[sec:stop2\]. We consider two numerical examples, the first one is the elliptic problem on the “cross” domain described in Section \[sec:cross\] and the second one is a convection-dominated problem having some anisotropic features, it is defined in Section \[sec:carpio\]. Both problems are discretized again by the SIPG method and the meshes are adapted by the technique from [@ESCO-18]. Each mesh $\Thm$ is generated from the computed primal and dual solutions on the previous one. These solutions are interpolated to the actual mesh, hence we have relatively very good initial approximations $\bx_{m,0}$ and $\by_{m,0}$ for the solution of .

Convection-dominated problem {#sec:carpio}
----------------------------

The second example is taken from [@FormaggiaPerotto04], see also [@ESCO-18; @CarpioPrietoBermejo_SISC13]. We solve the convection-diffusion equation $$\begin{aligned}
  \label{carpio}
  - \ve \Delta u + \nabla \cdot (\bkb\, u) &= 0  \qquad
  \text{in } \Om:=[0, 4]\times [0, 4] \setminus [0, 2] \times [0, 2],\end{aligned}$$ where $\ve = 10^{-3}$, the convection field $\bkb = (x_2,-x_1)$ and $\nabla \cdot$ is the divergence operator. We prescribe the Dirichlet and Neumann boundary conditions $$\begin{aligned}
  \label{carpio:BC}
  u=1 &\quad  \mbox{ on } \{x_1 = 0\}, \\
  \nabla u \cdot \bkn = 0 & \quad \mbox{ on } \Gamma_1:=\{\gom; x_1 = 4\} \cup \Gamma_2:=\{\gom; x_2 = 0\},  \notag \\
  u=0 & \quad \mbox{ elsewhere}. \notag\end{aligned}$$ The solution $u$ exhibits boundary layers as well as two circular-shaped internal layers. We consider the functional $  J(u) = \int_{\Gamma_1} \bkb\cdot\bkn\, u\dS$ with the reference value $J(u) = 0.07408122\pm 10^{-8}$.

Comparison of the stopping criteria
-----------------------------------

In order to demonstrate the robustness of the proposed stopping criteria, we solve the elliptic problem and the convection-diffusion problem by the mesh adaptation Algorithm \[alg:AMAhp\]. For both problems, we employ the anisotropic $h$-mesh adaptation using $P_3$ approximation and anisotropic $hp$-mesh adaptation, for details we refer to [@ESCO-18]. In order to observe the effect with not sufficiently resolved algebraic systems , we stop the adaptation process when $$\begin{aligned}
  \label{STOP}
  \frac12\left(\etaSmk + \etaDSmk\right)\leq \TOL,\end{aligned}$$ where $\etaSmk$ and $\etaDSmk$ are given by and we put $\TOL=10^{-8}$ for the elliptic problem and $\TOL=10^{-10}$ for the convection-diffusion problem .

We test Algorithm \[alg:AMAhp\] with the following iterative solvers and the stopping criteria.

1.  GMRES with preconditioned residual stopping criterion with tolerances $\TOLA=10^{-9}$, $10^{-6}$ and $10^{-3}$,

2.  BiCG with preconditioned residual stopping criterion with tolerances $\TOLA=10^{-9}$, $10^{-6}$ and $10^{-3}$,

3.  GMRES with algebraic goal-oriented stopping criterion with tolerances $\cA=10^{-2}$, $10^{-1}$ and $10^{0}$,

4.  BiCG with algebraic goal-oriented stopping criterion with tolerances $\cA=10^{-2}$, $10^{-1}$ and $10^{0}$,

5.  BiCG with $\sigma$-stopping criterion with tolerances $\cA=10^{-2}$, $10^{-1}$ and $10^{0}$.

For all solvers we use the block-ILU(0) preconditioner which is suitable for discontinuous Galerkin method, cf. [@lin_algeb]. For GMRES, we solve first the primal problem and then the dual problem, for BiCG, both problems are solved at once, of course. Moreover, GMRES is restarted always after 45 iterations. In case (c), the stopping criterion is tested after one restart of GMRES for the primal as well as the dual problems. In case (d), the stopping criterion is tested after 100 BiCG iterations.

The results are presented in Tables \[tab:CRh\] – \[tab:JBhp\], where we show the number of degrees of freedom $\DoF$ on the last mesh, the total error $\J(\u-\uhmk)$, the final estimator of the discretization and algebraic errors $\etaSmk$ and $\etaAmk$, respectively and their sum $\etamk$. Moreover, these tables contains the total number of (GMRES or BiCG) iterations on all mesh levels $m=0,1,\dots$ (iters) and the total computational time in seconds. If $\etaAmk \ge \etaSmk$ then the line ends with the character ’!’.

We remind that the error estimator $\etamk$ does not give the upper bound, see Remark \[rem:guaranteed\]. There are small differences for each of four examples but we can state the following observations.

-   the use of the preconditioned residual stopping criterion for GMRES as well as BiCG with too large $\TOLA$ leads to the dominance of the algebraic error and consequently $J(u-\uhmk)$ and $\etamk$ are much larger then the tolerance $\TOL$ from and , cf. (a) and (b),

-   GMRES requires more iterations then BiCG since the primal and dual problems are solver separately, therefore BiCG is faster, compare (a) vs. (b) and (c) vs. (d),

-   the goal-oriented stopping criterion allows to control the algebraic (and therefore also the total) error as expected for GMRES and BiCG, see (c) and (d),

-   the $\sigma$-stopping criterion works efficiently, the algebraic error estimator $\etaAmk$ is smaller than the discretization one $\etaSmk$, the computational time is substantially reduced in comparison with the goal-oriented criterion since the higher-order reconstructions is performed less frequently, see the last paragraph in Section \[sec:stop1\].

Summary of the results and outlook {#sec:sum}
==================================

We developed an efficient technique for the numerical solution of primal and dual algebraic systems arising in the goal-oriented error estimation and mesh adaptation. Both algebraic systems are solved simultaneously by BiCG method which allows to control the algebraic error during the iterative process. The proposed $\sigma$-stopping criterion is cheap for the evaluation and significantly reduce the computational costs. Moreover, it guarantees that the algebraic error estimate bounded by the discretization one.

Further natural step is to extend this approach for the solution of nonlinear partial differential equations. The the dual problem has to be build on a linearization of the primal one. However, employing a Newton-like method for the solution of the discretized primal problem, an approximate solution of the dual problem is available at each Newton step and the technique developed in this paper can be employed. However, it is necessary to balance the linear algebraic errors, the non-linear algebraic errors and the discretization errors. This is the subject of the further work.

#### Acknowledgements

The authors are thankful to their colleagues from the Charles University, namely M.Kub[í]{}nov[á]{}, T. Gergelits and F. Roskovec for a fruitful discussion.

[^1]: This work was supported by grant No. 17-04150J of the Czech Science Foundation.

[^2]: After the restart, the machine accuracy is achieved in few steps and the error estimators give vanishing values, therefore we do not show them.
