The well-known mortar finite element discretization method (see, e.g., [DDDiscretizationSolvers, 2005Mortar, 1999Mortar, 1993MortarDD, 2001MortarMultipliers]) penalizes jumps across adjacent elements via constraints. This couples degrees of freedom across two neighboring elements and, consequentially, the mortar method does not admit the element-by-element assembly property. In contrast, that property is intrinsic to conforming finite element formulations and it is useful, e.g., in “matrix-free” computations, since it reduces the coupling across elements. Moreover, the local structure, providing the element-by-element assembly, allows for the utilization of certain element-based coarsening methods, e.g., AMGe (element-based algebraic multigrid) methods [VassilevskiMG], such that an analogous assembly structure is also maintained on coarse levels. Here, based on the simple idea in [IPpreconditioner] of introducing a dedicated space on the element interfaces, the coupling across element boundaries is removed and a convenient local structure, admitting the element-by-element assembly property, is obtained for a modified mortar formulation.
The modification used in this paper is found upon the following generic idea. The mortar formulation originally employs interface jump constraints and respective Lagrangian multipliers , leading to a jump term , on each interface between every two adjacent elements and , in the resulting Lagrangian functional. The proposed modification is to replace this term by two alternative terms, introducing an additional interface unknown and other (“one-sided”) jump constraints (leading to Lagrangian multipliers and ), associated with each . This provides , where comes from the element and – from . Thus, represents a trace of the solution on the interface space associated with . Clearly, eliminating the unknowns , recovers the original jump constraints with . An important consequence of introducing the additional space of discontinuous (from face to face) functions , associated with the interfaces , is the ability to uniquely relate each of the two new terms with one of the neighboring elements; that is, with and with . Accordingly, cross-element coupling occurs only through these new interface unknowns.
The approach exploited in this work, to construct preconditioners for a given conforming discretization (i.e., an initial formulation with no jump terms involved), is to further replace and by subdomains and (e.g., unions of fine-scale finite elements) and – by the interface between and . The subdomains, , can be viewed as forming a coarse triangulation , when each is a union of elements from an initial fine-scale triangulation . Any such subdomain, , is referred to as an agglomerated element or an agglomerate, in short. The resulting modified mortar method employs a pair of discontinuous spaces: one space of functions , associated with the agglomerates , and a second space of functions , associated with the set of interfaces, , between any two neighboring and in . This, excluding any forcing terms, leads to a resulting Lagrangian functional with local terms , associated with each , where is the local, on , version of the original symmetric positive definite (SPD) bilinear form from the given conforming discretization. The bilinear form and the respective linear algebra equations, possessing the desired local structure, of the modified mortar method are obtained, as standard, from the problem of finding a saddle point of the Lagrangian functional. This bilinear form and a reduced (Schur complement) variant of it are utilized in the construction of preconditioners for the original conforming bilinear form . Importantly, as shown in this work, while the constrained mortar-based reformulation leads to an indefinite “saddle-point problem”, the obtained preconditioners are SPD, leading to more natural analysis of their properties and the application of the conjugate gradient method. Moreover, the Schur complement from the reduced, via static condensation, form is also SPD, allowing the utilization of the abundantly available solvers and preconditioners for systems with SPD matrices, like multigrid methods.
The main contribution of the present paper is the introduction and study of a modified mortar reformulation as a technique for obtaining preconditioners for the original conforming bilinear form, utilizing the auxiliary space approach going back to S. Nepomnyaschikh (see [2007DD]) and studied in detail by J. Xu [1996AuxiliarySpace]. Both additive and multiplicative variants of the auxiliary space preconditioners are studied in combination with generic smoothers, following the abstract theory in [VassilevskiMG, Theorem 7.18] by verifying the assumptions stated there. The modified mortar form admits static condensation. Namely, the unknowns and the Lagrangian multipliers can be eliminated, using that they are decoupled from each other across elements (or subdomains) in the modified formulation, obtaining a reduced problem for the interface unknowns, , only. A further advantage of the element-by-element, or subdomain-by-subdomain, assembly property of the condensed modified mortar formulation is the applicability of the spectral AMGe approach, cf. [2012SAAMGe], for building algebraic multigrid (AMG) preconditioners for the reduced (via static condensation on the interface space) mortar bilinear form, which can be viewed as a Schur complement of the full (modified) mortar form. These auxiliary space preconditioners are implemented and their theoretically shown mesh-independent performance is tested on model 2D and 3D scalar second order elliptic problems, including high order elements.
The rest of the paper is organised as follows. Section 2 outlines basic notions, finite element spaces, notation, and a model problem of interest. The modified mortar approach is introduced in Section 3, the resulting auxiliary space preconditioners are described in Section 4, and Section 5 is devoted to the analysis of these preconditioners, showing (Theorem 5.6 and Corollary 5.7) the general optimality of a fine-scale auxiliary space preconditioning approach utilizing the mortar reformulation. Section 6 presents the reduced (Schur complement), via static condensation, form and demonstrates (Theorem 6.2) its ability to provide an optimal preconditioning strategy. Numerical results are shown in Section 7. In the end, Section 8 provides conclusions and possible future work.
This section is devoted to providing foundations. Notation and abbreviations are introduced to simplify the presentation in the rest of the paper.
2.1. Mesh and agglomeration
A domain (of dimension ) with a Lipschitz boundary, a (fine-scale) triangulation of , and a finite element space on are given. The mesh provides a set of elements and respective associated faces, where a face is the interface, of dimension , between two adjacent elements. The focus of this paper is on consisting of continuous piecewise polynomial functions, equipped with the usual nodal dofs (degrees of freedom).
Let be a partitioning of into non-overlapping connected sets of fine-scale elements, called agglomerate elements or, simply, agglomerates; see Fig. 1. In general, can be obtained by partitioning the dual graph of , which is a graph with nodes the elements in , where two elements are connected by an edge if they share a face. That is, all are described in terms of the elements . In the rest of the paper, capitalization indicates agglomerate entities in , like element (short for “agglomerate element”), face, or entity, whereas regular letters indicate fine-scale entities in , like element, face, or entity.
Viewing the elements in as collections of respective fine-scale faces, an intersection procedure over these collections constructs the agglomerate faces in as sets of faces; cf. Fig. 2. Consequently, each face can be consistently recognized as the -dimensional surface that serves as an interface between two adjacent elements in . The set of obtained faces in is denoted by . Also, the respective sets of dofs that can be associated with elements, faces, elements, and faces are available.
For additional information on agglomeration and the topology of “coarse meshes” like , see [VassilevskiMG, 2002AMGeTopology].
2.2. Nonconforming spaces
A main idea in this work is to obtain discontinuous (nonconforming) finite element spaces and formulations on and . To that purpose, define the finite element spaces and via restrictions (or traces) of functions in onto and , respectively. Namely,
Note that and are fine-scale spaces despite the utilization of agglomerate mesh structures like and . Accordingly, the bases in and are derived via respective restrictions (or traces) of the basis in . The degrees of freedom in and are obtained in a corresponding manner from the dofs in , as illustrated in Fig. 3. For simplicity, “dofs” is reserved for the degrees of freedom in , whereas “edofs” and “bdofs” are reserved for and , respectively. Moreover, “adofs” designates the edofs and bdofs collectively and is associated with the product space . In more detail, edofs and bdofs are obtained by “cloning” all respective dofs for every agglomerate entity that contains the dofs. Hence, each entity receives and it is the sole owner of a copy of all dofs it contains and there is no intersection between entities in terms of edofs and bdofs, i.e., they are completely separated without any sharing, making and
spaces of discontinuous functions. Nevertheless, dofs, edofs, and bdofs are related via their common “ancestry” founded on the above “cloning” procedure. Thus, the restrictions (or traces) of vectors or finite element functions in one of the spaces and their representations as vectors or functions in some of the other spaces is seen and performed in purely “algebraic” context. For example, the meaning of, where is a vector in terms of the edofs of some and , is natural as a vector in terms of the bdofs of . This is unambiguous and should lead to no confusion as it only involves a subvector and an appropriate index mapping. In what follows, finite element functions are identified with vectors on the degrees of freedom in the respective spaces.
The portions of vectors corresponding to edofs and bdofs are respectively indexed by “” and “”, leading to the notation, for , , where and . By further splitting , it is obtained , where “” denotes the edofs in the interiors of all and “” are the edofs that can be mapped to some bdofs, based on the previously-described procedure of “cloning” dofs into edofs and bdofs. Analogously, the splitting , for , is introduced in terms of dofs, indexed “”, in the interiors of all and dofs, indexed “”, related to bdofs; see Fig. 3. Note that there is a clear difference between “” and “”, but also a clear relation. Locally, on , “” and “” can, in fact, be equated, but in a global setting, it is necessary to distinguish between “” and “” indices. This should not cause any ambiguity below.
Coarse subspaces and
can be constructed by, respectively, selecting linearly independent vectors, for every , and , for every , forming the bases for the coarse spaces. This is an “algebraic” procedure, formulating the coarse basis functions as linear combinations of fine-level basis functions, i.e., as vectors in terms of the fine-level degrees of freedom. The basis vectors are organized appropriately as columns of prolongation (or interpolation) matrices and , forming as . For consistency, the corresponding degrees of freedom (associated with the respective coarse basis vectors) in and are respectively called “edofs” and “bdofs”, and the collective term “adofs” is associated with . In essence, this is based on the ideas in AMGe (element-based algebraic multigrid) methods [VassilevskiMG, 2002AMGeTopology, 2003AMGe, 2007AMGe, 2008AMGe, 2012SAAMGe].
2.3. Model problem
The model problem considered in this paper is the second order scalar elliptic partial differential equation (PDE)
where , , is a given permeability field, is a given source, and is the unknown function. For simplicity of exposition, the boundary condition on , the boundary of , is considered, i.e., . The ubiquitous variational formulation,
of (2.1) is utilized, providing the weak form
where denotes the inner products in both and , and , for . Consider the fine-scale finite element space defined on . Using the finite element basis in , (2.3) induces the following linear system of algebraic equations:
for the global SPD stiffness matrix . Moreover, the local, on agglomerates, symmetric positive semidefinite (SPSD) stiffness matrices , for , are obtainable, such that (s.t.) (the summation involves an implicit local-to-global mapping).
3. Constrained mortar-based formulation
The preconditioners, proposed in this paper, are based on the ideas of a mortar method [DDDiscretizationSolvers, DDTheoryAlgorithms]. The space pair with its degrees of freedom (adofs) is employed, together with the ability to construct subspace pairs by selecting basis functions as vectors expressed in terms of the adofs in ; see the end of Section 2.2. Having the spaces determined, a mortar-based approach is introduced here, where the jumps across interfaces of elements are penalized via equality constraints. This provides a generic modified mortar reformulation of the original problem, which is utilized as a preconditioner in an auxiliary space framework. In this section, the mortar-based formulation is presented and discussed. In the sections that follow, it is further applied to construct auxiliary space preconditioners for (2.4), their properties are addressed, and a block-preconditioning technique for the constrained mortar-type problem, based on static condensation, is described.
Using the prolongation operators defined in the end of Section 2.2 and the local (on ) versions, , of the (fine-scale) matrix in (2.4), consider the discrete nonconforming constrained quadratic minimization reformulation of (2.2)
for . Here, are the local, on , versions of in (2.4), and are the -orthogonal projections onto the local, on , spaces spanned by the vectors , associated with and constituting the basis of , where is the restriction of the diagonal, , of the global , in (2.4), onto the bdofs of . The problem is posed directly on a subspace of , requiring the explicit use of the prolongators and , for generality and to avoid over-constraining the formulation. Assume that contains the local, on faces, constants and , for each , is a proper subspace of the trace space, on , of the functions in . Thus, the constraints provide that the jumps vanish only in a subspace on each face.
Formulation (3.1) induces the, respectively global and local, on , (modified) mortar matrices
Here, is associated with the bilinear form , where are respectively test and trial vectors. Also, represents
where is a test Lagrangian multiplier vector for the constraint in (3.1), while is associated with
where is a trial vector for the interface traces. Here, denotes the local, on , version of , i.e., it is the matrix with for columns. The Lagrangian multipliers are associated with the pairs of elements and corresponding faces as they enforce equalities involving “one-sided” traces. Observe that , , and are simply sub-matrices of , , and , respectively, since the assembly of from involves only copying without any summation. This is due to the fact that all coupling is through the constraints on the faces, which is represented by the off-diagonal blocks in , and no other connections across elements and faces exist. Clearly, and are generally symmetric indefinite matrices, where and are SPSD, while is square and SPD, and has a full column rank.
Furthermore, denote the leading block sub-matrices of and in (3.2), respectively, by
where, as noted above, . The following basic, but useful, result is obtained.
The matrices and in (3.3) are invertible and the block of is symmetric negative semidefinite (SNSD).
As long as (3.1) is not over-constrained (provided by the spaces utilized here), has a full column rank. Thus, any nonzero vector in the null space of , involves a nonzero local , i.e., for some , such that , on all , and , where denotes the local, on , version of , i.e., it is the matrix with for columns. Thus, is singular if and only if has a nonzero null vector in with vanishing -projections on all faces, which is not the case here, since the null space of is spanned by the constant vector and contains the piecewise constants. That is, and do not share a common nonzero null vector (); cf. [2005SaddleProblems]. Hence, and are invertible. Finally, owing to [2005SaddleProblems, formula (3.8)], it holds that the block of is SNSD. Particularly, if (which is always the case when ), the block of is SND (symmetric negative definite). ∎
In view of the full rank of , in (3.2) is invertible if and only if .
For simplicity, the argument in Lemma 3.1 makes use of the properties of the particular model problem (2.1). Namely, it utilizes that the possible null space of is spanned by a constant, which cannot vanish on any portion of the boundary of the element. In general, the result in Lemma 3.1 hold whenever cannot possess nonzero null-space vectors that vanish on the respective faces. This is always the case when considering problems coming from PDEs, for which Dirichlet-type boundary conditions on portions of the boundary lead to nonsingular problems.
where the notation in (3.3) is used. Observe that is obtainable by performing only local, on , computations and can be assembled element by element from . This is an important property (employed in Section 6), resulting from the utilization of interface spaces, like and (Section 2.2), and the particular formulation (3.1). Now, it is not difficult to establish the following corollary.
It holds that is SPSD, is SPD, and in (3.2) is invertible.
The SPSD property of follows immediately from (3.4) and the SNSD property of the block of in Lemma 3.1. Counting on the presence of essential boundary conditions for the PDE (2.1), at least one (cf. (2.4)) is nonsingular. Hence, in view of the full rank of , at least one is SPD and the assembly property provides that is SPD. Finally, the invertibility of and implies the invertibility of . ∎
4. Auxiliary space preconditioners
The application of the (modified) mortar formulation for building auxiliary space preconditioners is addressed now.
Let be a linear transfer operator, defined in detail below. Identifying finite element functions with vectors allows to view as a matrix and obtain . In order to define the action of , recall that the relation between dofs, on one side, and edofs, on the other, is known and unambiguous. Therefore, it is reasonable to define the action of as taking the arithmetic average, formulated in terms of edofs that correspond to a particular dof, of the entries of a given vector in and obtaining the respective entries of a mapped vector defined on dofs. Namely, for any dof , let be the set of corresponding edofs (respective “cloned” degrees of freedom) and consider a vector , defined in terms of edofs. Then,
Clearly, all row sums of equal 1 and each column of has exactly one nonzero entry. Assuming that is a regular (non-degenerate) mesh [BrennerFEM], is bounded, independently of . That is,
for a constant , which depends only on the regularity of , but not on . Indeed, let be the global maximum number of elements, in , that a dof can belong to, which, in turn, is bounded by the global maximum number of elements, in , that a dof can belong to.
Define the transfer operator as
where the linear mapping takes the form .
Let be a “smoother” for in (2.4), such that is SPD, and – a symmetric (generally, indefinite) preconditioner for . Define the additive auxiliary space preconditioner for
and the multiplicative auxiliary space preconditioner for
where is the symmetrized (in fact, SPD) version of . In case is symmetric, in can be replaced by . The action of is obtained via a standard “two-level” procedure:
Given , is computed by the following steps:
residual transfer to auxiliary space: ;
auxiliary space correction: ;
correction transfer from auxiliary space: ;
A particular smoother , for , which is a part of the auxiliary space preconditioners employed in this paper, is shortly described now. Particularly, a polynomial smoother, based on the Chebyshev polynomial of the first kind, is utilized. For a given integer , consider the polynomial of degree on
satisfying , where is the Chebyshev polynomial of the first kind on . Then, is defined as
Equivalently, , where is a parameter satisfying , for all , and is the diagonal of (or another appropriate diagonal matrix). Note that is SPD and the action of such a polynomial smoother is computed via Jacobi-type iterations, using the roots of the polynomial, which makes it convenient for parallel computations; see [MSthesis, Section 4.2.2]. In practice, in (4.5) can be replaced by a diagonal weighted -smoother like , where . Such a choice allows setting . More information on the subject can be found in [VassilevskiMG, 2011Smoothers, 2012SAAMGe, 2012ConvSAAMG, 1999TGElasticity, 2012xinv].
Properties of the preconditioners of the type in (4.3) and (4.4) are studied next, showing their optimality in a fine-scale setting. Consider the auxiliary space preconditioners, involving the exact inversion of ,
The seemingly apparent auxiliary space here is . However, owing to (3.1) and the structure of in (4.2), a subspace of (more precisely, a subspace of ) can be more accurately viewed as the auxiliary space and the properties of and depend on the properties of the “subactions” of and on that subspace.
Clearly, can be eliminated from (3.1) by replacing the constraints with , for all , where denote the respective elements adjacent to . These altered conditions pose direct constraints on the jumps across faces of . The modified minimization problem is equivalent to (3.1) in the sense that it has an identical set of minimizers in . The constraints can be implicitly imposed by employing the constrained subspace defined as
Obviously, any solution to (3.1) is in . Consequently, the unconstrained quadratic minimization
is SPSD. Here, for the convenience of maintaining consistent vector notation, the basis of and its edofs are employed to represent functions in the constrained subspace . The invertibility of (Corollary 3.4) in (3.2) implies the existence of a unique minimizer in of (3.1) and, equivalently, of (5.2), which in turn implies that is SPD on .
Notice that the matrix associated with , with respect to the full basis of , is in (3.2). That is,
The discussion above demonstrates that, while is SPSD on the entire , it is SPD, when restricted to the constrained subspace , i.e., for all , since the constraints filter out its nonzero null vectors (). Hence, the “inversion operator” is well-defined as follows: for any , is the unique function (vector) that satisfies
That is, when , provides the “least-squares solution” associated with a minimization of the type in (5.2). Then, owing to (4.2) and the equivalence between (3.1) and (5.2), it holds that . Thus, is SPSD and the following proposition is obtained.
The preconditioners and in (5.1) are SPD.
Therefore, the desired “subactions” of , , and , are respectively provided by , , and on the auxiliary space , equipped with the norm induced by , i.e., the -norm. In this context, the preconditioners in (5.1) can be expressed as
Particularly, when , i.e., no coarsening of is employed and , then and is SPD.
For the rest of this section, the fine-scale case is studied, where the respective is consistently defined as
i.e., the (coarse) interface subspace, , is still employed for the jumps. The analysis follows a similar pattern to [IPpreconditioner]. Define the operator , for , via
for each edof , where , for the corresponding dof . This describes a procedure that appropriately copies the entries of , so that the respective finite element functions, corresponding to and , can be viewed as coinciding in . That is, is an injection (embedding) of into . As matrices, has the fill-in pattern of , with all nonzero entries replaced by 1.
Clearly, , the identity on , implying that is surjective, i.e., it has a full row rank. Consequently, for any , (exactly) approximates in the sense
This is to be expected, since, in a sense, “includes” .
Showing the continuity of in terms of the respective “energy” norms is more challenging and requires to satisfy certain properties. Using the indexing notation in Section 2.2, introduce the following splittings, for ,