Optimal preconditioning for the symmetric and nonsymmetric coupling of adaptive finite elements and boundary elements
Abstract.
We analyze a multilevel diagonal additive Schwarz preconditioner for the adaptive coupling of FEM and BEM for a linear 2D Laplace transmission problem. We rigorously prove that the condition number of the preconditioned system stays uniformly bounded, independently of the refinement level and the local meshsize of the underlying adaptively refined triangulations. Although the focus is on the nonsymmetric JohnsonNédélec oneequation coupling, the principle ideas also apply to other formulations like the symmetric FEMBEM coupling. Numerical experiments underline our theoretical findings.
Key words and phrases:
FEMBEM coupling, preconditioner, multilevel additive Schwarz, adaptivity2010 Mathematics Subject Classification:
65N30, 65N38, 65F081. Introduction
There exist plenty of works on preconditioning of FEMBEM coupling equations, covering mainly the symmetric coupling with quasiuniform meshes, see [CKL98, FS09, HPPS03, HMS99, HS98, KS02, MS98] and the references therein. In contrast to that, only little is known on preconditioning of the nonsymmetric JohnsonNédélec coupling, see e.g. [Med98], and also on preconditioning of adaptive FEMBEM couplings. It is the main goal of this paper to close this gap and to extend the existing analysis to the case of the (adaptive) symmetric as well as nonsymmetric JohnsonNédélec coupling [JN80]. For the symmetric coupling [Cos88, Han90], the approach of Bramble & Pasciak [BP88] applies, which guarantees positive definiteness and symmetry of the Galerkin matrix with respect to a special inner product [CKL98, HPPS03, KS02]. Therefore, efficient iterative solvers designed for symmetric and positive definite matrices are applicable. However, due to the nonsymmetry, such an approach may not work for the JohnsonNédélec coupling in general. In [Med98], it was assumed that the coupling boundary is smooth. Hence, the doublelayer integral operator is compact. The system matrix can therefore be split into a symmetric part plus a compact perturbation part (the Galerkin matrix of the doublelayer integral operator). Preconditioning is done only on the symmetric part with the theory of [BP88], and convergence results for iterative solvers can then be obtained by compact perturbation theory assuming that the meshsize of the coarsest mesh is sufficiently small. In general, however, the coupling boundary is not smooth. Therefore, the preconditioner theory must not rely on compactness of .
To state the contributions of the current work, we consider the nonsymmetric stiffness matrix of the (stabilized) JohnsonNédélec coupling, which reads in blockform
(1) 
see Section 2 below. Here, denotes an appropriate stabilization vector, which ensures positive definiteness of . The matrix block is the (positive semidefinite) Galerkin matrix of the FEM part, and the matrix block is the Galerkin matrix of the simplelayer integral operator . As in [FS09, MS98], we deal with blockdiagonal preconditioners of the form
(2) 
Here, the appropriate operator induces a coercive, symmetric, and bounded bilinear form . The symmetric and positive definite matrices resp. are spectrally equivalent to the symmetric and positive definite Galerkin matrices resp. , i.e.
(3a)  
(3b) 
where is strongly related to the FEM block of the FEMBEM system (1), see (23)–(25) below. Inspired by [MS98], we prove that the condition number of as well as the number of iterations to reduce the relative residual by a factor in the preconditioned GMRES algorithm with inner product depends only on .
Usually, the condition number of Galerkin matrices and on adaptively refined meshes hinges on the global meshratio as well as on the number of degrees of freedom. Therefore, the construction of optimal preconditioners for iterative solvers is a necessary task. Here, optimality is understood in the sense that the condition number of the preconditioned matrix is independent of the meshsize and the degrees of freedom.
Very recently, it was proven in [XCH10] for 2D FEM with energy space that a local multilevel additive Schwarz preconditioner is optimal, i.e. the constants in (3) are independent of the meshsize, the number of degrees of freedom, and the number of levels. Here, “local” means, that scaling at each level is done only on newly created nodes plus neighbouring nodes, where the associated basis functions have changed. An analogous result for 2D and 3D hypersingular integral equations with energy space has been derived by the authors [FFPS13]. In [FFPS13, XCH10], the proofs rely on a stable space decomposition of the discrete subspaces in resp. . Alternatively, [XCN09] provides stable subspace decompositions in for higher order elements in any dimension on bisection grids.
In the present work, we prove the optimality of some local multilevel additive Schwarz preconditioner for 2D weaklysingular integral equations with energy space . The proof is derived by postprocessing of the corresponding result for the hypersingular integral equation [FFPS13]. Combining this with the result of [XCH10], we prove optimality of for the FEMBEM coupling.
Notation. Throughout the work, we explicitly state all constants and their dependencies in all statements of results. In proofs, however, we use the notation to abbreviate with a constant which is clear from the context. Moreover, abbreviates . Furthermore, the entries of a vector or a matrix are denoted by resp. . By , we denote the duality brackets between a Hilbert space and its dual .
Outline. The remainder of this work is organized as follows: In Section 2, we recall the basic facts on the JohnsonNédélec coupling and define the admissible meshrefinement strategies. We also formulate the preconditioned GMRES Algorithm 3, which is required to state the main result (Theorem 4). In Section 3, we prove the spectral estimates (3). Section 4 adapts the analysis for the symmetric coupling [MS98] and contains the proof of the main result (Theorem 4). The short Section 5 deals with extensions of the developed theory to the symmetric coupling and the oneequation BielakMacCamy coupling. Numerical examples from Section 6 underline our theoretical predictions and conclude the work.
2. JohnsonNédélec coupling and main result
2.1. Model problem and analytical setting
Let be a polygonal and simply connected domain with boundary . We consider the following Laplace transmission problem in free space:
(4a)  
(4b)  
(4c)  
(4d)  
(4e) 
Here, denotes the outer normal on , and satisfies with uniform bounds on the maximal resp. minimal eigenvalue
(5) 
With resp. and its dual , we denote the usual Sobolev spaces on resp. . For given data , , and , it is wellknown that the model problem (4) admits a unique solution , with finite energy , if we impose the compatibility condition
(6) 
to ensure the radiation condition (4e). Here, stands for the inner product, whereas denotes the extended inner product.
2.2. JohnsonNédélec coupling
For the formulation of the JohnsonNédélec coupling [JN80], the exterior solution (4b) is formulated by Green’s third formula. The latter gives rise to the simplelayer and doublelayer integral operator
(7a)  
(7b)  
where denotes the fundamental solution of the 2D Laplacian and is the normal derivative. Note that boundedness holds for all . In addition, our analysis requires the hypersingular integral operator  
(7c) 
where the integral is understood as finite part integral. It is known that and are symmetric in the sense of
(8) 
for all and . Moreover, the assumption ensures ellipticity of ,
(9) 
Finally, is semielliptic with kernel being the constant functions,
(10) 
The constants depend only on .
With the definitions (7) of the layer integral operators, the model problem (4) is equivalently recast by the JohnsonNédélec coupling: Given , , and , find such that
(11a)  
(11b) 
for all . Setting in (11), we see that the (nonstabilized) linear operator associated to the lefthand side of (11) is indefinite. However, let denote the contraction constant of the doublelayer potential [SW01]. Following the analysis in [OS13, Say09, Ste11], also (stabilized) Galerkin formulations of (11) admit unique solutions, if the ellipticity constant from (5) satisfies
(12) 
and if the equation is either explicitly stabilized [OS13, Ste11] or if the discrete subspace of contains the constant functions [Say09]. In [AFF13a], the result of [Say09] is reproduced with a new proof. Introducing the notion of implicit stabilization, an equivalent elliptic operator equation of (11) is derived which fits in the frame of the LaxMilgram lemma and thus leads to nonsymmetric, but positive definite Galerkin matrices. The main result of [AFF13a] reads as follows (and also holds for a strongly monotone, but nonlinear material tensor ):
Lemma 1.
Let (12) be satisfied. For , define
(13a)  
as well as  
(13b) 
Let and be arbitrary closed subspaces with . Then, the pair solves the Galerkin formulation of the JohnsonNédélec coupling (11)
(14a)  
(14b) 
for all , if and only if it solves the operator formulation
(15) 
for all . The operator is nonsymmetric, but linear, continuous, and elliptic, and the constants
(16) 
satisfy and depend only on and from (5) as well as on . Moreover, is a continuous linear functional on . In particular, (15) (and hence also (14)) admits a unique solution , and there holds the Céatype estimate
(17) 
where denotes the unique solution of the JohnsonNédélec coupling (11).∎
2.3. Adaptive meshrefinement and discrete spaces
Let be a given conforming initial triangulation of into compact and nondegenerate triangles. We suppose that a sequence of refined triangulations is obtained by newest vertex bisection, see Figure 1, where is the coarsest conforming mesh such that all marked elements have been bisected. For our analysis, the set of marked elements is arbitrary, but in practice obtained from local a posteriori refinement indicators, see e.g. [AFF13a]. We note that newest vertex bisection guarantees uniform shape regularity in the sense that
(18) 
where depends only on the initial mesh , see e.g. [Ver13, KPP13] and the references therein.
Let be a given initial partition of the coupling boundary into compact line segments. We suppose that a sequence of refined partitions is obtained by bisection, where the refined elements are refined into two sons of half length, i.e. with and where at least the marked elements are refined, i.e. . In addition, we suppose that the meshes are uniformly shape regular in the sense that
(19) 
where depends only on the initial partition . Possible choices include the 1D bisection algorithms from [AFF13b]. A further choice is to consider the partition of which is induced by the triangulation of . Formally, such a coupling of and is not required for the analysis, but simplifies the implementation and is therefore used in the numerical experiments of Section 6.
In this work, we consider lowestorder Galerkin elements. We approximate functions by functions and functions by functions , where
(20)  
(21) 
Let denote the set of nodes of the triangulation , and let denote the set of nodes of the triangulation . For resp. , we define the patch resp. . For the construction of optimal multilevel preconditioners on adaptively refined triangulations, we need the following subsets of resp. :
(22a)  
(22b) 
The sets resp. thus consist of the new nodes plus the old nodes, where the corresponding patches have changed. A visualization of is given in Figure 2. For each node , denotes the associated hatfunction with for all , where is Kronecker’s delta.
2.4. Galerkin system and blockdiagonal preconditioning
Let with denote the nodal basis of , and let with denote a basis of , where denotes the characteristic function on . The Galerkin matrix of the operator has the form
where the block matrices , , , and as well as the stabilization (column) vector are given by
(23a)  
(23b)  
(23c)  
(23d)  
(23e) 
We stress that as well as are sparse, whereas is dense. Note that the number of nonzeros in the matrix is bounded by . Moreover, the application of the rank stabilization matrix can be implemented efficiently with complexity for use with an iterative solver.
The discrete variational formulation (15) is equivalent to solving the following linear system of equations: Find such that
(24) 
where the righthand side vector reads
To formulate our blockdiagonal preconditioner, we require an appropriate operator which is related to the FEMdomain part of . The next lemma follows from a Rellich compactness argument, since for all constants . Details are analogous to, e.g., [AFF13a, Lemma 10] and therefore left to the reader.
Lemma 2.
For , define
(25a)  
Then, the operator is linear, symmetric, continuous, and elliptic, and the constants  
(25b) 
satisfy and depend only on and from (5) as well as on .∎
In this work, we investigate blockdiagonal preconditioners of the form
(26) 
where is a “good” approximation of the Galerkin matrix corresponding to the operator from Lemma 2 with respect to the nodal basis of , and is a “good” approximation of the Galerkin matrix . Our construction below ensures that and hence are symmetric and positive definite.
2.5. Preconditioned GMRES algorithm
Let denote a symmetric and positive definite matrix and let denote a (possibly) nonsymmetric, but positive definite matrix. Let, denote the standard unit vector with entries . The preconditioned GMRES algorithm reads as follows.
Algorithm 3 (Gmres).
Input: Matrices , righthand side vector , initial guess , relative tolerance , and maximum number of iterations with .

Allocate memory for the matrix , the vectors , , , , and .

Compute initial residual and .

Set counter , and initialize for all .
Iterate the following steps (i)–(vii):

Compute .

For all compute

Compute .

Define the submatrix with entries for and compute

Compute and .

Stop iteration if or .

Otherwise, compute , update counter , and goto (iv).
Output: , , and .
Note that for being the identity matrix, Algorithm 3 is the usual GMRES algorithm with inner product , see e.g. [SS86]. The main memory consumption is given by the vectors , while the matrix in step (iv) is a subblock of the matrix and thus does not need to be stored explicitly.
As is often the case for multilevel preconditioners, the application of on a vector is known, whereas the application of is unknown. We therefore note that the GMRES Algorithm 3 can be implemented without using to compute the inner products and the norms . To this end, one replaces the computation of by , where . Then, . In step (i), we compute instead of . Step (ii) is replaced by and . Instead of step (iii), we then compute . Note that in step (vi). Lastly, in step (vii) we replace the computation of by and .
2.6. Local multilevel preconditioner and main result
For both the FEM part and BEM part in (26), we will use local multilevel preconditioners which are optimal in the sense that the condition numbers of the preconditioned systems are independent of the number of levels and the meshsize .
For the preconditioner , we use an additive Schwarz multilevel diagonal preconditioner similar to the one in [XCH10]. Recall from (22) Define
(29) 
Let denote the canonical embedding with matrix representation and . Furthermore, let denote the diagonal of the Galerkin matrix with respect to the local set of nodes , i.e. for and . Then, our local multilevel diagonal preconditioner is defined via
(30) 
From the definition, we see that this preconditioner corresponds to a diagonal scaling on each level, where scaling is done on the local subset only.
For all boundary nodes , let
(31) 
denote the boundary hatfunction with for all . To construct an efficient preconditioner for the weaklysingular integral operator in 2D, we use the Haarbasis functions for all . Recall from (22). Let . Define the local subspaces
(32) 