Optimal preconditioning for the symmetric and non-symmetric coupling of adaptive finite elements and boundary elements

# Optimal preconditioning for the symmetric and non-symmetric 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 mesh-size of the underlying adaptively refined triangulations. Although the focus is on the non-symmetric Johnson-Nédélec one-equation coupling, the principle ideas also apply to other formulations like the symmetric FEM-BEM coupling. Numerical experiments underline our theoretical findings.

###### 2010 Mathematics Subject Classification:
65N30, 65N38, 65F08
Acknowledgment: The research of the authors MF, TF, and DP is supported through the FWF research projectÂ Adaptive Boundary Element Method, see http://www.asc.tuwien.ac.at/abem/, funded by the Austrian Science Fund (FWF) under grant P21732, as well as through the Innovative Projects Initiative of Vienna University of Technology. Moreover, MF and DP are partially funded through the FWF doctoral school Dissipation and Dispersion in Nonlinear PDEs, funded under grant W1245. This support is thankfully acknowledged.

## 1. Introduction

There exist plenty of works on preconditioning of FEM-BEM coupling equations, covering mainly the symmetric coupling with quasi-uniform 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 non-symmetric Johnson-Nédélec coupling, see e.g. [Med98], and also on preconditioning of adaptive FEM-BEM 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 non-symmetric Johnson-Né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 non-symmetry, such an approach may not work for the Johnson-Nédélec coupling in general. In [Med98], it was assumed that the coupling boundary is smooth. Hence, the double-layer 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 double-layer 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 mesh-size 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 non-symmetric stiffness matrix of the (stabilized) Johnson-Nédélec coupling, which reads in block-form

 (1) AL:=(AA−MT12M−KAV)+SST∈R(N+M)×(N+M),

see Section 2 below. Here, denotes an appropriate stabilization vector, which ensures positive definiteness of . The matrix block is the (positive semi-definite) Galerkin matrix of the FEM part, and the matrix block is the Galerkin matrix of the simple-layer integral operator . As in [FS09, MS98], we deal with block-diagonal preconditioners of the form

 (2) PL=(PA00PV).

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) dA⟨PAX,X⟩2 (3b) dV⟨PVΦ,Φ⟩2 ≤⟨AVΦ,Φ⟩2≤DV⟨PVΦ,Φ⟩2for all Φ∈RM,

where is strongly related to the FEM block of the FEM-BEM 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 mesh-ratio 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 mesh-size 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 mesh-size, 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 weakly-singular 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 FEM-BEM 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 Johnson-Nédélec coupling and define the admissible mesh-refinement 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 one-equation Bielak-MacCamy coupling. Numerical examples from Section 6 underline our theoretical predictions and conclude the work.

## 2. Johnson-Né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) −div(A∇u) =f in Ω, (4b) −Δuext =0 in Ωext:=R2∖¯¯¯¯Ω, (4c) u−uext =u0 on Γ, (4d) (A∇u−∇uext)⋅n =ϕ0 on Γ, (4e) |uext(x)| =O(1/|x|) as |x|→∞.

Here, denotes the outer normal on , and satisfies with uniform bounds on the maximal resp. minimal eigenvalue

 (5) 0

With resp. and its dual , we denote the usual Sobolev spaces on resp. . For given data , , and , it is well-known that the model problem (4) admits a unique solution , with finite energy , if we impose the compatibility condition

 (6) ⟨f,1⟩Ω+⟨ϕ0,1⟩Γ=0

to ensure the radiation condition (4e). Here, stands for the inner product, whereas denotes the extended inner product.

### 2.2. Johnson-Nédélec coupling

For the formulation of the Johnson-Nédélec coupling [JN80], the exterior solution (4b) is formulated by Green’s third formula. The latter gives rise to the simple-layer and double-layer integral operator

 (7a) V∈L(H−1/2+s(Γ);H1/2+s(Γ)), Vϕ(x):=∫ΓG(x,y)ϕ(y)dy, (7b) K∈L(H1/2+s(Γ);H1/2+s(Γ)), Kv(x):=∫Γ∂n(y)G(x,y)v(y)dy, where G(x,y):=−12πlog|x−y| denotes the fundamental solution of the 2D Laplacian and ∂n(⋅) is the normal derivative. Note that boundedness holds for all −1/2≤s≤1/2. In addition, our analysis requires the hypersingular integral operator (7c) W∈L(H1/2+s(Γ);H−1/2+s(Γ)), Wv(x):=−∂n(x)∫Γ∂n(y)G(x,y)v(y)dy,

where the integral is understood as finite part integral. It is known that and are symmetric in the sense of

 (8) ⟨ϕ,Vψ⟩Γ=⟨ψ,Vϕ⟩Γand⟨Wv,w⟩Γ=⟨Ww,v⟩Γ

for all and . Moreover, the assumption ensures -ellipticity of ,

 (9) cV∥ψ∥2H−1/2(Γ)≤⟨ψ,Vψ⟩Γfor all ψ∈H−1/2(Γ).

Finally, is semi-elliptic with kernel being the constant functions,

 (10) W1=0andcW∥v∥2H1/2(Γ)≤⟨Wv,v⟩Γfor all % v∈H1/2(Γ) with ∫Γvdx=0.

The constants depend only on .

With the definitions (7) of the layer integral operators, the model problem (4) is equivalently recast by the Johnson-Nédélec coupling: Given , , and , find such that

 (11a) ⟨A∇u,∇v⟩Ω−⟨ϕ,v⟩Γ =⟨f,v⟩Ω+⟨ϕ0,v⟩Γ, (11b) ⟨ψ,(12−K)u+Vϕ⟩Γ =⟨ψ,(12−K)u0⟩Γ

for all . Setting in (11), we see that the (non-stabilized) linear operator associated to the left-hand side of (11) is indefinite. However, let denote the contraction constant of the double-layer 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) cA≥cK/4

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 Lax-Milgram lemma and thus leads to non-symmetric, 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) ⟨L(u,ϕ),(v,ψ)⟩:=⟨A∇u,∇v⟩Ω−⟨ϕ,v⟩Γ+⟨ψ,(12−K)u+Vϕ⟩Γ+⟨1,(12−K)u+Vϕ⟩Γ⟨1,(12−K)v+Vψ⟩Γ as well as (13b) ⟨F,(v,ψ)⟩:=⟨f,v⟩Ω+⟨ϕ0,v⟩Γ+⟨ψ,(12−K)u0⟩Γ+⟨1,(12−K)u0⟩Γ⟨1,(12−K)v+Vψ⟩Γ.

Let and be arbitrary closed subspaces with . Then, the pair solves the Galerkin formulation of the Johnson-Nédélec coupling (11)

 (14a) ⟨A∇uℓ,∇vℓ⟩Ω−⟨ϕℓ,vℓ⟩Γ =⟨f,vℓ⟩Ω+⟨ϕ0,vℓ⟩Γ, (14b) ⟨ψℓ,(12−K)uℓ+Vϕℓ⟩Γ =⟨ψℓ,(12−K)u0⟩Γ

for all , if and only if it solves the operator formulation

 (15) ⟨L(uℓ,ϕℓ),(vℓ,ψℓ)⟩=⟨F,(vℓ,ψℓ)⟩

for all . The operator is non-symmetric, but linear, continuous, and elliptic, and the constants

 (16) cL:=inf0≠(u,ϕ)∈H⟨L(u,ϕ),(u,ϕ)⟩∥(u,ϕ)∥2HandCL:=sup0≠(u,ϕ)∈H∥L(u,ϕ)∥H∗∥(u,ϕ)∥H

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éa-type estimate

 (17) ∥(u,ϕ)−(uℓ,ϕℓ)∥H≤CLcLinf(vℓ,ψℓ)∈Hℓ∥(u,ϕ)−(vℓ,ψℓ)∥H,

where denotes the unique solution of the Johnson-Nédélec coupling (11).∎

### 2.3. Adaptive mesh-refinement and discrete spaces

Let be a given conforming initial triangulation of into compact and non-degenerate 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) supT∈TΩℓdiam(T)|T|1/2≤γ<∞for all ℓ∈N0,

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) diam(T)diam(T′)≤γ<∞for all T,T′∈TΓℓ with T∩T′≠∅ and all ℓ∈N0,

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 lowest-order Galerkin elements. We approximate functions by functions and functions by functions , where

 (20) Xℓ :={v∈C(Ω):v|T is affine for all T∈TΩℓ}, (21) Yℓ :={ψ∈L2(Γ):ψ|T is constant for all T∈TΓℓ}.

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) ˜NΩ0 :=NΩ0, ˜NΩℓ :=NΩℓ∖NΩℓ−1∪{z∈NΩℓ−1:ωΩℓ(z)⫋ωΩℓ−1(z)}for ℓ≥1, (22b) ˜NΓ0 :=NΓ0 ˜NΓℓ :=NΓℓ∖NΓℓ−1∪{z∈NΓℓ−1:ωΓℓ(z)⫋ωΓℓ−1(z)}for ℓ≥1.

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 hat-function with for all , where is Kronecker’s delta.

### 2.4. Galerkin system and block-diagonal 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

 AL=(AA−MT12M−KAV)+SST,

where the block matrices , , , and as well as the stabilization (column) vector are given by

 (23a) (AA)jk =⟨A∇ηℓzk,∇ηℓzj⟩Ωfor j,k=1,…,Nℓ, (23b) (AV)jk =⟨ψTj,VψTk⟩Γfor j,k=1,…,Mℓ, (23c) (M)jk =⟨ψTj,(ηℓzk)|Γ⟩Γfor j=1,…,Mℓ,k=1,…,Nℓ, (23d) (K)jk =⟨ψTj,K(ηℓzk)|Γ⟩Γfor j=1,…,Mℓ,k=1,…,Nℓ, (23e) (S)j=⟨1,(12−K)(ηℓzj)|Γ⟩Γfor j=1,…,Nℓ,(S)j+Nℓ=⟨1,VψTj⟩Γfor j=1,…,Mℓ.

We stress that as well as are sparse, whereas is dense. Note that the number of non-zeros 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) ALU=F∈RNℓ+Mℓ,

where the right-hand side vector reads

 (F)j =⟨F,(ηℓzj,0)⟩for j=1,…,Nℓ, (F)j+Nℓ =⟨F,(0,ψTj)⟩for j=1,…,Mℓ.

To formulate our block-diagonal preconditioner, we require an appropriate operator which is related to the FEM-domain 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) ⟨Au,v⟩:=⟨A∇u,∇v⟩Ω+⟨1,(12−K)u⟩Γ⟨1,(12−K)v⟩Γ. Then, the operator A:H1(Ω)→(H1(Ω))∗ is linear, symmetric, continuous, and elliptic, and the constants (25b) cA:=inf0≠u∈H1(Ω)⟨Au,u⟩∥u∥2H1(Ω)andCA:=sup0≠u∈H1(Ω)∥Au∥(H1(Ω))∗∥u∥H1(Ω)

satisfy and depend only on and from (5) as well as on .∎

In this work, we investigate block-diagonal preconditioners of the form

 (26) PL=(PA00PV),

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.

Instead of (24), we solve the preconditioned system

 (27) P−1LALU=P−1LF∈RNℓ+Mℓ.

For this non-symmetric system of linear equations, we use a preconditioned GMRES algorithm [SS86], which will be discussed in the following subsection. The preconditioned Galerkin matrix reads in block form

 (28)

### 2.5. Preconditioned GMRES algorithm

Let denote a symmetric and positive definite matrix and let denote a (possibly) non-symmetric, but positive definite matrix. Let, denote the standard unit vector with entries . The preconditioned GMRES algorithm reads as follows.

###### Algorithm 3 (Gmres).

Input: Matrices , right-hand 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 sub-matrix with entries for and compute

 Yk←argminY∈Rk∥∥R0∥PE1−¯¯¯¯¯HkY∥2.
• 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 sub-block 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 mesh-size .

For the preconditioner , we use an additive Schwarz multilevel diagonal preconditioner similar to the one in [XCH10]. Recall from (22) Define

 (29) ˜Xℓ:=span{ηℓz:z∈˜NΩℓ}⊆Xℓ.

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) P−1A:=(PLA)−1:=L∑ℓ=0˜Iℓ(˜DℓA)−1(˜Iℓ)T.

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) ζℓz∈Zℓ:={v∈C(Γ):v|T is affine for all T∈TΓℓ}

denote the boundary hat-function with for all . To construct an efficient preconditioner for the weakly-singular integral operator in 2D, we use the Haar-basis functions for all . Recall from (22). Let . Define the local subspaces

 (32) ˜Yℓ:=span{χℓz:z∈