Natural domain decomposition algori thms for the solution of time-harmonic elastic waves

# Natural domain decomposition algori thms for the solution of time-harmonic elastic waves

R. Brunet Department of Mathematics and Statistics, University of Strathclyde, Glasgow, UK, E-mail: Romain.Brunet@strath.ac.uk.    V. Dolean Department of Mathematics and Statistics, University of Strathclyde, Glasgow, UK, and Laboratoire J.A. Dieudonné, CNRS, University Côte d’Azur, Nice, France. E-mail: work@victoritadolean.com.    M. J. Gander Université de Genève, 2-4 rue du Lièvre, Genève. E-mail: martin.gander@unige.ch.
###### Abstract

We study for the first time Schwarz domain decomposition methods for the solution of the Navier equations modeling the propagation of elastic waves. These equations in the time harmonic regime are difficult to solve by iterative methods, even more so than the Helmholtz equation. We first prove that the classical Schwarz method is not convergent when applied to the Navier equations, and can thus not be used as an iterative solver, only as a preconditioner for a Krylov method. We then introduce more natural transmission conditions between the subdomains, and show that if the overlap is not too small, this new Schwarz method is convergent. We illustrate our results with numerical experiments, both for situations covered by our technical two subdomain analysis, and situations that go far beyond, including many subdomains, cross points, heterogeneous materials in a transmission problem, and Krylov acceleration. Our numerical results show that the Schwarz method with adapted transmission conditions leads systematically to a better solver for the Navier equations than the classical Schwarz method.

Key words. Domain decomposition methods, Schwarz preconditioners, time-harmonic elastic waves, Navier equations.

AMS subject classifications. 65N55, 65N35, 65F10

## 1 Introduction

Time harmonic problems are difficult to solve by iterative methods in the medium to high frequency regime, see [18] for the case of the Helmholtz equation, which is the prototype of such time harmonic problems with oscillatory solutions. Domain decomposition methods are a natural choice as iterative solvers for such problems, since they are by construction parallel and can still locally use direct solvers without convergence problems. To obtain good domain decomposition convergence for time harmonic problems, adapted transmission conditions are however needed between subdomains. Such transmission conditions were first studied for the Helmholtz equation by Desprès in [10, 11], and later optimized variants were introduced and analyzed by Chevalier in his PhD thesis [7], see also Chevalier and Nataf [8], the work by Collino, Delbue, Joly and Piacentini [9], and Gander et al. [24, 23, 25]. Very similar in nature to the Helmholtz equations, high-frequency time-harmonic Maxwell’s equations are also very difficult to solve iteratively, and the design of efficient domain decomposition methods for the intermediate to high frequency regime is even harder. First optimized transmission conditions both for the first and second order formulations of Maxwell’s equations can already be found in the PhD thesis of Chevalier [7, section 4.7] and Collino et al. [9], but were then more systematically developed by Alonso -Rodriguez and Gerardo-Giorda [1], and especially in Dolean et al. [14, 13, 15], see also Peng, Rawat and Lee [29] , and references therein. The Analytic Incomplete LU factorization (AILU) [19], the sweeping preconditioner [16, 17], the source transfer domain decomposition [5, 6], the method based on single layer potentials [32], and the method of polarized traces [35], are all methods in this same class of domain decomposition methods with more effective transmission conditions, which became known under the name optimized Schwarz methods, see [20, 21] for an introduction, and [26] and references therein for a thorough treatment when applied to time harmonic wave propagation problems.

To the best of our knowledge, the use of Schwarz methods for time-harmonic elastic waves modeled by the Navier equations has not been studied so far, and our goal is to investigate classical Schwarz methods, and also a new variant that uses more natural transmission conditions between the subdomains when applied to the Navier equations. To do so, we study the Schwarz methods at the continuous level, for a simplified decomposition as it has become standard with two subdomains, to gain insight into the effect of transmission conditions on the performance of the method. To test the method, we then discretize the problems and implement the Schwarz methods using Restricted Additive Schwarz (RAS) introduced by Cai and Sarkis in [4], which represents a faithful implementation of the continuous parallel Schwarz method of Lions, see [21]. This is especially important when more natural transmission conditions are used, see [31] for Optimized RAS (ORAS).

Our paper is structured as follows: in Section LABEL:sec:classical, we present and analyze the classical Schwarz algorithm applied to the Navier equations. We prove for a simplified two subdomain setting at the continuous level that the Schwarz algorithm is not a convergent iterative method in this case. We then introduce new transmission conditions in Section LABEL:sec:optimal and show first that there exist transmission conditions which make the Schwarz method converge in a finite number of steps. These transmission conditions involve however non local operators, and we thus introduce a local, low frequency approximation for the Navier equations, for which we prove convergence of the new Schwarz method provided the overlap is not too small. In Section LABEL:sec:numerical we study these new Schwarz methods numerically, first for a two subdomain decomposition covered by our analysis, but then also for the case of many subdomains with cross points and material heterogeneities. Our numerical results show that the new Schwarz method performs much better than the classical one when used as a preconditioner for a Krylov method.

## 2 Classical Schwarz algorithm for the Navier Equations

We are interested in solving the Navier equations in the frequency domain,

 −(Δe+ω2ρ)u=fin Ω, \hb@xt@.01(2.1)

where the operator is defined by . To study the basic (non)-convergence properties of the Schwarz algorithm applied to the Navier equations (LABEL:NavierEq), we consider the domain and decompose it into two unbounded overlapping subdomains and , with overlap parameter . The classical parallel Schwarz algorithm then starts with an initial guess on subdomain , , and solves for iteration index

 −(Δe+ω2ρ)un1=fin Ω1,−(Δe+ω2ρ)un2=fin Ω2,un1=un−12at x=δ,un2=un−11at x=0. \hb@xt@.01(2.2)

To study the convergence properties of this algorithm, we use a Fourier transform in the direction. We denote by the Fourier parameter and the Fourier transformed solution,

 ^u(x,k)=∫∞−∞e−ikyu(x,y)dy,u(x,y)=12π∫∞−∞eiky^u(x,k)dk.

The convergence factor for each Fourier mode of (LABEL:ClassicalSchwarz) is given in

###### Lemma 2.1 (Convergence factor of classical Schwarz)

For a given initial guess , , the classical Schwarz algorithm (LABEL:ClassicalSchwarz) with overlap multiplies at each iteration the error in each Fourier mode with the convergence factor

 ρcla(k,ω,Cp,Cs,δ)=max{|r+|,|r−|}, \hb@xt@.01(2.3)

where the eigenvalues of the iteration matrix are

 r±=X22+e−δ(λ1+λ2)±12√X2(X2+4e−δ(λ1+λ2)),X=k2+λ1λ2k2−λ1λ2(e−λ1δ−e−λ2δ)\definecolorpgfstrokecolorrgb0,0,0\pgfsys@color@gray@stroke0\pgfsys@color@gray@fill0, \hb@xt@.01(2.4)

and are given by

 λ1=√k2−ω2C2s,λ2=√k2−ω2C2p,Cp=√λ+2μρ,Cs=√μρ. \hb@xt@.01(2.5)

Proof. The convergence factor can be obtained by a direct computation working on the error equations, as it is shown in the short publication [3].

We show in Figure LABEL:Fig3a a plot of the modulus of the convergence factor (LABEL:RhoClassical) as function of the Fourier mode for an example of the parameters in the Navier equations.

We see that the classical Schwarz method converges for high frequencies, , diverges for intermediate frequencies, , and stagnates for low frequencies . We prove in the next theorem that this behavior holds for all choices of parameters in the Navier equations, and thus the classical Schwarz method is not an effective iterative solver for these equations.

###### Theorem 2.2 ((Non-) Convergence of the overlapping classical Schwarz method)

The convergence factor (LABEL:ClassicalSchwarz) of the overlapping classical Schwarz method (LABEL:ClassicalSchwarz) applied to the Navier equations (LABEL:NavierEq) satisfies

 |ρcla(k,ω,Cp,Cs,δ)|=1,k∈[0,ωCp]∪{ωCs},|ρcla(k,ω,Cp,Cs,δ)|>1,k∈(ωCp,ωCs),|ρcla(k,ω,Cp,Cs,δ)|<1,k∈(ωCs,∞), \hb@xt@.01(2.6)

where the last two results are shown to hold for overlap small.

Proof. The proof is quite technical. To simplify the notation, we define for the case when the roots in (LABEL:lambda12) are complex the quantities

 i¯λ1:=λ1=i√ω2C2s−k2,i¯λ2:=λ2=i√ω2C2p−k2. \hb@xt@.01(2.7)

We have to treat five cases: three intervals for , and two values separating the intervals: in the first interval , , and the eigenvalues (LABEL:r+r-) become

 r±=X22+e−iδ(¯λ1+¯λ2)±12√X2(X2+4e−iδ(¯λ1+¯λ2)),X=k2−¯λ1¯λ2k2+¯λ1¯λ2(e−i¯λ1δ−e−i¯λ2δ).

The square of their modulus is given by

 |r±|2=1+√A2+B2+(x2+y2)24+er(x2−y2)+2xyeiPart1±√22× \hb@xt@.01(2.8) ⎛⎜ ⎜ ⎜⎝x2−y2+2er2(√A2+B2+A)12+csgn(B−iA)(xy+ei)(√A2+B2−A)12Part2⎞⎟ ⎟ ⎟⎠,

where the complex sign is defined as

 csgn(x)={100orR(x)=0&I(x)>0,

and we introduced the quantities , , and ,

 er :=−sin(δ(¯λ1+¯λ2)),ei:=cos(δ(¯λ1+¯λ2)), x :=R(X)=k2−¯λ1¯λ2k2+¯λ1¯λ2(cos(¯λ1δ)−cos(¯λ2δ)), y

The terms and appearing in the square root are real and defined by , which gives after some computations

 A =(x2−y2)2−4x2y2−8eixy+4er(x2−y2), B =4(xy+ei)(x2−y2)+8erxy.

Then we obtain by a direct computation that

 √A2+B2=(x2+y2)√(x2+y2)2+8er(x2−y2)+16eixy+16 =16(k2−¯λ1¯λ2)2sin2(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)2 ×⎛⎝1−2(k2−¯λ1¯λ2)2sin2(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)2+(k2−¯λ1¯λ2)4sin4(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)4⎞⎠ =16sin2(δ2(¯λ1−¯λ2))((k2−¯λ1¯λ2)2cos2(δ2(¯λ1−¯λ2))+4¯λ1k2¯λ2)(k2−¯λ1¯λ2)−2(k2+¯λ1¯λ2)4,

and

 x2−y2 =−4(k2−¯λ1¯λ2)2sin2(δ2(¯λ1−¯λ2))cos(δ(¯λ1+¯λ2))(k2+¯λ1¯λ2)2, x2+y2 =4(k2−¯λ1¯λ2)2sin2(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)2, xy =2(k2−¯λ1¯λ2)2sin(δ(¯λ1+¯λ2))sin2(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)2.

We now show that in (LABEL:rcase1) vanishes identically: we get on the one hand

 (x2+y2)24=4(k2−¯λ1¯λ2)4sin4(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)4, \hb@xt@.01(2.9)

and on the other hand, we have

 √A2+B24 \hb@xt@.01(2.10) er(x2−y2) =−4(k2−¯λ1¯λ2)2sin2(δ2(¯λ1−¯λ2))cos2(δ(¯λ1+¯λ2))(k2+¯λ1¯λ2)2, 2eixy

and we obtain by adding the three terms from (LABEL:part1b) to each other

 −4(k2−¯λ1¯λ2)4sin4(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)4. \hb@xt@.01(2.11)

This leads, by adding (LABEL:part1a) and (LABEL:part1c) indeed to . We next show that also in (LABEL:rcase1) vanishes identically: we get

 x2−y22+er=cos(δ(¯λ1+¯λ2))⎛⎝1−2(k2−¯λ1¯λ2)2sin2(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)2⎞⎠, xy+ei=−sin(δ(¯λ1+¯λ2))⎛⎝1−2(k2−¯λ1¯λ2)2sin2(δ2(¯λ1−¯λ2))(k2+¯λ1¯λ2)2⎞⎠,

and for the term involving and

 √√A2+B2±A =4k2−¯λ1¯λ2(k2+¯λ1¯λ2)2sin(δ2(¯λ1−¯λ2))√1∓cos(2δ(¯λ1+¯λ2)) ×√(k2−¯λ1¯λ2)2cos2(δ2(¯λ1−¯λ2))+4k2¯λ1¯λ2.

By analyzing the signs of the different terms, we obtain for the complex sign

 csgn(B−iA)=sg(cos(δ(¯λ1+¯λ2))sin(δ(¯λ1+¯λ2))),

and after a lengthy computation we obtain

 Part2=Ck× (√1+cos(2δ(¯λ1+¯λ2))sin(δ(¯λ1+¯λ2))

where is a complicated factor depending on . A direct computation for the second factor of shows that independently of the value of , we get . We can thus conclude from (LABEL:rcase1) that and therefore the algorithm stagnates in the first interval , see the first interval in Figure LABEL:Fig3a.

At the boundary between the first and second interval, where , we have that and , and therefore the eigenvalues in (LABEL:r+r-) become

 r±=12(1+e−2i¯λ1δ)±12√(1−e−2i¯λ1δ)2,X=e−i¯λ1δ−1,

and being positive we have equivalently

 r+=1, r−=e−2i¯λ1δ⟹ρcla(ωCp,ω,Cp,Cs,δ)=max{|r+|,|r−|}=1,

and hence the algorithm stagnates also when the first interval is closed on the right, i.e. for .

In the second interval, , we have that and , and hence the eigenvalues in (LABEL:r+r-) become

 r±=X22+e−δ(i¯λ1+λ2)±12√X2(X2+4e−δ(i¯λ1+λ2)),X=k2+i¯λ1λ2k2−i¯λ1λ2(e−i¯λ1δ−e−λ2δ).

We compute the modulus of the eigenvalues and expand them for overlap parameter small to find

 |r+|=1+2ω2λ2¯λ21C2p(k4+¯λ21λ22)δ+O(δ2),|r−|=1−2ω2λ2k2C2s(k4+¯λ21λ22)δ+O(δ2).

We thus obtain that is bigger than one for small and the method diverges, see the middle interval in Figure LABEL:Fig3a111Numerically we observe that also for a large overlap, the algorithm diverges, see Figure LABEL:Fig3a, but this seems to be difficult to prove..

Between the second and third interval, where , we have that and , and hence the eigenvalues in (LABEL:r+r-) become

 r±=12(1+e−2λ2δ)±12√(1−e−2λ2δ)2.

We thus obtain

 r+=1, r−=e−2λ2δ⟹ρcla(ωCs,ω,Cp,Cs,δ)=max{|r+|,|r−|}=1,

and the algorithm stagnates for .

In the last interval, , and by expanding from (LABEL:r+r-) for small, we get

 r+=1−2λ2ω2C2s(k2−λ1λ2)δ+O(δ2)<1,r−=1−2λ1ω2C2p(k2−λ1λ2)δ+O(δ2)<1,

since . We can thus conclude that

 ρcla(k,ω,Cp,Cs,δ)=max{|r+|,|r−|}<1,

see the last interval in Figure LABEL:Fig3a, where we also see that , since all the real exponentials involved in the expressions of are decreasing to as increases.

We see from Theorem LABEL:ClassicalSchwarzTheorem that the classical Schwarz method with overlap can not be used as an iterative solver to solve the Navier equations, since the method stagnates for low frequencies and even diverges for intermediate frequencies; only high frequencies are converging. A precise estimate for how fast the classical Schwarz method diverges depending on the overlap is given in the short publication [3].

## 3 New Transmission Conditions for the Schwarz algorithm

A remedy for the divergence problems of the classical Schwarz method is to introduce different transmission conditions, and to consider the new Schwarz method

 −(Δe+ω2ρ)un1=fin Ω1,−(Δe+ω2ρ)un2=fin Ω2,(\@fontswitchT1+\@fontswitchS1)un1=(\@fontswitchT1+\@fontswitchS1)un−12x=δ,(\@fontswitchT2+\@fontswitchS2)un2=(\@fontswitchT2+\@fontswitchS2)un−11,x=0, \hb@xt@.01(3.1)

where the traction operators , , are defined by , and the operators are two by two matrix valued operators one can choose to obtain better convergence. The traction operators play for the Navier equations the role the Neumann condition plays for the Poisson equation. Like we obtained the convergence factor of the classical Schwarz algorithm using a Fourier transform in Lemma LABEL:th:convclas, we can obtain the convergence factor in the case where more general transmission operators with Fourier symbols are used.

###### Lemma 3.1

For a given initial guess , , the general Schwarz algorithm with overlap (LABEL:GeneralOptimizedSchwarz) has for each Fourier mode the convergence factor

 ρopt(k,ω,Cp,Cs,δ)=(max{|r+|,|r−|})12,r±=X22+Y±12√X2(X2+4Y), \hb@xt@.01(3.2)

with

 X=e−λ1δb11−e−λ2δb22,Y=b11b22−b12b21eλ1δeλ2δ,[b11b12b21b22]:=B−12B1, \hb@xt@.01(3.3)

where

 B1=⎡⎢ ⎢⎣ˆ\@fontswitchS2(1,1)−2λ1ρC2s−iλ1ˆ\@fontswitchS2(1,2)kˆ\@fontswitchS2(1,2)+i2k2ρC2s−λ2ˆ\@fontswitchS2(1,1)−ρω2kˆ\@fontswitchS2(2,1)−i2k2C2sρ−λ1ˆ\@fontswitchS2(2,2)−ρω2kˆ\@fontswitchS2(2,2)+2C2sρλ2+iλ2ˆ\@fontswitchS2(2,1)k⎤⎥ ⎥⎦, \hb@xt@.01(3.4)
 B2=⎡⎢ ⎢⎣ˆ\@fontswitchS2(1,1)+2λ1ρC2s+iλ1ˆ\@fontswitchS2(1,2)kˆ\@fontswitchS2(1,2)+i2k2ρC2s+λ2ˆ\@fontswitchS2(1,1)−ρω2kˆ\@fontswitchS2(2,1)−i2k2C2sρ+λ1ˆ\@fontswitchS2(2,2)−ρω2kˆ\@fontswitchS2(2,2)−2C2sρλ2−iλ2ˆ\@fontswitchS2(2,1)k⎤⎥ ⎥⎦, \hb@xt@.01(3.5)

and are given in (LABEL:lambda12).

Proof. This result is obtained by a direct calculation, replacing the solutions in Fourier space into the transmission conditions of the general Schwarz algorithm (LABEL:GeneralOptimizedSchwarz), for details, see the PhD thesis [2, Lemma 2.3].

### 3.1 An Optimal Schwarz Method

The new transmission conditions in (LABEL:GeneralOptimizedSchwarz) are a very powerful tool to fix convergence problems of the classical Schwarz method, and are used in many modern domain decomposition methods for time harmonic wave propagation, like the sweeping preconditioner, source transfer and the method of polarized traces, which are all variants of the so called optimized Schwarz methods [20, 21]; for a review, see [26]. To see how powerful this idea is, we start by introducing the best possible choice, namely transparent boundary conditions (TBC) as transmission conditions in (LABEL:GeneralOptimizedSchwarz), which leads to what is called an optimal Schwarz method222Optimal here is not used in the sense of scalability, but really means faster convergence is not possible!:

###### Theorem 3.2 (Convergence of the optimal Schwarz algorithm.)

If one chooses in the new Schwarz algorithm (LABEL:GeneralOptimizedSchwarz) the operators with the Fourier symbols

 ˆ\@fontswitchS1(1,1)=ρλ1ω2k2−λ1λ2,ˆ\@fontswitchS1(1,2)=+ikρ(2C2s−ω2k2−λ1λ2),ˆ\@fontswitchS1(2,1)=−ikρ(2C2s−ω2k2−λ1λ2),ˆ\@fontswitchS1(2,2)=ρλ2ω2k2−λ1λ2,ˆ\@fontswitchS2(1,1)=ˆ\@fontswitchS1(1,1),ˆ\@fontswitchS2(1,2)=−ˆ\@fontswitchS1(1,2),ˆ\@fontswitchS2(2,1)=−ˆ\@fontswitchS1(2,1),ˆ\@fontswitchS2(2,2)=ˆ\@fontswitchS1(2,2), \hb@xt@.01(3.6)

where and are given in (LABEL:lambda12), the resulting algorithm converges in two iterations, and this for all values of the overlap , even without overlap, .

Proof. If we replace ) defined in (LABEL:OptimalChoice) into (LABEL:OptiB1), the convergence factor obtained vanishes identically and the algorithm thus converges in two iterations, independently of any initial guess and the overlap .

To use the optimal choice (LABEL:OptimalChoice) as transmission operators in practice, one needs to back transform the associated TBC into the physical domain, and the corresponding are non local operators, because of the inverse transform with square root terms at the interfaces, like it is the case for many TBCs. It is therefore of interest to design local approximations for the optimal transmissions conditions, like in the development of absorbing boundary conditions (ABCs), which will lead to a new class of practical, so called optimized Schwarz algorithms. We approximate the symbols of the optimal transmission conditions in (LABEL:OptimalChoice) by polynomial symbols in which correspond to derivatives after the Fourier backtransform, and are thus local operators.

### 3.2 Optimized Schwarz Methods

We have seen in Section LABEL:sec:classical that the classical Schwarz method converges well for high frequency error components, large, but stagnates for low frequency components and even diverges for intermediate range frequencies, see Figure LABEL:Fig3a. It is therefore natural to approximate the operators in the transmission conditions using a low frequency expansion in the Fourier variable of the optimal choice given in Theorem LABEL:OptimalSchwarzThm. This leads to the so called Taylor transmission conditions (TTC), which have the symbols

 ˆ\@fontswitchS1(1,1)=iρωCp+iρ