Optimal Preconditioners for Finite Element Approximations of Convection-Diffusion Equations on structured meshes ††thanks: The work of the first author was partially supported by MIUR, grant number and the work of the second author and third author was partially supported by MIUR, grant number 2006017542 and 20083KLJEZ.
The paper is devoted to the spectral analysis of effective preconditioners for linear systems obtained via a Finite Element approximation to diffusion-dominated convection-diffusion equations. We consider a model setting in which the structured finite element partition is made by equi-lateral triangles. Under such assumptions, if the problem is coercive, and the diffusive and convective coefficients are regular enough, then the proposed preconditioned matrix sequences exhibit a strong clustering at unity, the preconditioning matrix sequence and the original matrix sequence are spectrally equivalent, and the eigenvector matrices have a mild conditioning. The obtained results allow to show the optimality of the related preconditioned Krylov methods. The interest of such a study relies on the observation that automatic grid generators tend to construct equi-lateral triangles when the mesh is fine enough. Numerical tests, both on the model setting and in the non-structured case, show the effectiveness of the proposal and the correctness of the theoretical findings.
Key words. Matrix sequences, clustering, preconditioning, non-Hermitian matrix, Krylov methods, Finite Element approximations
AMS subject classifications. 65F10, 65N22, 15A18, 15A12, 47B65
The paper is concerned with the spectral and computational analysis of effective preconditioners for linear systems arising from Finite Element approximations to the elliptic convection-diffusion problem
with domain of . We consider a model setting in which the structured finite element partition is made by equi-lateral triangles. The interest of such a partition relies on the observation that automatic grid generators tend to construct equi-lateral triangles when the mesh is fine enough.
The analysis is performed having in mind two popular preconditioned Krylov methods. More precisely, we analyze the performances of the Preconditioned Conjugate Gradient (PCG) method in the case of the diffusion problem and of the Preconditioned Generalized Minimal Residual (PGMRES) in the case of the convection-diffusion problem.
We define the preconditioner as a combination of a basic (projected) Toeplitz matrix times diagonal structures. The diagonal part takes into account the variable coefficients in the operator of (LABEL:eq:modello), and especially the diffusion coefficient , while the (projected) Toeplitz part derives from a special approximation of (LABEL:eq:modello) when setting the diffusion coefficient to and the convective velocity field to . Under such assumptions, if the problem is coercive, and the diffusive and convective coefficients are regular enough, then the proposed preconditioned matrix sequences have a strong clustering at unity, the preconditioning matrix sequence and the original matrix sequence are spectrally equivalent, and the eigenvector matrices have a mild conditioning. The obtained results allow to show the optimality of the related preconditioned Krylov methods. It is important to stress that interest of such a study relies on the observation that automatic grid generators tend to construct equi-lateral triangles when the mesh is fine enough. Numerical tests, both on the model setting and in the non-structured case, show the effectiveness of the proposal and the correctness of the theoretical findings.
The outline of the paper is as follows. In Section LABEL:sez:fem we report a brief description of the FE approximation of convection-diffusion equations and the preconditioner definition. Section LABEL:sez:clustering is devoted to the spectral analysis of the underlying preconditioned matrix sequences, in the case of structured uniform meshes. In Section LABEL:sez:numerical_tests, after a preliminary discussion on complexity issues, selected numerical tests illustrate the convergence properties stated in the former section and their extension under weakened assumption or in the case of unstructured meshes. A final Section LABEL:sez:conclusions deals with perspectives and future works.
2 Finite Element approximation and Preconditioning Strategy
Problem (LABEL:eq:modello) can be stated in variational form as follows:
where is the space of square integrable functions, with square integrable weak derivatives vanishing on . We assume that is a polygonal domain and we make the following hypotheses on the coefficients:
The previous assumptions guarantee existence and uniqueness for problem (LABEL:eq:formulazione_variazionale) and hence the existence and uniqueness of the (weak) solution for problem (LABEL:eq:modello).
For the sake of simplicity, we restrict ourselves to linear finite element approximation of problem (LABEL:eq:formulazione_variazionale). To this end, let be a usual finite element partition of into triangles, with and . Let be the space of linear finite elements, i.e.
The finite element approximation of problem (LABEL:eq:formulazione_variazionale) reads:
For each internal node of the mesh , let be such that , and if . Then, the collection of all ’s is a base for . We will denote by the number of the internal nodes of , which corresponds to the dimension of . Then, we write as and the variational equation (LABEL:eq:formulazione_variazionale_fe) becomes an algebraic linear system:
According to these notations and definitions, the algebraic equations in (LABEL:eq:modello_discreto) can be rewritten in matrix form as the linear system
where and represent the approximation of the diffusive term and approximation of the convective term, respectively. More precisely, we have
where suitable quadrature formula are considered in the case of non constant coefficient functions and .
As well known, the main drawback in the linear system resolution is due to the asymptotical ill-conditioning (i.e. very large for large dimensions), so that preconditioning is highly recommended. Hereafter, we refer to a preconditioning strategy previously analyzed in the case of FD/FE approximations of the diffusion problem [12, 15, 16, 18, 19, 17] and recently applied to FD/FE approximations [3, 4, 10] of (LABEL:eq:modello) with respect to the Preconditioned Hermitian and Skew-Hermitian Splitting (PHSS) method [2, 3]. More precisely, the preconditioning matrix sequence is defined as
where , i.e., the suitable scaled main diagonal of and clearly equals .
The computational aspects of this preconditioning strategy with respect to Krylov methods will be discussed later in section LABEL:sez:complexity_issues. Here, preliminarily we want to stress as the preconditioner is tuned only with respect to the diffusion matrix : in other words, we are implicity assuming that the convection phenomenon is not dominant, and no stabilization is required in order to avoid spurious oscillations into the solution.
Moreover, the spectral analysis is performed in the non-Hermitian case by referring to the Hermitian and skew-Hermitian (HSS) decomposition of (that can be performed on any single elementary matrix related to by considering the standard assembling procedure).
According to the definition, the HSS decomposition is given by
since by definition, the diffusion term is a
Hermitian matrix and does not contribute to the skew-Hermitian
part of . Notice also that if
More in general, the matrix is symmetric and positive definite
whenever . Indeed, without the condition
, the matrix does not have a definite sign in general:
in fact, a specific analysis of the involved constants is required in order to guarantee the nonnegativity of
the term . Moreover,
the Lemma below allows to obtain further information regarding such a spectral assumption,
where indicate both the usual vector norms and the induced matrix
 Let be the matrix sequence defined
according to (LABEL:eq:def_E_nb).
Under the assumptions in (LABEL:eq:ipotesi_coefficienti), then we find
with absolute positive constant only depending on and .
The claim holds both in the case in which the matrix elements in
(LABEL:eq:def_psi) are evaluated exactly and whenever a quadrature formula
with error is considered for approximating the involved integrals.
Hereafter, we will denote by , , the matrix sequence associated to a family of meshes , with decreasing finesse parameter . As customary, the whole preconditioning analysis will refer to a matrix sequence instead to a single matrix, since the goal is to quantify the difficulty of the linear system resolution in relation to the accuracy of the chosen approximation scheme.
3 Spectral analysis and clustering properties in the case of structured uniform meshes
The aim of this section is to analyze the spectral properties of the preconditioned matrix sequences in the case of some special domains partitioned with structured uniform meshes, so that spectral tools derived from Toeplitz theory [6, 11] can be successfully applied. The applicative interest of the considered type of domains will be motivated in the short.
Indeed, let be a variate Lebesgue integrable function defined over . By referring to the Fourier coefficients of this function , called generating function,
with , one can build the sequence of Toeplitz matrices . The matrix is said to be the Toeplitz matrix of order , , generated by .
Now, the spectral properties of the matrix sequence are completely understood and characterized in terms of the underlying generating functions (see, e.g.,  for more details). In the following, with respect to the latter claim, we will refer to the notion of equivalent generating functions. In such a respect, we claim that two nonnegative function and defined over a domain are equivalent if and only if and , where means that there exists a pure positive constant such that almost everywhere on .
The main motivation of this paper lies in the analysis of the template case reported in Figure LABEL:fig:esagono_esagoni.a, where we consider a partition into equilateral triangles. It is self-evident that such a problem represents just an academic example. However, it is a fact that a professional mesh generator will locally produce a partitioning, which is “asymptotically” similar to the one reported in the figure.
Our goal is to prove the PCG optimality with respect to the diffusion problem approximation, i.e., the number of iterations for reaching the solution, within a fixed accuracy, can be bounded from above by a constant independent of the dimension . Some additional results about PGMRES convergence properties in the case of the convection-diffusion problem approximation are also reported, both in the case of the hexagonal domain with a structured uniform mesh as in Figure LABEL:fig:esagono_esagoni.a and in the case with a structured uniform mesh as in Figure LABEL:fig:esagono_esagoni.b, so extending previous results proved in  with respect to the PHSS method .
The spectral analysis makes reference to the following definition.
 Let be a sequence of Hermitian matrices of increasing dimensions . The sequence is clustered at in the eigenvalue sense if for any , . The sequence is properly (or strongly) clustered at if for any the number of the eigenvalues of not belonging to can be bounded by a pure constant eventually depending on , but not on . In the case of non-Hermitian matrices, the same definition can be extended by considering the complex disk , instead of the interval .
3.1 Diffusion Equations
We start by considering the simpler diffusion problem, i.e., the analysis concerns the Hermitian matrix sequences , .
Due to the very special choice of the domain , the matrices arising in the case with structured meshes as in Figure LABEL:fig:esagono_esagoni.a result to be a proper projection of Toeplitz matrices generated by the function , , related to a bigger parallelogram shaped domain containing (see Figure LABEL:fig:esagono_esagoni_parallelogrammi), i.e., , with number of the internal nodes. Here, , is a proper projection matrix simply cutting those rows (columns) referring to nodes belonging to , but not to . Thus, the matrix is full rank, i.e., , and .
In the same way, we can also consider the Toeplitz matrix with the same generating function arising when considering a smaller parallelogram shaped domain contained in , i.e., , with , proper projection matrix and such that , and .
These embedding arguments allow to bound, according to the min-max principle (see  for more details on the use of this proof technique), the minimal eigenvalue of the matrices as follows
A first technical step in our spectral analysis concerns relationships between this generating function and the more classical generating function arising in the case of FE approximations on a square with Friedrichs-Keller meshes (see Figure LABEL:fig:esagono_esagoni.b), or standard FD discretizations.
It can be easily observed that these two function are equivalent, in the sense of the previously reported definition, since we have
Thus, because is a matrix-valued linear positive operator (LPO) for every , the Toeplitz matrix sequences generated by this equivalent functions result to be spectrally equivalent, i.e., for any
furthermore, due to the strict positivity of every (see  for the precise definition), since is not identically zero, we have that
is positive definite, and since is not identically zero, we find that
is also positive definite. Here, we are referring to the standard ordering relation between Hermitian matrices, i.e., the notation , with and Hermitian matrices, means that is nonnegative definite.
An interesting remark pertains to the fact that the function is the most natural one from the FE point of view, since no contribution are lost owing to the gradient orthogonality as instead in the case related to ; nevertheless its relationships with the function can be fully exploited in performing the spectral analysis. More precisely, from (LABEL:eq:lmin_proiezione) and (LABEL:eq:rel_toeplitz) and taking into account (see e.g. ) that
Following the very same reasoning, we also find that
Since, by the embedding argument, both and are asymptotic to , it follows that . Finally, following the same analysis for the maximal eigenvalue, we find
where, by , we know that
and hence the spectral condition number of grows as i.e.
where the constant hidden in the previous relation is mild and can be easily estimated.
It is worth stressing that the same matrix can also be considered in a more general setting. In fact, the matrix sequence can also be defined as , since again each internal node in is a vertex of the same constant number of triangles. Therefore, by referring to projection arguments, the spectral analysis can be equivalently performed both on the matrix sequence and .
No matter about this choice, we make use of a second technical step, which is based on standard Taylor’s expansions.
Let , with and hexagonal domain partitioned as in Figure LABEL:fig:esagono_esagoni.a. For any such that there exists a proper such that the Taylor’s expansions centered at a proper have the form
where and are constants independent of .
In the cases at hand, the validity of this claim is just a direct check: the key point relies in the symmetry properties induced by the structured uniform nature of the considered meshes, both in the case of and . Hereafter, we give some detail with respect to the case of the hexagonal domain partitioned as in Figure LABEL:fig:esagono_esagoni.a, see  for the proof in the case . Thus, let’s consider the Taylor’s expansion centered at a proper . By calling and similarly denoting the derivatives of the diffusion coefficient, it holds that
with and denoting the barycenter of the triangle . By referring to the assembling procedure, we will choose as center of the related Taylor’s expansion in the case of Figure LABEL:fig:simmetrie.b-LABEL:fig:simmetrie.d the point marked by . The symmetric position of the involved triangle barycenters (marked by ), with respect to allows to claim (LABEL:eq:exp.a). The same holds true for (LABEL:eq:exp.b) by considering each pertinent point in Figure LABEL:fig:simmetrie.a and, analogously, for (LABEL:eq:exp.c).
 Let’s considering the following notation of scaled matrix of a given matrix . Under the assumptions of Lemma LABEL:lemma:Taylor_expansion_esagoni the following representation holds true
where and are uniformly bounded in spectral norm and have the same pattern as .
We are now ready to prove the PCG optimality.
Let and be the Hermitian positive definite matrix sequences defined according to (LABEL:eq:def_A) and (LABEL:eq:def_P) in the case of the hexagonal domain partitioned as in Figure LABEL:fig:esagono_esagoni.a. Under the assumptions in (LABEL:eq:ipotesi_coefficienti), the sequence is properly clustered at . Moreover, for any all the eigenvalues of belong to an interval well separated from zero [Spectral equivalence property].
Proof. Since , with such that , and , we have
Since is full rank, it is evident that the spectral behavior of
is in principle better than the one of , to which we can address the spectral analysis. Thus, the proof technique refers to [17, 19] and the very key step is given by the relations outlined in (LABEL:eq:lmin_toeplitz_bottom) and (LABEL:eq:lmin_toeplitz_top).
It is worth stressing that the previous claims can also be proved by considering the sequence , instead of , simply by directly referring to the quoted asymptotical expansions. The interest may concern the analysis of the same uniform mesh on more general domains.
Moreover, the same spectral properties has been proved in the case of uniform structured meshes as in Figure LABEL:fig:esagono_esagoni.b in .
3.2 Convection-Diffusion Equations
The natural extension of the claim in Theorem LABEL:teo:clu+se_A_esagoni, in the case of the matrix sequence with , can be proved under the additional assumptions of Lemma LABEL:lemma:normaE, in perfect agreement with Theorem 5.3 in .
Let and be the Hermitian positive definite matrix sequences defined according to (LABEL:eq:def_ReA) and (LABEL:eq:def_P) in the case of the hexagonal domain partitioned as in Figure LABEL:fig:esagono_esagoni.a. Under the assumptions in (LABEL:eq:ipotesi_coefficienti), the sequence is properly clustered at . Moreover, for any all the eigenvalues of belong to an interval well separated from zero [Spectral equivalence property].
The claim holds both in the case in which the matrix elements in (LABEL:eq:def_theta) and (LABEL:eq:def_psi) are evaluated exactly and whenever a quadrature formula with error is considered to approximate the involved integrals.
Proof. The proof can be done verbatim as in , since the key point is proved by referring to relations in (LABEL:eq:lmin_toeplitz_bottom) and (LABEL:eq:lmin_toeplitz_top).
Let and be the
Hermitian matrix sequences defined according to
(LABEL:eq:def_ImA) and (LABEL:eq:def_P) in the case
of the hexagonal domain partitioned as in Figure
Under the assumptions in (LABEL:eq:ipotesi_coefficienti), the sequence is spectrally bounded and properly clustered at with respect to the eigenvalues. The claim holds both in the case in which the matrix elements in (LABEL:eq:def_psi) are evaluated exactly and whenever a quadrature formula with error is considered to approximate the involved integrals.
Proof. This result has been proved in  with respect to Friedrichs-Keller triangulations, but it can easily be extended to the case of matrix sequences arising in the case of a FE partitioning as in Figure LABEL:fig:esagono_esagoni.a by using the same arguments.
On the basis of these two splitted spectral results, we can easily obtain the spectral description of the whole preconditioned matrix sequence , according to the theorem below. The proof technique refers to an analogous result proved in  with respect to FD discretizations of convection-diffusion equations.
Let and be the matrix sequences defined according to (LABEL:eq:def_A) and (LABEL:eq:def_P), both in the case of the hexagonal domain with a structured uniform mesh as in Figure LABEL:fig:esagono_esagoni.a and in the case with a structured uniform mesh as in Figure LABEL:fig:esagono_esagoni.b
Under the assumptions in (LABEL:eq:ipotesi_coefficienti), the sequence is properly clustered at with respect to the eigenvalues. In addition, these eigenvalues all belong to a uniformly bounded rectangle with positive real part, well separated from zero.
The claim holds both in the case in which the matrix elements in (LABEL:eq:def_theta) and (LABEL:eq:def_psi) are evaluated exactly and whenever a quadrature formula with error is considered to approximate the involved integrals.
Proof. A localization result for the eigenvalues of the sequence can be stated simply by referring to the properties of the field of values of the preconditioned matrix: in fact, for any
Since for any
the claimed localization result is obtained as a consequence of those previously proved in Theorems LABEL:teo:clu+se_ReA_esagoni and LABEL:teo:clu+sb_ImA. Moreover, the same localization results together with the claimed clustering properties allow to apply Theorem 4.3 in  and the proof is concluded.
The subsequent result, giving estimates on the condition number of the eigenvector matrix of the preconditioned structure, is of interest in the study of the convergence speed of the GMRES. In particular a logarithmic growth in number of iteration is associated to a polynomial bound in the spectral condition number: such a bound is precisely given below.
Let and constant functions and let and be the matrix sequences defined as in Theorem LABEL:teo:cluA
can be diagonalized as
where the matrix of the eigenvectors can be chosen such that .
Proof. Under the quoted assumption it holds that and is skew-symmetric. Thus, we have
with skew-symmetric matrix. Therefore, the matrix is a normal normal, so that , where is diagonal and is unitary. Consequently,
with eigenvector matrix. Due to the fact that is unitary, we have and finally by invoking the key relation (LABEL:eq:cond-an1).
4 Numerical tests
The section is divided into two parts. In the first we briefly discuss the complexity features of our preconditioning proposals. In the second part we report and critically discuss few numerical experiments in which the meshes are both structured and unstructured, while the regularity features required by the theoretical analysis are somehow relaxed.
4.1 Complexity issues
First of all, we report some remarks about the computational costs of the proposed iterative procedure.
The main idea is that
such a technique is easily applicable. In fact, a Krylov method is
considered, so simply requiring a matrix vector routine for sparse
matrices and a solver for the chosen preconditioner.
Therefore, since the preconditioner is defined as , where , the solution of the linear system in (LABEL:eq:modello_discreto) with matrix is reduced to computations involving diagonals and the matrix ( and for the diffusion problem, respectively).
As well known, whenever the domain is partitioned by considering a uniform structured mesh this latter task can be efficiently performed by means of fast Poisson solvers, among which we can list those based on the cyclic reduction idea (see e.g. [7, 8, 21]) and several specialized algebraic or geometric multigrid methods (see e.g. [9, 22, 14]). In addition, the latter can be efficiently considered also in more general mesh settings. The underlying idea is that the main effort in devising efficient algorithms must be devoted only to this simpler problem with constant coefficient, instead of the general one.
Now, the effectiveness of the proposed method is measured by referring to the optimality definition below.
 Let be a given sequence of linear systems of increasing dimensions. An iterative method is optimal if
the arithmetic cost of each iteration is at most proportional to the complexity of a matrix-vector product with matrix ,
the number of iterations for reaching the solution within a fixed accuracy can be bounded from above by a constant independent of .
In such a respect Theorem LABEL:teo:clu+se_A_esagoni proves the optimality of the PCG method: the iterations number for reaching the solution within a fixed accuracy can be bounded from above by a constant independent of the dimension and the arithmetic cost of each iteration is at most proportional to the complexity of a matrix-vector product with matrix .
Moreover, Theorem LABEL:teo:cluA proves that all the eigenvalues of the preconditioned matrix belong in a complex rectangle , with , independent of the dimension . It is worth stressing that the existence of a proper eigenvalue cluster and the aforementioned localization results in the preconditioned spectrum can be very important for fast convergence of preconditioned GMRES iterations (see, e.g., ).
Finally, we want to give notice that the PCG/PGMRES numerical
performances do not get worse in the case of unstructured meshes,
despite the lack of a rigorous proof.
4.2 Numerical results
Before analyzing in detail selected numerical results, we wish to give technical information on the performed numerical experiments. We apply the PCG or the PGMRES method, in the symmetric and non-symmetric case, respectively, with the preconditioning strategy described in section LABEL:sez:fem, to FE approximations of the problem (LABEL:eq:modello). Whenever required, the involved integrals have been approximated by means of the barycentric quadrature rule (the approximation by means of the nodal quadrature rule gives rise to similar results, indeed both are exact when applied to linear functions).
The domains of integration are those reported in Figure LABEL:fig:esagono_esagoni and we assume Dirichlet boundary conditions. All the reported numerical experiments are performed in Matlab, by employing the available pcg and gmres library functions; the iterative solvers start with zero initial guess and the stopping criterion is considered. The case of unstructured meshes is also discussed and compared, together with various types of regularity in the diffusion coefficient.
In fact, we consider the case of a coefficient function satisfying the assumptions (LABEL:eq:ipotesi_coefficienti). More precisely, the second columns in Table LABEL:tab:IT-ES123esagono_triangoliequilateri_M5 report the number of iterations required to achieve the convergence for increasing values of the coefficient matrix size when considering the FE approximation with structured uniform meshes as in Figure LABEL:fig:esagono_esagoni and with template function , .The subsequent meshes are obtained by means of a progressive refinement procedure, consisting in adding new nodes corresponding to the middle point of each edge. The numerical experiments plainly confirm the previous theoretical analysis in section LABEL:sez:clustering: the convergence behavior does not depend on the coefficient matrix dimension .
As anticipated, despite the lack of corresponding theoretical results, we want to test the convergence behavior also in the case in which the regularity assumption on in Theorems LABEL:teo:clu+se_A_esagoni and LABEL:teo:cluA are not satisfied. The analysis is motivated by favorable known numerical results in the case of FD approximations (see, for instance, [18, 19, 3]) or FE approximation with only the diffusion term . More precisely, we consider as template the function , the function , with or .
The number of required iterations is listed in the remaining columns in Table LABEL:tab:IT-ES123esagono_triangoliequilateri_M5. Notice that also in these cases with a or diffusion function, the iteration count does not depend on the coefficient matrix dimension . The same results are reported in Table LABEL:tab:IT-ES123quadrato_triangolirettangoli_M1 with respect to structured uniform meshes on the domain .
Furthermore, we want to test our proposal in the case of some unstructured meshes generated by triangle  with a progressive refinement procedure. The first meshes in the considered mesh sequences are reported in Figures LABEL:fig:Mesh_M6newnew and LABEL:fig:Mesh_M4, respectively.
Tables LABEL:tab:IT-ES123esagono_unstructured_M6newnew and LABEL:tab:IT-ES123quadrato_unstructured_M4 report the number of required iterations in the case of the previous template functions. Negligible differences in the iteration trends are observed for increasing dimensions . All these remarks are in perfect agreement with the outliers analysis of the matrices , with respect to a cluster at (see some examples in Figure LABEL:fig:outlier_analysis).
To conclude the section we take into consideration the more realistic setting in which the meshes are generated by a specialized automatic procedure. According to this point of view, we have applied a spectral approximation of the matrix sequence in terms of product of low-cost matrix structures, i.e., and , that carry the “structural” content of the variational problem and the specific “informative” content contained in the diffusion coefficient function , respectively. However, the matrix can be still not easy to handle. Therefore, we go beyond in such idea by defining a preconditioner in which the matrix related to the unstructured mesh is replaced by a suitable projection of the Toeplitz matrix , , (see Section LABEL:sez:diffusion_equation).
The numerical results are reported in Table LABEL:tab:IT_perturbed_mesh. As expected, the number of iterations is a constant with respect to the dimension when the preconditioner is applied. Nevertheless, for large enough, the same seems to be true also for the preconditioner . Indeed, for increasing dimensions , the unstructured partitioning is more and more similar to the one given by equilateral triangles sketched in Fig. LABEL:fig:esagono_esagoni.a.
A theoretical ground supporting these observation is still missing and would be worth in our opinion to be studied and developed.