Adaptive solution of linear systems of equations based on a posteriori error estimators

Abstract

In this paper, we discuss a new adaptive approach for iterative solution of sparse linear systems arising from partial differential equations (PDEs) with self-adjoint operators. The idea is to use the a posteriori estimated local distribution of the algebraic error in order to steer and guide the solve process in such way that the algebraic error is reduced more efficiently in the consecutive iterations. We first explain the motivation behind the proposed procedure and show that it can be equivalently formulated as constructing a special combination of preconditioner and initial guess for the original system. We present several numerical experiments in order to identify when the adaptive procedure can be of practical use.

Appendix: A posteriori error estimates

In this section, we describe briefly one of the ways how to estimate the local distribution of the algebraic error (6) using flux reconstructions following [8, 10, 17] and references therein. We start by introducing the basic techniques of these a posteriori error estimates, then we detail how to get a sharp computable upper bound of the algebraic error. Note that this section recalls existing techniques and results and we only adapt the a posteriori error estimate of [17] to our chosen model problem.

1.1 Basic a posteriori error estimates

Assuming that a Galerkin solution u_h is available, we start by bounding the energy norm of the error $\underline {u} - \underline {u}_{h}$ represented as

$$ \lVert{\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u} - \underline{u}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} = \sup_{\underline{v} \in V, \lVert{\nabla \underline{v}}\rVert=1} (\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u} - \underline{u}_{h}), \underline{\mathbf{K}}^{\frac12}\nabla \underline{v}). $$

(38)

Note that following (2), then

$$ \lVert{\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u} - \underline{u}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} = \sup_{\underline{v} \in V, \lVert{\nabla \underline{v}}\rVert=1}\{(\underline{f},\underline{v})-(\underline{\mathbf{K}}^{\frac12}\nabla \underline{u}_{h},\underline{\mathbf{K}}^{\frac12}\nabla \underline{v})\}. $$

(39)

The key ingredient of our estimate is a reconstructed flux 𝜃_h which is a piecewise polynomial function in the Raviart–Thomas–Nédélec subspace $\textbf {RTN}(\mathcal {T}_{h})$ of the infinite-dimensional space H(div,Ω) which is reconstructed in order to mimic the continuous flux ${\boldsymbol \theta } := -\underline {\mathbf {K}} \nabla \underline {u}$. In other words, 𝜃_h is reconstructed to satisfy

$$ \nabla \cdot {\boldsymbol \theta}_{h} = \underline{f}. $$

(40)

Recall from Section 2 that $\underline {f}$ is assumed to be piecewise constant with respect to the mesh $\mathcal {T}_{h}$. Now, we use the Green and the Cauchy–Schwarz inequality together with (39) and we follow [17, Section 4.1], to write

$$ \begin{array}{@{}rcl@{}} \lVert{\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u} - \underline{u}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} &=& \underset{\nabla \cdot {\boldsymbol \sigma} = \underline{f}}{\underset{{\boldsymbol \sigma} \in \textbf{H}(\text{div},{\Omega})}{\inf}} \underset{{\lVert{\nabla \underline{v}}\rVert=1}}{\underset{\underline{v} \in V}{\sup}}\{(\underline{f}-\nabla \cdot {\boldsymbol\sigma} ,\underline{v}) \\ & & \qquad \qquad \qquad -(\underline{\mathbf{K}}^{\frac12}\nabla \underline{u}_{h} + \underline{\mathbf{K}}^{-\frac12} {\boldsymbol \sigma},\underline{\mathbf{K}}^{\frac12}\nabla \underline{v})\} \\ &=& \underset{\nabla \cdot {\boldsymbol \sigma} = \underline{f}}{\underset{{\boldsymbol \sigma} \in \textbf{H}(\text{div},{\Omega})}{\inf}} \underset{\lVert{\nabla \underline{v}}\rVert=1}{\underset{\underline{v} \in V}{\sup}} \{-(\underline{\mathbf{K}}^{\frac12}\nabla \underline{u}_{h} + \underline{\mathbf{K}}^{-\frac12} {\boldsymbol \sigma},\underline{\mathbf{K}}^{\frac12}\nabla \underline{v})\} \\ &=& \underset{\nabla \cdot {\boldsymbol \sigma} = \underline{f}}{\underset{{\boldsymbol \sigma} \in \textbf{H}(\text{div},{\Omega})}{\inf}} \lVert{\underline{\mathbf{K}}^{\frac12}\nabla \underline{u}_{h} + \underline{\mathbf{K}}^{-\frac12} {\boldsymbol \sigma}}\rVert_{\mathrm{L}^{2}({\Omega})} \\ &\leq& \lVert{\underline{\mathbf{K}}^{\frac12}\nabla \underline{u}_{h} + \underline{\mathbf{K}}^{-\frac12} {\boldsymbol \theta}_{h}}\rVert_{\mathrm{L}^{2}({\Omega})}. \end{array} $$

This gives a guaranteed upper bound on the discretization error. Note that the obtained estimate relies only on the weak formulation and the reconstructed flux 𝜃_h. Finally, it is important to mention that the needed reconstructed flux 𝜃_h can be easily reconstructed for various discretization schemes like finite elements, nonconforming finite elements, discontinuous Galerkin, finite volumes, and mixed finite elements (see [8] for more details).

1.2 Upper bound on the algebraic error

In this section, we suppose that we use an iterative solver to obtain an approximate solution $\underline {u}^{(i)}_{h}$ of (3) after running i iterations. In order to estimate the algebraic error we first introduce a representation of the algebraic residual vector which will be a function $\underline {s}_{h}^{(i)}\in \mathrm {L}^{2}({\Omega })$ satisfying

$$ (\underline{s}_{h}^{(i)},\varphi_{j}) = \mathbf{r}^{(i)}_{j}, \qquad 1\leq j\leq n. $$

(41)

Details about the reconstruction of $\underline {s}_{h}^{(i)}$ with two different examples can be found in [17, Section 5.1]. Note that from (41) and the definition of the algebraic residual vector r⁽ⁱ⁾ one can write

$$ (\underline{s}_{h}^{(i)},\varphi_{j}) = (\underline{f},\varphi_{j}) -(\underline{\mathbf{K}}^{\frac12} \nabla \underline{u}^{(i)}_{h}, \underline{\mathbf{K}}^{\frac12} \nabla \varphi_{j}), \qquad 1\leq j\leq n. $$

(42)

Consequently,

$$ (\underline{s}_{h}^{(i)},\underline{v}_{h}) = (\underline{f},\underline{v}_{h}) -(\underline{\mathbf{K}}^{\frac12} \nabla \underline{u}^{(i)}_{h},\underline{\mathbf{K}}^{\frac12} \nabla \underline{v}_{h}). $$

(43)

Using (3), then (43) gives

$$ (\underline{s}_{h}^{(i)},\underline{v}_{h}) = (\underline{\mathbf{K}}^{\frac12} (\nabla \underline{u}_{h} - \nabla \underline{u}^{(i)}_{h}), \underline{\mathbf{K}}^{\frac12}\nabla \underline{v}_{h}). $$

(44)

This representation of the algebraic residual vector plays a key role in the estimation of the algebraic error. Actually, by applying the Cauchy–Schwarz inequality together with the Friedrichs inequality on (44), one gets

$$ \begin{array}{@{}rcl@{}} (\underline{\mathbf{K}}^{\frac12} (\nabla \underline{u}_{h} - \nabla \underline{u}^{(i)}_{h}) ,\underline{\mathbf{K}}^{\frac12}\nabla \underline{v}_{h}) &=& (\underline{s}_{h}^{(i)},\underline{v}_{h}) \\ &\leq& \|\underline{s}_{h}^{(i)}\|_{\mathrm{L}^{2}({\Omega})} \lVert{\underline{v}_{h}}\rVert_{\mathrm{L}^{2}({\Omega})}\\ &\leq& \|\underline{s}_{h}^{(i)}\|_{\mathrm{L}^{2}({\Omega})} \left( C_{\mathrm{F}} h_{\Omega} \lVert{\nabla\underline{v}_{h}}\rVert_{\mathrm{L}^{2}({\Omega})}\right)\\ &\leq& \|\underline{s}_{h}^{(i)}\|_{\mathrm{L}^{2}({\Omega})} \left( C_{\mathrm{F}} h_{\Omega}{\lambda}^{-\frac12}_{\underline{\mathbf{K}}} \lVert{\underline{\mathbf{K}}^{\frac12}\nabla\underline{v}_{h}}\rVert_{\mathrm{L}^{2}({\Omega})}\right), \end{array} $$

where 0 < C_F ≤ 1 is the constant from the Friedrichs inequality, h_Ω the diameter of the domain Ω, and ${\lambda }_{\underline {\mathbf {K}}}$ the smallest eigenvalue of $\underline {\mathbf {K}}$. First computable upper bound is then obtained for the algebraic error as

$$ \lVert{\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u}_{h} - \underline{u}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} \leq C_{\mathrm{F}} h_{\Omega}{\lambda}^{-\frac12}_{\underline{\mathbf{K}}} \|\underline{s}_{h}^{(i)}\|_{\mathrm{L}^{2}({\Omega})}. $$

(45)

However, this upper bound yields often a significant overestimation (see, e.g., [17, Sections 3.1 and 4.2]). An improvement of the upper bound (45) can be obtained by using flux reconstruction techniques and additional algebraic iterations. Following Section 6 and [17, Section 5.3], we consider a reconstructed flux ${\boldsymbol \theta }^{(i)}_{h} \in \textbf {RTN}(\mathcal {T}_{h})$ satisfying $\nabla \cdot {\boldsymbol \theta }^{(i)}_{h} = \underline {f} - \underline {s}_{h}^{(i)}$. Then, after ν > 0 additional iterations we similarly construct from the algebraic residual vector r^(i+ν) a representation $\underline {s}_{h}^{(i+\nu )}$, and another reconstructed flux ${\boldsymbol \theta }^{(i+\nu )}_{h} \in \textbf {RTN}(\mathcal {T}_{h})$ satisfying $\nabla \cdot {\boldsymbol \theta }^{(i+\nu )}_{h} = \underline {f} - \underline {s}_{h}^{(i+\nu )}$. With these different reconstructions, we have

$$ \underline{s}_{h}^{(i)} = \underline{s}_{h}^{(i+\nu)} + \nabla \cdot {\boldsymbol \theta}^{(i+\nu)}_{h} - \nabla \cdot {\boldsymbol \theta}^{(i)}_{h}, $$

and therefore (44) gives

$$ (\underline{\mathbf{K}}^{\frac12} (\nabla \underline{u}_{h} - \nabla \underline{u}^{(i)}_{h}) ,\underline{\mathbf{K}}^{\frac12}\nabla \underline{v}_{h}) = (\underline{\mathbf{K}}^{-\frac12} ({\boldsymbol \theta}^{(i+\nu)}_{h} - {\boldsymbol \theta}^{(i)}_{h}), \underline{\mathbf{K}}^{\frac12} \nabla\underline{v}_{h} ) +(\underline{s}_{h}^{(i+\nu)},\underline{v}_{h}) . $$

Consequently,

$$ \lVert{\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u}_{h} - \underline{u}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} \leq \lVert{\underline{\mathbf{K}}^{- \frac12} ({\boldsymbol \theta}^{(i+\nu)}_{h} - {\boldsymbol \theta}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} + C_{\mathrm{F}} h_{\Omega}{\lambda}^{-\frac12}_{\underline{\mathbf{K}}} \lVert{\underline{s}_{h}^{(i+\nu)}}\rVert_{\mathrm{L}^{2}({\Omega})} $$

(46)

The idea of using additional algebraic iterations is very useful in practice [9]. In fact, for a sufficiently large value of ν, one can assume that there exists γ > 0 such that

$$ C_{\mathrm{F}} h_{\Omega}{\lambda}^{-\frac12}_{\underline{\mathbf{K}}} \lVert{\underline{s}_{h}^{(i+\nu)}}\rVert_{\mathrm{L}^{2}({\Omega})} \leq \gamma \lVert{\underline{\mathbf{K}}^{- \frac12} ({\boldsymbol \theta}^{(i+\nu)}_{h} - {\boldsymbol \theta}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})}, $$

so that

$$ \lVert{\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u}_{h} - \underline{u}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} \leq (1+\gamma)\lVert{\underline{\mathbf{K}}^{- \frac12} ({\boldsymbol \theta}^{(i+\nu)}_{h} - {\boldsymbol \theta}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})}, $$

which gives an upper bound easily computed and cheaply evaluated in practice even for complex problems (see [21] for details).

Remark 5 (Local indicators for the algebraic error)

In order to estimate the local distribution of the algebraic error (6) using flux reconstructions, one can use the local indicators $\eta ^{(i)}_{\text {alg},K} \hspace {-0.1cm} := \hspace {-0.1cm} \lVert {\underline {\mathbf {K}}^{- \frac 12} ({\boldsymbol \theta }^{(i+\nu )}_{h} \hspace {-0.1cm} - {\boldsymbol \theta }^{(i)}_{h})}\rVert _{\mathrm {L}^{2}(K)} + C_{\mathrm {F}} h_{\Omega }{\lambda }^{-\frac 12}_{\underline {\mathbf {K}}} \lVert {\underline {s}_{h}^{(i+\nu )}}\rVert _{\mathrm {L}^{2}(K)}$ . Relying on the previous discussion, in practice with a sufficiently large ν one can use the local algebraic indicator $\lVert {\underline {\mathbf {K}}^{- \frac 12} ({\boldsymbol \theta }^{(i+\nu )}_{h} - {\boldsymbol \theta }^{(i)}_{h})}\rVert _{\mathrm {L}^{2}(K)}$ which can be the ingredient of our adaptive procedure.

Remark 6 (A posteriori error estimates for the total error)

A computable upper bound can be obtained on the energy norm of the total error $\underline {u} - \underline {u}^{(i)}_{h}$ following the same ideas of Section 6 and Section 6:

$$ \lVert{\underline{\mathbf{K}}^{\frac12}\nabla(\underline{u} - \underline{u}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} \leq \eta^{(i)}_{\text{dis}} + \eta^{(i)}_{\text{alg}} $$

with

$$ \eta^{(i)}_{\text{dis}} := \lVert{\underline{\mathbf{K}}^{\frac12}\nabla \underline{u}^{(i)}_{h} + \underline{\mathbf{K}}^{-\frac12} {\boldsymbol \theta}^{(i)}_{h}}\rVert_{\mathrm{L}^{2}({\Omega})} $$

and

$$ \eta^{(i)}_{\text{alg}} := \lVert{\underline{\mathbf{K}}^{- \frac12} ({\boldsymbol \theta}^{(i+\nu)}_{h} - {\boldsymbol \theta}^{(i)}_{h})}\rVert_{\mathrm{L}^{2}({\Omega})} + C_{\mathrm{F}} h_{\Omega}{\lambda}^{-\frac12}_{\underline{\mathbf{K}}} \lVert{\underline{s}_{h}^{(i+\nu)}}\rVert_{\mathrm{L}^{2}({\Omega})}, $$

see [17] for the full demonstration.

Remark 7 (A posteriori error estimates in the multilevel setting)

There is also another way how to construct upper bounds for the total and algebraic errors without the need of running additional iterations of the algebraic solver. The construction assumes the existence of the hierarchy of meshes, with a global solve on the coarsest mesh (see [16] for more details).