Quasi-optimal mesh sequence construction through Smoothed Adaptive Finite Element Method

TL;DR

This paper introduces S-AFEM, a smoothing-based adaptive finite element method that achieves similar refinement patterns to classical AFEM but with significantly reduced computational cost, validated through analysis and experiments.

Contribution

The paper presents a novel S-AFEM algorithm that replaces algebraic solutions with smoothing iterations, offering efficiency gains while maintaining accuracy.

Findings

Refinement patterns are similar to classical AFEM.

S-AFEM reduces computational cost significantly.

Numerical experiments confirm efficiency and speedup.

Abstract

We propose a new algorithm for Adaptive Finite Element Methods (AFEMs) based on smoothing iterations (S-AFEM), for linear, second-order, elliptic partial differential equations (PDEs). The algorithm is inspired by the ascending phase of the V-cycle multigrid method: we replace accurate algebraic solutions in intermediate cycles of the classical AFEM with the application of a prolongation step, followed by the application of a smoother. Even though these intermediate solutions are far from the exact algebraic solutions, their a-posteriori error estimation produces a refinement pattern that is substantially equivalent to the one that would be generated by classical AFEM, at a considerable fraction of the computational cost. We provide a qualitative analysis of how the error propagates throughout the algorithm, and we present a series of numerical experiments that highlight the efficiency and the computational speedup of S-AFEM.

Figures

Tables

Table 1. Table 1 : Comparison of the computational cost of the solution stage for ten cycles of adaptive refinement using classical AFEM and S-AFEM on bi-linear elements, for the pure diffusion case.

	Peak 2D	L-shaped 2D	Peak 3D	Fichera 3D
First and last solve (same for AFEM and S-AFEM)	0.0187s	0.0601s	32s	101s
AFEM Intermediate solves (CG)	0.0663s	0.219s	76.4s	185s
S-AFEM Intermediate smoothing steps (Richardson)	0.005s	0.00892s	0.252s	0.426s
S-AFEM intermediate speedup	13.26	24.6	303.7	434.3
S-AFEM total speedup	3.59	4.045	3.361	2.819

Equations

$$V_1\subset V_2\subset \cdots \subset V_{\bar{k}}\equiv V_h,$$

$$-\Delta u=f\text{in}\Omega \text{and}u=0\text{on}\partial \Omega ,$$

$$(\nabla u,\nabla v)=(f,v)\forall v\in H_0^1(\Omega ).$$

$$(\nabla u_h,\nabla v_h)=(f,v_h)\forall v_h\in V_h.$$

$$(\nabla (u-u_h),\nabla v_h)=0\forall v_h\in V_h.$$

$$A\mathbf{u}=\mathbf{f}\text{in}{\mathbb{R}}^N,$$

$$(\nabla u_k,\nabla v_k)=(f,v_k)\forall v_k\in V_k,$$

$$A_k{\mathbf{u}}_k={\mathbf{f}}_k$$

$$\Vert \mathcal{R}\{u_h\}{\Vert }_{\star }:=\underset{v\in H_0^1(\Omega )\setminus \{0\}}{\sup }\frac{\mathcal{R}\{u_h\}(v)}{|v{|}_1}=\underset{v\in H_0^1(\Omega )\setminus \{0\}}{\sup }\frac{(f,v)-a(u_h,v)}{|v{|}_1}.$$

$$\Vert e_h\Vert \le C_{\text{rel}}\eta (u_h)+h.o.t{.}_{\text{rel}}$$

$$\eta (u_h)\le C_{\text{eff}}\Vert e_h\Vert +h.o.t{.}_{\text{eff}},$$

$$J_E(v_h)$$

$$J(v_h)$$

$$|u-u_h{|}_1\le C^{\star }({\text{osc}}^2+J^2(u_h){)}^{1/2},$$

$$\underset{\text{total error}}{\underbrace{u-u_h^{\ell }}}=\underset{\text{discretization error}}{\underbrace{(u-u_h)}}+\underset{\text{algebraic error}}{\underbrace{(u_h-u_h^{\ell })}}.$$

$$|u-v_h{|}_1^2\le 2C^2(J^2(v_h)+{\text{osc}}^2)+2C_{\text{intp}}^2|u_h-v_h{|}_1^2,$$

$${\mathbf{u}}_k^{(i+1)}={\mathbf{u}}_k^{(i)}+{\omega }_k({\mathbf{f}}_k-A_k{\mathbf{u}}_k^{(i)})\text{for}i=0,1,\dots ,\ell$$

$$J^2(v_h)\le C_{\star }(|u-v_h{|}_1^2+{\text{osc}}^2)\forall v_h\in V_h,$$

$$J^2(v_h)\le C_1{J}^2(u_h)+C_2|u_h-v_h{|}_1^2+C_3{\text{osc}}^2\forall v_h\in V_h.$$

$$J^2(v_h)$$

$$=C_{\star }(|u-u_h{|}_1^2+|u_h-v_h{|}_1^2+{\text{osc}}^2)$$

$$\le C_{\star }({C^{\star }}^2({\text{osc}}^2+J^2(u_h))+|u_h-v_h{|}_1^2+{\text{osc}}^2)$$

$$=C_{\star }{C^{\star }}^2{J}^2(u_h)+C_{\star }|u_h-v_h{|}_1^2+C_{\star }({C^{\star }}^2+1){\text{osc}}^2$$

$$=C_1{J}^2(u_h)+C_2|u_h-v_h{|}_1^2+C_3{\text{osc}}^2,$$

$${\mathbf{a}}_{k+1}:={\mathbf{u}}_{k+1}-I_k^{k+1}{\mathbf{u}}_k\in {\mathbb{R}}^{N_{k+1}},k=1,\dots ,\bar{k}-1$$

$${\mathbf{e}}_{k+1}^{(\ell )}={M_{k+1}}^{\ell }({\mathbf{a}}_{k+1}+I_k^{k+1}{\mathbf{e}}_k^{(\ell )}),\text{ for }k=1,\dots ,\bar{k}-1.$$

$${\mathbf{u}}_{k+1}^{(0)}$$

$$=I_k^{k+1}{\mathbf{u}}_k-I_k^{k+1}{\mathbf{e}}_k^{(\ell )},$$

$${\mathbf{e}}_{k+1}^{(0)}$$

$$={\mathbf{u}}_{k+1}-I_k^{k+1}{\mathbf{u}}_k+I_k^{k+1}{\mathbf{e}}_k^{(\ell )}$$

$$={\mathbf{a}}_{k+1}+I_k^{k+1}{\mathbf{e}}_k^{(\ell )}.$$

$${\mathbf{e}}_{k+1}^{(\ell )}$$

$$={M_{k+1}}^{\ell }({\mathbf{a}}_{k+1}+I_k^{k+1}{\mathbf{e}}_k^{(\ell )}),$$

$$-\Delta u+\mathbf{\beta }\cdot \nabla u=f$$

$$u=u_g$$

Full text

We provide a qualitative analysis of how the error propagates throughout the algorithm, and we present a series of numerical experiments that highlight the efficiency and the computational speedup of S-AFEM.

{AMS}

65N15, 65N22, 65N30, 65N50, 65N55

1 Introduction

The efficient numerical simulation of complex real-world phenomena requires the use of computationally affordable discrete models. The adaptive finite element method (AFEM) is one such a scheme for the numerical solution of partial differential equations (PDEs) in computational sciences and engineering. In finite element simulations (FEM), the domain of the PDE is discretized into a large set of small and simple domains (the cells or elements) depending on a size parameter $h>0$, and the PDE is transformed into an algebraic system of equations. Rigorous analysis of the numerical method allows one to estimate the discretization error both a priori (giving global bounds on the total discretization error that depend on a global size parameter h), and a posteriori (providing a local distribution of the error on the discretized mesh in terms of known quantities). Classical AFEM consists of successive loops of the steps $\mathcal{S}\!olve\longrightarrow\mathcal{E}\!stimate\longrightarrow\mathcal{M}ark\longrightarrow\mathcal{R}efine$ to decrease the total discretization error, by repeating the FEM solution process ($\mathcal{S}\!olve$) on a mesh that has been refined ($\mathcal{R}efine$) where the a-posteriori analysis ($\mathcal{E}\!stimate$) has shown that the error is larger ($\mathcal{M}ark$).

Intermediate solution steps are instrumental for the construction of the finally adapted grid, and play no role in the final solution, which is the only one that is retained for analysis and processing.

In this work we present a simple yet effective algorithm to reduce the overall computational cost of the AFEM algorithm, by providing a fast procedure for the construction of quasi-optimal mesh sequences that do not require the exact solution of the algebraic problem in the intermediate loops of AFEM.

We refer to this new algorithm as Smoothed Adaptive Finite Element (S-AFEM) method: the $\mathcal{S}\!olve$ step of AFEM in all intermediate loops is replaced by the application of a prolongation step ($\mathcal{P}\!rolongate$), followed by the application of a smoother ($\mathcal{S}\!mooth$):

S-AFEM takes its inspiration from the ascending phase of the V-cycle multigrid method where a sequence of prolongation and smoothing steps is applied to what is considered an algebraically exact solution at the coarsest level. In the multigrid literature, this procedure is used to transfer the low frequency information contained in the coarse solution to a finer –nested– grid, where some steps of a smoothing iteration are applied in order to improve the accuracy of the solution in the high frequency range (see, for example, the classical references [20, 21, 35, 13, 12]). The iteration of this procedure is very effective in providing accurate algebraic solutions in $O(N)$ time, where $N$ is the dimension of the final algebraic system. Even a small number of smoothing iterations is sufficient to eliminate the high frequency error, while the prolongation from coarser grids guarantees the convergence in the low frequency regime, resulting in an overall accurate solution, also when local refinement is present (see, for example, [22]).

The classical AFEM algorithm generates nested grids and subspaces111During AFEM, the grids remain nested if no de-refinement occurs. In the following we will work under this assumption., but the construction of the next (unknown) grid in the sequence still requires an exact algebraic solution on the current grid to trigger the $\mathcal{E}\!stimate$-$\mathcal{M}ark$-$\mathcal{R}efine$ steps. In this paper we show, however, that in many practical situations it is not necessary to use a fully resolved solution in the intermediate steps in order to obtain a good refinement pattern: numerical evidence shows that it is sufficient for the high frequencies of the error to be dumped, in order to identify the next grid in the sequence through the $\mathcal{E}\!stimate$-$\mathcal{M}ark$ steps. In this context, the construction of a grid with excellent approximation properties may require as little as a single ascending step of a V-cycle multigrid method.

It is still not clear how to explain rigorously why the sequence of meshes constructed with S-AFEM is close to the one obtained by classical AFEM. In this work we examine the numerical behaviour of S-AFEM in a family of linear second order elliptic problems, and provide a qualitative analysis, based on classical results from the AFEM and multigrid literature, to point the reader towards currently open questions, and towards possible paths to complete the analysis.

In particular, we consider conforming discretizations of a Poisson problem, and of a class of drift-diffusion problems in 2D and 3D, and we show that

• $\star$

the a-posteriori error estimator applied to the outcome of a single ascending phase of the V-cycle multigrid method triggers a $\mathcal{M}ark$ step where the refinement pattern is substantially equivalent to the one that would be generated by a classical $\mathcal{S}\!olve$ step, at a considerable fraction of the computational cost;

• $\star$

even if the final grid is not exactly identical to the one that would be obtained with the classical AFEM, the accuracy of the final solution is comparable in most cases;

• $\star$

the S-AFEM algorithm is robust with respect to different smoothers, and with respect to different discretization degrees.

The article is organised as follows: we start by describing the general S-AFEM algorithm in Section 2, where the main algorithm is exposed. The connection with the multilevel framework and with classical a posteriori error analysis is made in Section 3, where a qualitative analysis of S-AFEM is presented. Section 4 is dedicated to the numerical validation of the algorithm, and presents a detailed campaign of simulations that shows when the S-AFEM algorithm can be used successfully. Finally, in Section 5 we provide some conclusions and perspectives for future works.

2 The S-AFEM algorithm

We consider a class of linear elliptic, second-order, boundary value problems (BVPs), whose variational formulation reads: seek $u\in V$ s.t. $\mathcal{A}u=f$ in $V$ under suitable boundary conditions, where $(V,\|\bullet\|)$ is a normed Hilbert space defined over a Lipschitz bounded domain $\Omega$, the linear operator $\mathcal{A}:V\rightarrow V^{\star}$ is a second order elliptic operator, and $f\in V^{\star}$ is a given datum. The finite element method provides numerical solutions to the above problem in a finite dimensional solution space $V_{h}\subset V$, typically made up by continuous and piecewise polynomial functions, and transforms the continuous problem above in a discrete model of type $\mathcal{A}_{h}u_{h}=f_{h}$ in $V_{h}$ under suitable boundary conditions, where, e.g., $\mathcal{A}_{h}=\mathcal{A}\mid_{V_{h}}.$ The overall procedure leads to the solution of a (potentially very large) linear algebraic system of equations of type $A\mathbf{u}=\mathbf{f}$ in $\mathbb{R}^{N}$, where $N=\dim(V_{h})$.

Given an initial (coarse) triangulation $\mathcal{T}_{1}$, we consider a (a priori unknown) nested sequence of shape regular triangulations $\mathcal{T}_{k}$, for $k=1,\dots,\bar{k}$, which induces a nested sequence of finite element spaces

[TABLE]

on which we define standard prolongation operators, considering the canonical embedding $i_{k}^{k+1}:V_{k}\hookrightarrow V_{k+1}$ that embeds functions $u_{k}\in V_{k}$ in the space $V_{k+1}$. We denote by $I_{k}^{k+1}:\mathbb{R}^{N_{k}}\rightarrow\mathbb{R}^{N_{k+1}}$ the corresponding discrete matrices and we let $N_{k}:=\dim(V_{k}),$ for $k=1,\dots,\bar{k}$.

The sequence of grids and the solution on the finest grid is computed with the S-AFEM algorithm, defined as:

The motivation behind the strategy at the base of S-AFEM is the numerical observation that classical residual-based a posteriori error estimators [34] used in the $\mathcal{E}\!stimate$ step are mostly insensitive to low frequencies in the solution, as shown in Figure 1 for a benchmark example.

Their application to very inaccurate approximate solutions in the intermediate loops – only capturing high frequency oscillations through a smoother – produces an equally good grid refinement pattern at each step of the S-AFEM algorithm, with an accuracy on the final approximation step that is comparable to the one obtained with the classical AFEM algorithm, at a fraction of the computational cost.

In the S-AFEM algorithm, we capture the smoothest (i.e. less oscillatory) part of the discrete approximation in the first loop ($k=1$), by solving the discrete system exactly on the coarsest level. As the mesh is locally refined from one level to the other, we increase the higher portion of the spectrum of the matrix $A_{k}$. Thanks to the structure of the refinement in typical finite element methods, mostly high frequencies are added to the system, while low frequencies are substantially left unaltered.

Even though the distance between the algebraic approximation $u_{h}^{\ell}$ (coming from S-AFEM algorithm at step $k$) and the exact solution $u$ drifts away during the various steps of the algorithm (curve $|u-u_{h}^{\ell}|_{1}$ in Figure 1), the error estimator evaluated on $u_{h}^{\ell}$ remains substantially attached to the error estimator evaluated on an algebraically exact solution $u_{h}$ computed on the same mesh (comparison between curves Est($u_{h}^{\ell}$) and Est($u_{h}$) in Figure 1).

The first and last loops of the S-AFEM algorithm coincide with those of the classical AFEM. In intermediate S-AFEM loops, however, the solution $u_{h}^{\ell}$ is far from the exact algebraic solution $u_{h}$. These intermediate solutions serve solely to the construction of the final grid, and find no other use in the final computations, therefore their inexactness is irrelevant, provided that the finally adapted grid provides a good approximation. Their role is instrumental in triggering the $\mathcal{E}\!stimate-\mathcal{M}ark-\mathcal{R}efine$ steps.

In our setting, intermediate steps are only computed through a fixed number of smoothing iterations, and have non-negligible algebraic errors. This is in contrast with the common practical assumption made in AFEM, where it is assumed that the $\mathcal{S}\!olve$ step produces the exact solution of the algebraic system. Recent developments dedicated a great deal of effort to account for inexactness of the algebraic approximations and introduce stopping criteria based on the interplay between discretization and algebraic computation in adaptive FEM. Among others, we mention the seminal contributions [10, 23, 4, 5, 31, 26, 30, 29, 17, 28, 25, 16].

Nevertheless, most of this literature focuses on ways to estimate the algebraic error, without really exploiting the other side of the coin: inexact approximate solutions, with large algebraic error, offer large computational savings when used in the correct way. S-AFEM provides a good strategy to exploit this fact.

3 Qualitative analysis of S-AFEM

To fix the ideas, in this section we apply S-AFEM to a model Poisson problem with Richardson iteration as a smoother and discuss the interplay between the algebraic solution in intermediate steps, and the classical a-posteriori error estimator theory.

A larger selection of problem types and smoother algorithms is tested in Section 4 where we compare the application of a fixed number of Richardson iterations or of the conjugate gradient (CG) method for symmetric systems, and of the generalized minimal residual method (GMRES) method for non-symmetric systems coming from the discretization of drift-diffusion problems.

Let $\Omega\subset\mathbb{R}^{d}$ $(d=1,2,3)$ be a bounded, polygonal domain (an open and connected set with polygonal boundary) with Lebesgue and Sobolev spaces $L^{2}(\Omega)$ and $H^{1}_{0}(\Omega)$. We look for the solution $u\in H^{1}_{0}(\Omega)$ such that

[TABLE]

where $f\in L^{2}(\Omega)$ is a given source term. We use standard notations for norms and scalar products in Sobolev spaces (cf. [1]): for $u\in H^{1}_{0}(\Omega)$ we write $|u|_{1}:=(\int_{\Omega}|\nabla u|^{2})^{1/2}$ and denote by $(\cdot,\cdot)$ the $L^{2}(\Omega)$- scalar product with corresponding norm $\|\cdot\|$. The weak form of (2) is to find $u\in H^{1}_{0}(\Omega)$ s.t.

[TABLE]

We consider a shape regular family of triangulations $\{\mathcal{T}_{h}\}_{h}$ of $\Omega$ in the sense of Ciarlet [15], depending on a parameter $h>0$ with shape regularity parameter $C_{\mathcal{T}_{h}}$ consisting of cells $T$ that are convex quadrilaterals in two dimensions, and convex hexahedrons in three dimensions.

The set of all edges/faces $E$ of the cells is denoted by $\mathcal{E}_{h}$ and similarly, $\mathcal{E}_{h,int}:=\mathcal{E}_{h}\setminus\partial\Omega$ is the set of internal edges/faces. We use the Courant finite element space $V_{h}:=\{\varphi\in C^{0}(\Omega)\text{ s.t. }\varphi|_{T}\in P^{k}(T)\quad\forall T\in\mathcal{T}_{h},\varphi=0\text{ on }\partial\Omega\}\subset H^{1}_{0}(\Omega).$ The Galerkin solution $u_{h}\in V_{h}$ is obtained by solving the discrete system

[TABLE]

In exact arithmetic, the discretization error $e_{h}:=u-u_{h}$ satisfies the standard orthogonality condition

[TABLE]

Let $N=\dim(V_{h}),$ the discrete system (4) leads to a linear algebraic system of type

[TABLE]

where $A$ denotes the symmetric positive definite (SPD) stiffness matrix with entries $a_{ij}:=(\nabla\varphi_{j},\nabla\varphi_{i})$ $\forall\,i,j=1,..,N,$ $\mathbf{u}=[u_{1},\dots,u_{N}]^{T}$ denotes the coefficients vector in $\mathbb{R}^{N}$ of the discrete approximation $u_{h}=\sum_{j=1}^{N}u_{j}\varphi_{j}\in V_{h}$ and $\mathbf{f}=[f_{1},\dots,f_{N}]^{T}$ is the vector with entries $f_{j}=(f,\varphi_{j})\,\,\forall j=1,..,N.$

In particular, we assume that an initial (coarse) triangulation $\mathcal{T}_{1}$ is given and we consider a (a priori unknown) nested sequence of shape regular triangulations $\mathcal{T}_{k}$, for $k=1,\dots,\bar{k}$, which induces a nested sequence of finite element spaces $V_{1}\subset V_{2}\subset\dots\subset V_{\bar{k}}$. We let $N_{k}:=\dim(V_{k}),$ for $k=1,\dots,\bar{k}$. By construction, the inequalities $N_{1}<N_{2}<\dots<N_{\bar{k}}$ hold true. The associated discrete systems for each level $k=1,2,\dots,\bar{k}$ read

[TABLE]

and they generate linear systems of type

[TABLE]

of respective dimensions $N_{k}$.

The sequence of meshes in the classical AFEM algorithm is constructed through the $\mathcal{S}\!olve$-$\mathcal{E}\!stimate$-$\mathcal{M}ark$-$\mathcal{R}efine$ steps, where the $\mathcal{S}\!olve$ step should compute an algebraically exact solution of (8). In the $\mathcal{E}\!stimate$ step, standard residual-based a posteriori error estimators are the most widely used. They were first introduced in the context of FEM by Babuška and Rheinboldt in [8] and they have been thereafter widely studied in the literature (for a review, see [34, 2]).

Their derivation is based on the residual functional associated to the Galerkin solution, which is defined as $\mathcal{R}\{u_{h}\}:H^{1}_{0}(\Omega)\longrightarrow\mathbb{R}$, $\mathcal{R}\{u_{h}\}:=(f,\bullet)-a(u_{h},\bullet)$ with corresponding dual norm

[TABLE]

The identity $|e_{h}|_{1}=\|\mathcal{R}\{u_{h}\}\|_{\star}$ leads to reliable and efficient residual-based a posteriori bounds for the discretization error via estimation of the residual function, i.e.,

[TABLE]

and

[TABLE]

where the multiplicative constants $C_{\text{rel}}$ and $C_{\text{eff}}$ are independent on the mesh size and h.o.t. denote oscillations of the right-hand side $f$, which are generally negligible w.r.t. $\|e_{h}\|$.

We use standard residual-based a posteriori error estimators which are locally defined through the jump of the gradient of the discrete approximation across the edges/faces $E$ of the cells, i.e., for a given function $v_{h}\in V_{h},$ define for $E\in\mathcal{E}_{h}$ and $T\in\mathcal{T}_{h}$

[TABLE]

where $\left[\bullet\right]$ indicates the jump of a piecewise continuous function across the edge/face $E$ in normal direction $n_{E}$. A classical upper bound on the discretization error using $J$ is given by the estimate

[TABLE]

where the constant $C^{\star}>0$ depends on the shape of the triangulation, on $\Omega$ and on $\Gamma$, but it is independent of $f$ and of the mesh-sizes $h_{T}$, and osc is an oscillatory term (see [14] for the exact definition of osc and for a proof of (13)).

In S-AFEM, we do not solve the linear systems (8) for $k=2,\ldots,\bar{k}-1$, but only apply a smoother (typically, $\ell$ steps of a smoothing iteration) by taking as an initial guess the prolongation of the approximation from the previous level, obtaining an algebraically inexact approximation $u_{h}^{\ell}$ of $u_{h}$ in all intermediate steps. In this case, the total error in intermediate steps can be written as the sum of two contributions

[TABLE]

A vast literature is dedicated to the extension of the standard residual-based a posteriori error estimator theory to incorporate in some way the algebraic error deriving from an inexact $\mathcal{S}\!olve$ step. We refer to the seminal and investigative paper by Pape$\check{z}$ and Strako$\check{s}$ [30] and the references therein for various approaches.

In particular, in [30] the authors give a detailed proof of a (worst case scenario) residual-based upper bound on the energy norm of the total error

[TABLE]

for arbitrary $v_{h}\in V_{h}$. In order to provide sharper bounds, one should exploit the fact that $v_{h}$ is actually not arbitrary, but it originates from the $\mathcal{S}\!mooth$ step of S-AFEM, i.e., it is the result of a smoothing iteration.

For simplicity of exposition, in this section we use a fixed number of Richardson iterations as a smoother, but other choices are possible (see, for example, the reviews in [12, 19, 35]).

Given $\omega_{k}\in\mathbb{R}$ a fixed parameter and $\mathbf{u}_{k}^{(0)}\in\mathbb{R}^{N_{k}}$ an initial guess, a Richardson smoothing iteration for (8) takes the form

[TABLE]

After $i+1$ iterations, the error $\mathbf{e}_{k}^{(i+1)}:=\mathbf{u}_{k}-\mathbf{u}_{k}^{(i+1)}$ satisfies the error propagation formula $\mathbf{e}_{k}^{(i+1)}=M_{k}\mathbf{e}_{k}^{(i)}=\dots=M_{k}^{i+1}\mathbf{e}_{k}^{(0)}$, where $M_{k}:=Id_{N_{k}}-\omega_{k}A_{k}$ is the Richardson iteration matrix. When using Richardson method as an iterative solver, convergence takes place for $0<\omega_{k}<2/\rho(A_{k}),$ where $\rho(A_{k})$ is the spectral radius of $A_{k}$. The optimal choice of the parameter $\omega_{k}$ is in this case $\omega_{k}=1/\rho(A_{k})$ (see, e.g., [20, 32]).

The high frequency components of the error are reduced by a factor which is close to zero, while the low frequency components of the error are left substantially unaltered. The high frequencies are also the ones that have a greater influence on classical a-posteriori error estimators, justifying the S-AFEM algorithm.

To show the effect of the smoothing iterations on the a-posteriori error estimator, we consider as an example case the Peak Problem in two dimensions as described in Subsection 4.1, and we apply ten cycles of the classical AFEM algorithm using non-preconditioned Richardson iterations for the algebraic resolution of the system with initial guess given by the prolongation of the previous approximation for each cycle.

In Figure 2

we plot the $\ell^{2}$-norm of the residual $\mathbf{r}^{(\ell)}_{k}:=A_{k}\mathbf{e}_{k}^{(\ell)}$ and the value of the a-posteriori error estimator $J(u^{\ell}_{k})$ for all cycles as the Richardson iteration count $\ell$ increases from $1$ to $30$.

The same behaviour is present in every refinement cycle: the first few Richardson iterations induce a rapid drop in the residual norm (due to convergence of the highly oscillatory terms in the solution), while the second part of the iterations converge very slowly, corresponding to the convergence speed of the low frequencies in the solution (typical of Richardson iteration). The estimator, on the other hand, stagnates after very few Richardson iterations (around two or three), suggesting that $J(u_{h}^{\ell})$ is almost the same as $J(u_{h})$ for $\ell\geq 3$, i.e., the error estimator is mainly affected by the highly oscillatory components of the discrete algebraic solution $u_{h}^{\ell}$.

In this respect, classical results of a-posteriori error analysis do not provide sharp bounds for the estimator evaluated on the smoothed algebraic approximation $J(u_{h}^{\ell})$, and this closeness remains an open problem.

One can combine the upper bound (13) with a global lower bound [11] of the estimator evaluated on a generic $v_{h}$, i.e.,

[TABLE]

to relate $J(u_{h}^{\ell})$ with $J(u_{h})$. However, the result (proved in the following Lemma) remains a worst case estimate (similar to (15)), and fails to capture the behaviour that we observe, for example in Figure 1, where $J(u_{h})$ and $J(u_{h}^{\ell})$ are substantially equivalent.

Lemma 3.1**.**

There exist positive constants $C_{1},C_{2},C_{3}$ that only depend on the minimum angle of the triangulation, on $\Omega$, and on $\Gamma$, that are independent of $f,u,u_{h}$, and of the mesh-sizes $h_{T}$ such that

[TABLE]

Proof 3.2**.**

For a given function $v_{h}\in V_{h},$ we decompose $u-v_{h}=(u-u_{h})+(u_{h}-v_{h})$ and we apply the equality $|u-v_{h}|^{2}_{1}=|u-u_{h}|^{2}_{1}+|u_{h}-v_{h}|^{2}_{1}$ (see, e.g. [24]) to the lower bound (17)

[TABLE]

*where we have used the upper bound (13) in (19). *

If we apply Lemma 3.1 with $v_{h}=u^{\ell}_{h}$, we obtain an upper bound on $J^{2}(u^{\ell}_{h})$ in terms of $J^{2}(u_{h})$ and of the algebraic error. A similar result involving the full estimator $\eta(u^{\ell}_{h})$ can be found in [4].

It is still unclear how to improve the upper bound (18) to explain why $J(u_{h}^{\ell})$ and $J(u_{h})$ are as close as the numerical evidence suggests. What remains to be proved is that a sharper estimate on the constant $C_{2}$ of (18) may be obtained, that depends on the frequency content of $v_{h}$, showing that $C_{2}$ is small when $v_{h}=u_{h}^{\ell}$ stems from a smoothing iteration, bringing $J^{2}(u_{h})$ close to $J^{2}(u_{h}^{\ell})$ even though $|u_{h}^{\ell}-u_{h}|^{2}_{1}$ is in fact not small at all.

What can be done, however, is an estimate of the evolution of $|u_{h}^{\ell}-u_{h}|_{1}$ in the intermediate steps of S-AFEM, exploiting classical results of multigrid analysis. Next theorem provides such result when the smoothing iteration is performed using Richardson method.

Theorem 3.3** (Algebraic Error propagation in S-AFEM).**

Let $\mathbf{e}^{(\ell)}_{k}:=\mathbf{u}_{k}-\mathbf{u}^{\ell}_{k}$ denote the algebraic error after $\ell$ smoothing iterations at step $k$ of S-AFEM for $k=2,\dots,\bar{k}-1$. Let $\mathbf{a}_{1}=0$ and define

[TABLE]

*denote the difference between the exact algebraic solution $\mathbf{u}_{k+1}$ at level $k+1$ and the prolongation to level $k+1$ of the exact algebraic solution $\mathbf{u}_{k}$ at level $k$. Then, the following error propagation formula holds true

[TABLE]

Proof 3.4**.**

In the $\mathcal{P}\!rolongate$ step of S-AFEM, the outcome $\mathbf{u}_{k}^{(\ell)}$ of the $\mathcal{S}\!mooth$ procedure at step $k$ is prolongated to step $k+1$, and used as an initial guess in Richardson iteration at step $k+1$. We can write

[TABLE]

and therefore express the initial error at step $k+1$ as

[TABLE]

By applying directly the property of the error propagation formula for Richardson iterations, we obtain the final algebraic error at level $k+1$:

[TABLE]

*which proves the recursive formula. *

Theorem 3.3 shows that the nature of the algebraic error in S-AFEM is the result of $\ell$ smoothing iterations applied to a vector that accumulates the (smoothed) prolongation of the exact algebraic solution coming from step zero.

The rationale behind S-AFEM is then that in the first step we perform a full $\mathcal{S}\!olve$, resulting in a negligible algebraic error, and the only components of the error that we are introducing when prolongating from step $k$ to step $k+1$ are high frequency errors (introduced by local refinement). These, however, are reduced very quickly by $\ell$ steps of smoothing iterations.

The residual algebraic error that persists as S-AFEM proceeds (clearly visible in Figure 1) seems not to have a detrimental effect on $J(u_{h}^{\ell})$. Such an algebraic error is probably confined on medium to low frequencies, and shows no noticeable effect on the $\mathcal{E}\!stimate$ and $\mathcal{M}ark$ steps of S-AFEM.

Although the value we plot in Figure 2 for the estimator is a global one, and gives no information on the distribution of the local estimator on the grid, it is a good hint that the overall behaviour of such distribution will not be changing too much after the first few Richardson iterations. We show some numerical evidence that this is actually the case in the numerical validation provided in Section 4.

To summarize the idea behind S-AFEM, we argue that in the intermediate AFEM cycles it is not necessary to solve exactly the discrete system. What matters instead is to capture accurately the highly oscillatory components of the discrete approximation. Low frequency components may have an influence on the error estimator, however, this is mostly a global influence, that has a small effect on the cells that will actually be marked for refinement in the $\mathcal{M}ark$ step. As an example, consider Figure 13, where a Peak problem in 3D is solved using both AFEM and S-AFEM. The estimator evaluated on $u_{h}^{\ell}$ in this case is farther away from the one evaluated on $u_{h}$ w.r.t. the same problem in two dimensions, but the convergence rate of the solution obtained with the sequence of grids constructed with S-AFEM is still the optimal one, and comparable to the one obtained with AFEM at a fraction of the computational cost.

4 Numerical validation

We consider a class of drift-diffusion problems of the following form

[TABLE]

in two and three dimensions.

We compare the classical AFEM algorithm with S-AFEM based on three different smoothing strategies: Richardson iterations (Richardson), conjugate gradient iterations (CG), and generalized minimal residual iterations (GMRES) in the symmetric case ($\boldsymbol{\beta}=0$), and GMRES alone in the non-symmetric case ($\boldsymbol{\beta}\neq 0$), and for different discretization degrees.

We test two classical experiments used to benchmark adaptive finite element methods when $\boldsymbol{\beta}$ is zero, and a simple drift-diffusion problem with constant transport term $\boldsymbol{\beta}$ to test S-AFEM in general drift-diffusion problems.

The numerical results presented in this paper were realized using a custom C++ code based on the deal.II library [3, 6, 7], and on the deal2lkit library [33]. In the $\mathcal{M}ark$ step, we use the classical Dörfler marking strategy [18]: for any level $k$ we mark for refinement the subset of elements

[TABLE]

where $L$ is a threshold error, defined as the largest value such that

[TABLE]

The parameter $\Theta$ is such that $0\leq\Theta\leq 1$, where $\Theta=1$ corresponds to an almost uniform refinement, while $\Theta=0$ corresponds to no refinement. In our numerical tests, unless otherwise stated, we set $\Theta=0.3$. The refinement strategy that we adopt in this work is based on the use of hanging nodes [9].

The results presented in this section are a subset of the full campaign of simulations presented in [27].

4.1 Two-dimensional examples: pure diffusion, bi-linear case, Richardson smoother

Smooth domain, peak right hand side

The first example we consider consists in solving the model problem with no transport term ($\boldsymbol{\beta}=0$) on a square domain, with a custom forcing term that contains a peak in a specified point in the domain, forcing the exact solution to be

[TABLE]

as shown in Figure 3.

L-shaped domain, smooth right hand side

In the second two-dimensional test case in pure diffusion, we consider a L-shaped domain, i.e., a square where the upper right corner is removed, and the reentrant corner coincides with the origin. No forcing term is added to the problem, but the boundary conditions are set so that the following exact solution is obtained (when expressed in polar coordinates)

[TABLE]

as shown in Figure 4.

In both cases, we apply ten cycles of classical AFEM and of S-AFEM, respectively. For the AFEM algorithm, we use the CG method as iterative solver, with an algebraic multigrid preconditioner (AMG), and we iterate until the $\ell^{2}$-norm of the residual is below a tolerance of $10^{-12}$ for each cycle. For S-AFEM, we modify the intermediate cycles and we only apply three Richardson iterations. For reference, we report a comparison between the cells marked for refinement by AFEM and S-AFEM after four cycles for the Peak Problem in Figure 5 and after nine cycles for the L-shaped domain Problem in Figure 6.

In both cases, the set of marked cells, although different in some areas, produces a refined grid that is very similar between the classical AFEM and the S-AFEM, and where the accuracy of the final solution is comparable.

In Figures 7 and 9

we compare the values of the global estimators $J(u_{h})$ and $J(u_{h}^{\ell})$ and of the $H_{1}$ semi-norm of the total errors for each cycle for the Peak Problem, and for the L-shaped domain Problem respectively, when using S-AFEM. For reference, Figures 8 and 10 show the error and the estimator in the classical AFEM algorithm for the two examples. Notice that the first step of AFEM and of S-AFEM are the same. The last step in the S-AFEM case shows comparable results as in the AFEM algorithm for both examples.

Notice that in S-AFEM the value of the global estimator is almost the same of the one that would be obtained by solving using CG preconditined with AMG ($J(u_{h})$ in Figures 7 and 9), showing that in the two dimensional setting the error estimator (12) is mainly affected by the high frequencies of the discrete solution, which are well captured with just a few Richardson iterations. On the other hand, the total error increases in the intermediate cycles, due to the algebraic error that has been accumulated by applying smoothing iterations instead of solving the algebraic problem until convergence, as quantified by Theorem 3.3. This error measures the distance between the exact algebraic solution and the smooth non-oscillatory components of the approximate solution that are not captured by Richardson iteration, and have little or no influence on the error estimator, and therefore on the generated grid. After ten cycles, we solve the algebraic problem until converge using CG and AMG, as in the first cycle, and we obtain a solution whose error is controlled by the estimator, as expected.

4.2 Three-dimensional examples: pure diffusion, bi-linear case, Richardson smoother

Smooth domain, peak right hand side

The first three-dimensional test case that we propose is a model problem on a cube domain, where the forcing term contains a peak in a specified point that forces the exact solution to be given by

[TABLE]

as shown in Figure 11.

Fichera corner domain, smooth right hand side

In the second three-dimensional example, we consider the classic Fichera corner domain, i.e., a cube where the upper right corner is removed, and the reentrant corner coincides with the origin. We set the exact solution to be

[TABLE]

and we add a forcing term that induces the above exact solution

as shown in Figure 12.

In both examples, the estimator applied to the algebraic solution after three smoothing steps (see Figures 13

and 15)

seems to be more sensitive to the low frequency content of $u_{h}^{\ell}$.

For reference, Figures 14 and 16 show the error and the estimator in the classical AFEM algorithm for the two examples. In the three-dimensional case Lemma 3.1 and Theorem 3.3 provide a sharper estimate, and we don’t observe the same behaviour as in the two-dimensional case (i.e., $J(u_{h}^{\ell})$ does not seem to remain close to $J(u_{h})$). Nonetheless, the difference in accuracy at the final step between AFEM and S-AFEM is negligible also in this case, showing that the differences in the refinement patterns between AFEM and S-AFEM remain small and do not hinder the final accuracy.

4.3 Robustness with respect to approximation degree, smoothing algorithms, and number of smoothing steps

We now consider different variants of our S-AFEM algorithm, where we apply a different number of smoothing iterations, and different smoother types in the intermediate steps (respectively, Richardson iteration, CG, and GMRES smoothers), for high order finite element discretizations of the pure diffusion case.

We apply both AFEM and S-AFEM to the two-dimensional Corner Problem (29) and we show a comparison for different fixed FEM degrees, as $deg=1,2,3$, and for different choices of smoothers for the intermediate cycles, respectively Richardson iteration, the CG method, and the GMRES method. For all cases, we plot the value of the error estimator $J$ and the value of the $|\cdot|_{1}$ semi-norm of the total error, as the number of smoothing iterations $\ell$ increases from $1$ to $5$ in Figures 17-25.

For bi-linear finite elements (cf. Figures 17–19), all considered smoothers (Richardson, CG, and GMRES) turned out to be good. In all cases, the estimator $J(u_{h}^{\ell})$ with $\ell=1,3,5$ exhibits the same behaviour (same order of convergence) of the estimator $J(u_{h})$, showing that one or two smoothing iterations would be enough for the intermediate cycles. When we look at the total error, the CG behaves better leading to less error accumulated at the intermediate levels as shown in Figure 18, while Richardson behaves the worst (cf. Figure 17). Nevertheless, in all cases the accuracy of the final approximation for the last cycle obtained by S-AFEM, is almost the same to the one that is generated by classical AFEM.

As the polynomial degree increases, we observed that the Dörfler marking strategy does not provide good refinement patterns for the different problems, unless a fine tuning is made on the marking parameter. Using the same value for $\Theta$ used for degree one, no cells are marked for refinement in higher order finite elements, making the choice for this parameter too much problem dependent and polynomial degree dependent. As an alternative marking strategy, we opted for a marking criterion where a fraction of $1/3$ of the cells with the largest error indicators are selected for refinement, leading to an increase of the number of degrees of freedom of roughly a factor of two in each refinement cycle, independently on the problem type.

For $deg=2$, all smoothers work well exhibiting a quasi optimal convergence order compared to classical AFEM, and the accuracy of the final approximation is almost the same, as shown in Figures 20, 21, and 22.

For $deg=3$ and higher (only $deg=3$ is shown here; see [27] for the full set of simulations), Richardson iteration turns out to be a bad smoother for S-AFEM, unless further tuning of the relaxation parameter $\omega$ is performed. Although $J(u_{h}^{\ell}$) seems to exhibit the same behaviour of $J(u_{h})$, as the smoothing iteration count $\ell$ increases, contrarily to what one might expect, the value of the estimator increases with increasing DOFs (Figure 23), showing that our selection of $\omega$ may not be correct for these problems, and we should estimate in a better way the spectral radius of the final matrix $A$, and modify $\omega$ accordingly.

On the other hand, both the CG method and the GMRES method turn out to be good smoothers for S-AFEM, without the need to tune any parameter, as evidenced in Figures 21, 22, 24, 25. In all cases, in fact, $J(u_{h}^{\ell})$, for $\ell\geq 2$, shows the same optimal convergence rate as the $J(u_{h})$ obtained by the classical AFEM, and although the total error at the intermediate cycles is evident, the accuracy of the final approximations is almost the same.

4.4 Two dimensional drift-diffusion problem

We proceed with testing the accuracy of S-AFEM for a class of drift-diffusion problems where the transport term $\boldsymbol{\beta}$ is non-negligible. We consider a two-dimensional problem, where $\boldsymbol{\beta}=(\beta,\beta)^{T}$, and the scalar parameter $\beta$ takes the values $1,10$ and $50$. In particular, we impose boundary conditions and forcing terms so that the exact solution is

[TABLE]

Since the linear system associated to Problem (25) when $\boldsymbol{\beta}\neq 0$ is not symmetric, we use as a solver and as a smoother the GMRES method (cf. eg. [24]). We apply classic AFEM and S-AFEM for both bi-linear and higher order finite element discretizations for $deg=1,2,3$. We plot the value of the estimator $J$ and the value of the $|\cdot|_{1}$ semi-norm of the total error, using a fixed number of GMRES iterations as a smoother, with $\ell=1,3,5$.

For all the cases where the transport term $\boldsymbol{\beta}=(1,1)^{T}$ is small, the behaviour of the estimator for S-AFEM is exactly the same as the one given by AFEM, for the case $deg=1$, as shown in Figure 26, while it approaches it as the GMRES smoothing iteration count increases for higher order FEM discretizations (i.e. for $deg\geq 2$), as shown in Figures 27, and 28. However, in all cases, the accuracy of the final approximations is almost the same.

For the choices of the transport term $\boldsymbol{\beta}=(10,10)^{T}$ (corresponding to moderate transport) and $\boldsymbol{\beta}=(50,50)^{T}$ (corresponding to large transport) again the behaviour of the estimator for S-AFEM is exactly the same as the one given by AFEM for the case $deg=1$, as shown in Figures 29 and 32, while it approaches it as the GMRES smoothing iteration count increases for higher order FEM discretizations, as evidenced in Figures 30, 31, 33, and 34. In all cases, however, the accuracy of the final approximations is almost the same, showing that S-AFEM turns out to be a good method also for drift-diffusion problems.

4.5 Computational costs

In the following table we show a comparison of the computational cost associated to the classical AFEM and to the smoothed AFEM in the pure diffusion case, for the first four examples we presented in this section.

The results were obtained on a 2.8 GHz Intel Core i7 with 4 cores and 16GB of RAM, using MPI parallelization on all 4 cores.

Table 1 only shows the comparison between AFEM and S-AFEM in the solve phase, where S-AFEM is always faster than AFEM, offering an average speedup of a factor three. In the table we compare the computational cost of all intermediate cycles in S-AFEM (Intermediate smoothing steps (Richardson) in the table), with the corresponding computational cost for standard AFEM (Intermediate solves (CG) in the table). The first and last solve are the same in the two algorithms, and are reported to provide a scaling with respect to the total computational cost of the solution phase in the program. Other phases (like graphical output, mesh setup, assembling setup, and error estimation) are not shown since they are identical in the two algorithms.

5 Conclusions

This work proposes a new smoothed algorithm for adaptive finite element methods (S-AFEM), inspired by multilevel techniques, where the exact algebraic solution in intermediate steps is replaced by the application of a prolongation step, followed by a fixed number of smoothing steps.

The main argument behind the S-AFEM algorithm is that the combined application of the $\mathcal{E}\!stimate$-$\mathcal{M}ark$ steps of AFEM is largely insensitive to substantial algebraic errors in low frequencies. Indeed, even though the intermediate solutions produced by S-AFEM are far from the exact algebraic solutions, we show that their a posteriori error estimation produces a refinement pattern for each cycle that is substantially equivalent to the one that would be generated by classical AFEM, leading roughly to the same set of cells marked for refinement.

Our strategy is based on solving exactly the problem at the coarsest level and at the finest level, and then applying the $\mathcal{E}\!stimate$-$\mathcal{M}ark$-$\mathcal{R}efine$ steps directly to the result of a smoother ($\mathcal{S}\!mooth$) in intermediate levels.

We provide numerical evidences that the S-AFEM strategy is competitive in cost and accuracy by considering some variants of our algorithm, where different smoothers are used in the intermediate cycles (respectively Richardson iteration, the CG method, and the GMRES method).

We conclude that, in general, CG and GMRES act as robust smoothers in S-AFEM also for high order approximations, and for non-symmetric problems, like, for example, drift-diffusion problems with dominant transport.

Our numerical evidences show that two or three smoothing iterations are enough for the two-dimensional case, while three-dimensional problems require from five up to seven smoothing iteration in order to produce good final approximations, independently on the polynomial degree of the finite element approximation.

We analyzed the error propagation properties of the S-AFEM algorithm, and provided a bound on the a-posteriori error estimator applied to the approximated algebraic solution. The results are not sharp, and do not provide a definitive answer on the convergence of the final S-AFEM solution to the AFEM one, but could be used as a ground state for further investigation, which is currently ongoing.

Acknowledgments

LH is partially supported by the Italian Ministry of Instruction, University and Research (MIUR), under the 2017 PRIN project NA-FROM-PDEs MIUR PE1, “Numerical Analysis for Full and Reduced Order Methods for the efficient and accurate solution of complex systems governed by Partial Differential Equations”. OM is thankful to Durham University for the hospitality during her visiting research there.