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

TL;DR

This paper introduces and analyzes new Schwarz domain decomposition algorithms with improved transmission conditions for solving time-harmonic elastic wave equations, demonstrating enhanced convergence and efficiency over classical methods.

Contribution

The paper proves classical Schwarz methods do not converge for Navier equations and proposes new transmission conditions that ensure convergence and better performance.

Findings

Classical Schwarz method is not convergent for Navier equations.

New Schwarz method with adapted transmission conditions converges when overlap is sufficient.

Numerical experiments show improved solver efficiency with the new method.

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.

Tables

Table 1. Table 1: Physical characteristics for the heterogeneous transmission problem.

Domain	E	$\nu$	$\rho$	$\mu$	$\lambda$	$C_p$	$C_s$	f	$\omega$
	${2.10}^{11}$	0.3	7800	${77.10}^9$	${12.10}^{10}$	5927	3142	${10}^4$	$2\pi {10}^4$
$0.5\le r\le 1$	${2.10}^{11}$	0.47	7800	${68.10}^9$	${11.10}^{11}$	12588	2952	${10}^4$	$2\pi {10}^4$

Equations

$$-\left(\Delta^{e}+\omega^{2}\rho\right)\mathbf{u}=\mathbf{f}\quad\mbox{in $\Omega$},$$

$$\begin{array}[]{rcllrcll}-\left(\Delta^{e}+\omega^{2}\rho\right)\mathbf{u}_{1}^{n}&=&\mathbf{f}&\mbox{in $\Omega_{1}$,}&-\left(\Delta^{e}+\omega^{2}\rho\right)\mathbf{u}_{2}^{n}&=&\mathbf{f}&\mbox{in $\Omega_{2}$,}\\ \mathbf{u}_{1}^{n}&=&\mathbf{u}_{2}^{n-1}&\mbox{at $x=\delta$},&\mathbf{u}_{2}^{n}&=&\mathbf{u}_{1}^{n-1}&\mbox{at $x=0$}.\end{array}$$

$$\hat{\mathbf{u}}(x,k)={\int }_{-\infty }^{\infty }{e}^{-\mathrm{i}ky}\mathbf{u}(x,y)\mathrm{d}y,\mathbf{u}(x,y)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{\mathrm{i}ky}\hat{\mathbf{u}}(x,k)\mathrm{d}k.$$

$${\rho }_{cla}\left(k,\omega ,C_p,C_s,\delta \right)=\max \{|r_{+}|,|r_{-}|\},$$

$$r_{\pm}=\frac{X^{2}}{2}+e^{-\delta(\lambda_{1}+\lambda_{2})}\pm\frac{1}{2}\sqrt{X^{2}\left(X^{2}+4e^{-\delta(\lambda_{1}+\lambda_{2})}\right)},\,X=\frac{k^{2}+\lambda_{1}\lambda_{2}}{k^{2}-\lambda_{1}\lambda_{2}}\left(e^{-\lambda_{1}\delta}-e^{-\lambda_{2}\delta}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}$$

$${\lambda }_1=\sqrt{k^2-\frac{{\omega }^2}{C_s^2}},{\lambda }_2=\sqrt{k^2-\frac{{\omega }^2}{C_p^2}},C_p=\sqrt{\frac{\lambda +2\mu }{\rho }},C_s=\sqrt{\frac{\mu }{\rho }}.$$

$$\begin{array}[]{rcll}|\rho_{cla}(k,\omega,C_{p},C_{s},\delta)|&=&1,&k\in[0,\frac{\omega}{C_{p}}]\cup\{\frac{\omega}{C_{s}}\},\\ |\rho_{cla}(k,\omega,C_{p},C_{s},\delta)|&>&1,&k\in(\frac{\omega}{C_{p}},\frac{\omega}{C_{s}}),\\ |\rho_{cla}(k,\omega,C_{p},C_{s},\delta)|&<&1,&k\in(\frac{\omega}{C_{s}},\infty),\end{array}$$

$$\mathrm{i}{\bar{\lambda }}_1:={\lambda }_1=\mathrm{i}\sqrt{\frac{{\omega }^2}{C_s^2}-k^2},\mathrm{i}{\bar{\lambda }}_2:={\lambda }_2=\mathrm{i}\sqrt{\frac{{\omega }^2}{C_p^2}-k^2}.$$

$$r_{\pm }=\frac{X^2}{2}+e^{-\mathrm{i}\delta ({\bar{\lambda }}_1+{\bar{\lambda }}_2)}\pm \frac{1}{2}\sqrt{X^2\left(X^2+4e^{-\mathrm{i}\delta \left({\bar{\lambda }}_1+{\bar{\lambda }}_2\right)}\right)},X=\frac{k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2}{k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2}\left(e^{-\mathrm{i}{\bar{\lambda }}_1\delta }-e^{-\mathrm{i}{\bar{\lambda }}_2\delta }\right).$$

$$|r_{\pm }{|}^2=1+\underset{{\mathrm{Part}}_1}{\underbrace{\frac{\sqrt{A^2+B^2}+(x^2+y^2{)}^2}{4}+e_r\left(x^2-y^2\right)+2xy{e}_i}}\pm \frac{\sqrt{2}}{2}\times$$

$$\left(\underset{{\mathrm{Part}}_2}{\underbrace{\frac{x^2-y^2+2e_r}{2}{\left(\sqrt{A^2+B^2}+A\right)}^{\frac{1}{2}}+\mathrm{csgn}\left(B-\mathrm{i}A\right)(xy+e_i){\left(\sqrt{A^2+B^2}-A\right)}^{\frac{1}{2}}}}\right),$$

$$\mathrm{csgn}(x)=\{\begin{array}{rl}1& 0<\mathfrak{R}(x)\text{mbox}or\mathfrak{R}(x)=0\&0<\mathfrak{I}(x),\\ -1& \mathfrak{R}(x)>0\text{mbox}or\mathfrak{R}(x)=0\&\mathfrak{I}(x)>0,\end{array}$$

$$e_r$$

$$x$$

$$y$$

$$A$$

$$B$$

$$\sqrt{A^2+B^2}=\left(x^2+y^2\right)\sqrt{{\left(x^2+y^2\right)}^2+8e_r\left(x^2-y^2\right)+16e_{i}xy+16}$$

$$=\frac{16{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2{\sin }^2\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)}{{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2}$$

$$\times \left(1-\frac{2{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2{\sin }^2\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)}{{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2}+\frac{{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^4{\sin }^4\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)}{{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^4}\right)$$

$$=\frac{16{\sin }^2\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)\left({\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2{\cos }^2\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)+4{\bar{\lambda }}_1{k}^2{\bar{\lambda }}_2\right)}{{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^{-2}{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^4},$$

$$x^2-y^2$$

$$x^2+y^2$$

$$xy$$

$$\frac{{\left(x^2+y^2\right)}^2}{4}=4\frac{{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^4{\sin }^4\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)}{{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^4},$$

$$\frac{\sqrt{A^2+B^2}}{4}$$

$$e_r\left(x^2-y^2\right)$$

$$2e_{i}xy$$

$$-\frac{4\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2{)}^4{\sin }^4\left(\frac{\delta }{2}({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)}{{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^4}.$$

$$\frac{x^2-y^2}{2}+e_r=\cos \left(\delta \left({\bar{\lambda }}_1+{\bar{\lambda }}_2\right)\right)\left(1-2\frac{{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2{\sin }^2\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)}{{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2}\right),$$

$$xy+e_i=-\sin \left(\delta \left({\bar{\lambda }}_1+{\bar{\lambda }}_2\right)\right)\left(1-2\frac{{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2{\sin }^2\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)}{{\left(k^2+{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2}\right),$$

$$\sqrt{\sqrt{A^2+B^2}\pm A}$$

$$\times \sqrt{{\left(k^2-{\bar{\lambda }}_1{\bar{\lambda }}_2\right)}^2{\cos }^2\left(\frac{\delta }{2}\left({\bar{\lambda }}_1-{\bar{\lambda }}_2\right)\right)+4k^2{\bar{\lambda }}_1{\bar{\lambda }}_2}.$$

$$\mathrm{csgn}(B-\mathrm{i}A)=\mathrm{sg}\left(\cos \left(\delta \left({\bar{\lambda }}_1+{\bar{\lambda }}_2\right)\right)\sin \left(\delta \left({\bar{\lambda }}_1+{\bar{\lambda }}_2\right)\right)\right),$$

$${\mathrm{Part}}_2=C_k\times$$

$$-\mathrm{csgn}\left(B-\mathrm{i}A\right)\cos \left(\delta \left({\bar{\lambda }}_1+{\bar{\lambda }}_2\right)\right)\sqrt{1-\cos \left(2\delta \left({\bar{\lambda }}_1+{\bar{\lambda }}_2\right)\right)}),$$

$$r_{\pm }=\frac{1}{2}(1+e^{-2\mathrm{i}{\bar{\lambda }}_1\delta })\pm \frac{1}{2}\sqrt{{\left(1-e^{-2\mathrm{i}{\bar{\lambda }}_1\delta }\right)}^2},X=e^{-\mathrm{i}{\bar{\lambda }}_1\delta }-1,$$

Full text

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 2, 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 3 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 4 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,

[TABLE]

where the operator $\Delta^{e}$ is defined by $\Delta^{e}\mathbf{u}=\mu\Delta\mathbf{u}+(\lambda+\mu)\nabla(\nabla\cdot\mathbf{u})$. To study the basic (non)-convergence properties of the Schwarz algorithm applied to the Navier equations (1), we consider the domain $\Omega:={\mathbb{R}}^{2}$ and decompose it into two unbounded overlapping subdomains $\Omega_{1}:=(-\infty,\delta)\times{\mathbb{R}}$ and $\Omega_{2}:=(0,\infty)\times{\mathbb{R}}$, with overlap parameter $\delta>0$. The classical parallel Schwarz algorithm then starts with an initial guess $\mathbf{u}_{j}^{0}$ on subdomain $\Omega_{j}$, $j=1,2$, and solves for iteration index $n=1,2,\ldots$

[TABLE]

To study the convergence properties of this algorithm, we use a Fourier transform in the $y$ direction. We denote by $k\in\mathbb{R}$ the Fourier parameter and $\hat{u}(x,k)$ the Fourier transformed solution,

[TABLE]

The convergence factor for each Fourier mode of (2) is given in

Lemma 1** (Convergence factor of classical Schwarz).**

For a given initial guess $\mathbf{u}_{j}^{0}\in(L^{2}(\Omega_{j}))^{2}$, $j=1,2$, the classical Schwarz algorithm (2) with overlap $\delta>0$ multiplies at each iteration the error in each Fourier mode $k$ with the convergence factor

[TABLE]

where the eigenvalues of the iteration matrix are

[TABLE]

and $\lambda_{1,2}\in\mathbb{C}$ are given by

[TABLE]

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 1 a plot of the modulus of the convergence factor (3) as function of the Fourier mode $k$ for an example of the parameters in the Navier equations.

We see that the classical Schwarz method converges for high frequencies, $|\rho_{cla}|<1$, diverges for intermediate frequencies, $|\rho_{cla}|>1$, and stagnates for low frequencies $|\rho_{cla}|=1$. 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** ((Non-) Convergence of the overlapping classical Schwarz method).**

The convergence factor (2) of the overlapping classical Schwarz method (2) applied to the Navier equations (1) satisfies

[TABLE]

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

Proof.

The proof is quite technical. To simplify the notation, we define for the case when the roots $\lambda_{1,2}$ in (5) are complex the quantities

[TABLE]

We have to treat five cases: three intervals for $k$, and two values $k\in\{\frac{\omega}{C_{p}},\frac{\omega}{C_{s}}\}$ separating the intervals: in the first interval $k\in(0,\frac{\omega}{C_{p}})$, $\lambda_{1,2}\in\mathrm{i}\mathbb{R_{+}}$, and the eigenvalues (4) become

[TABLE]

The square of their modulus is given by

[TABLE]

where the complex sign is defined as

[TABLE]

and we introduced the quantities $e_{r}$, $e_{i}$, $x$ and $y$,

[TABLE]

The terms $A$ and $B$ appearing in the square root are real and defined by $A+\mathrm{i}B:=X^{2}\left(X^{2}+4\operatorname{e}^{-\mathrm{i}\delta\left(\bar{\lambda}_{1}+\bar{\lambda}_{2}\right)}\right)$, which gives after some computations

[TABLE]

Then we obtain by a direct computation that

[TABLE]

and

[TABLE]

We now show that $\operatorname{Part}_{1}$ in (8) vanishes identically: we get on the one hand

[TABLE]

and on the other hand, we have

[TABLE]

and we obtain by adding the three terms from (10) to each other

[TABLE]

This leads, by adding (9) and (11) indeed to $\operatorname{Part}_{1}\equiv 0$. We next show that also $\operatorname{Part}_{2}$ in (8) vanishes identically: we get

[TABLE]

and for the term involving $A$ and $B$

[TABLE]

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

[TABLE]

and after a lengthy computation we obtain

[TABLE]

where C_{k}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{R^{*}}:=\mathbb{R}\backslash\left\{0\right\}} is a complicated factor depending on $k$. A direct computation for the second factor of $\operatorname{Part}_{2}$ shows that independently of the value of $\operatorname{csgn}(B-\mathrm{i}A)$, we get $\operatorname{Part}_{2}\equiv 0$. We can thus conclude from (8) that $\rho_{cla}\left(k,\omega,C_{p},C_{s},\delta\right)=\max\{|r_{+}|,|r_{-}|\}=|r_{+}|=|r_{-}|=1$ and therefore the algorithm stagnates in the first interval $k\in[0,\frac{\omega}{C_{p}})$, see the first interval in Figure 1.

At the boundary between the first and second interval, where $k=\frac{\omega}{C_{p}}$, we have that $\lambda_{2}=0$ and $\lambda_{1}\in\mathrm{i}\mathbb{R^{*}_{+}}$, and therefore the eigenvalues in (4) become

[TABLE]

and $\mathfrak{Re}(1-\operatorname{e}^{-2\mathrm{i}\bar{\lambda}_{1}\delta})=1-\cos(2\bar{\lambda}_{1}\delta)$ being positive we have equivalently

[TABLE]

and hence the algorithm stagnates also when the first interval is closed on the right, i.e. for $k\in[0,\frac{\omega}{C_{p}}]$.

In the second interval, $k\in(\frac{\omega}{C_{p}},\frac{\omega}{Cs})$, we have that $\lambda_{1}\in\mathrm{i}\mathbb{R^{*}_{+}}$ and $\lambda_{2}\in\mathbb{R^{*}_{+}}$, and hence the eigenvalues in (4) become

[TABLE]

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

[TABLE]

We thus obtain that $\rho_{cla}(k,\omega,C_{p},C_{s},\delta)=\max\{|r_{+}|,|r_{-}|\}$ is bigger than one for $\delta$ small and the method diverges, see the middle interval in Figure 1111Numerically we observe that also for a large overlap, the algorithm diverges, see Figure 1, but this seems to be difficult to prove..

Between the second and third interval, where $k=\frac{\omega}{C_{s}}$, we have that $\lambda_{1}=0$ and $\lambda_{2}=\frac{\omega\sqrt{C_{p}^{2}-C_{s}^{2}}}{C_{s}C_{p}}>0$, and hence the eigenvalues in (4) become

[TABLE]

We thus obtain

[TABLE]

and the algorithm stagnates for $k=\frac{\omega}{C_{s}}$.

In the last interval, $k\in\left(\frac{\omega}{Cs},\infty\right)$, $\lambda_{1,2}\in\mathbb{R^{*}_{+}}$ and by expanding $r_{\pm}>0$ from (4) for $\delta$ small, we get

[TABLE]

since $k^{2}-\lambda_{1}\lambda_{2}>0$. We can thus conclude that

[TABLE]

see the last interval in Figure 1, where we also see that $\lim_{k\rightarrow\infty}r_{\pm}=0$, since all the real exponentials involved in the expressions of $r_{\pm}$ are decreasing to [math] as $k$ increases. ∎

We see from Theorem 2 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

[TABLE]

where the traction operators ${\cal T}_{j}$, $j=1,2$, are defined by $T_{j}(\mathbf{u})=2\mu\frac{\partial\mathbf{u}}{\partial n_{j}}+\lambda\mathbf{n}_{j}\nabla\cdot\mathbf{u}+\mu\mathbf{n}_{j}\times\nabla\times\mathbf{u}$, and the operators ${\cal S}_{j}$ are two by two matrix valued operators one can choose to obtain better convergence. The traction operators ${\cal T}_{j}$ 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 1, we can obtain the convergence factor in the case where more general transmission operators ${\cal S}_{1,2}$ with Fourier symbols $\widehat{\cal S}_{1,2}$ are used.

Lemma 3**.**

For a given initial guess $\mathbf{u}_{j}^{0}\in(L^{2}(\Omega_{j}))^{2}$, $j=1,2$, the general Schwarz algorithm with overlap (12) has for each Fourier mode the convergence factor

[TABLE]

with

[TABLE]

where

[TABLE]

[TABLE]

and $\lambda_{1,2}\in\mathbb{C}$ are given in (5).

Proof.

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

3.1 An Optimal Schwarz Method

The new transmission conditions in (12) 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 (12), 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 4** (Convergence of the optimal Schwarz algorithm.).**

If one chooses in the new Schwarz algorithm (12) the operators ${\cal S}_{j}$ with the Fourier symbols

[TABLE]

where $\lambda_{1}$ and $\lambda_{2}$ are given in (5), the resulting algorithm converges in two iterations, and this for all values of the overlap $\delta\geq 0$, even without overlap, $\delta=0$.

Proof.

If we replace $(\widehat{\cal S}_{1},\widehat{\cal S}_{2}$) defined in (17) into (15), the convergence factor obtained vanishes identically and the algorithm thus converges in two iterations, independently of any initial guess and the overlap $\delta\geq 0$. ∎

To use the optimal choice (17) as transmission operators in practice, one needs to back transform the associated TBC into the physical domain, and the corresponding ${\cal S}_{j}$ 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 (17) by polynomial symbols in $\mathrm{i}k$ which correspond to derivatives after the Fourier backtransform, and are thus local operators.

3.2 Optimized Schwarz Methods

We have seen in Section 2 that the classical Schwarz method converges well for high frequency error components, $k$ large, but stagnates for low frequency components and even diverges for intermediate range frequencies, see Figure 1. It is therefore natural to approximate the operators ${\cal S}_{j}$ in the transmission conditions using a low frequency expansion in the Fourier variable $k$ of the optimal choice given in Theorem 4. This leads to the so called Taylor transmission conditions (TTC), which have the symbols

[TABLE]

and $\widehat{\cal S}_{2}$ with the same relation to $\widehat{\cal S}_{1}$ as for the optimal choice in Theorem 4. A zeroth order approximation would thus be

[TABLE]

which was also obtained as an ABC using a different argument in [33]. These ABCs happen to be exact for a particular combination of plane waves, and thus have a physical sense for this particular problem.

We show in Figure 2 the modulus of the convergence factor of the optimized Schwarz method with Taylor transmission conditions.

We see that the method now converges very well for low frequencies, and also for intermediate frequencies. For high frequencies, we see that without overlap, $\delta=0$, the method stagnates, since the convergence factor equals 1. Increasing the overlap leads to convergence of the very high frequencies, and when the overlap becomes big enough, the method seems to converge for all frequencies, except at the two points $k\in\{\frac{\omega}{C_{p}},\frac{\omega}{C_{s}}\}$. This is a very important improvement compared to the classical Schwarz method, see Figure 1, and while for Helmholtz equations there is one non-convergent frequency when using optimized transmission conditions [24, 23, 25], for the Navier equations there are two. We prove in the following theorem that the numerical observations in Figure 2 indeed hold for all parameter choices in the Navier equations in the non-overlapping case.

Theorem 5** (Convergence of the non-overlapping Schwarz algorithm with TTC).**

The new Schwarz method (12) with TTC (18) for non-overlapping decompositions converges for $k\in(0,\frac{\omega}{C_{s}})\backslash\{\frac{\omega}{C_{p}}\}$, and stagnates with the contraction factor being equal to $1$ for $k\in[\frac{\omega}{C_{s}},\infty)$.

Proof.

The proof is again quite technical: the eigenvalues of the iteration matrix are given by

[TABLE]

where the elements in the matrix $B$ are given by

[TABLE]

and $Z_{1}$, $Z_{2}$, $K$ and $D$ are defined by

[TABLE]

We define now $\bar{\lambda}_{j}\in\mathbb{R_{+}},\,j=1,2$ as in (7), and study the five cases for $k$ as in the proof of Theorem 2: if $k\in(0,\frac{\omega}{C_{p}})$ then $\lambda_{1,2}\in\mathrm{i}\mathbb{R_{+}}$, and using (20) we obtain

[TABLE]

A direct computation shows that $X^{2}+2Y>0$ and $X^{2}+4Y>0$, and hence $r_{+}>|r_{-}|>0$, so we just need to check that

[TABLE]

To show this second inequality, we compute

[TABLE]

and the last inequality can be checked by first setting $X=\widetilde{X}/{{D}}$ and $Y=\widetilde{Y}/{{D}}^{2}$, which leads to the condition

[TABLE]

where $C\in\mathbb{R^{*}}$ is a complicated factor depending on $C_{p}$, $C_{s}$, $\omega$, and $k$, and the other terms are positive. We thus conclude that in this case the algorithm is convergent.

If $k=\frac{\omega}{C_{p}}$ then $\lambda_{1}=\mathrm{i}\frac{\omega\sqrt{C_{p}^{2}-C_{s}^{2}}}{C_{s}C_{p}}$ and $\lambda_{2}=0$, and the elements in the matrix $B$ are

[TABLE]

and the eigenvalues $r_{\pm}$ are given by

[TABLE]

Since $C_{p}^{3}-4C_{p}C_{s}^{2}+4C_{s}^{3}>0$, we have $|r_{-}|<1$, and thus $\rho_{T_{0}}=1$.

If $k\in(\frac{\omega}{C_{p}},\frac{\omega}{C_{s}})$ then $\lambda_{1}\in\mathrm{i}\mathbb{R_{+}}$ and $\lambda_{2}\in\mathbb{R_{+}}$, and we obtain

[TABLE]

By computing their modulus, we get

[TABLE]

where

[TABLE]

An upper bound $M$ for the modulus of the eigenvalues is thus obtained choosing the plus sign,

[TABLE]

and it suffices to prove that $M<1$. To do so, it is sufficient to show that for the numerator in the first term of $M$, we have

[TABLE]

and for the numerator in the second term of $M$, we have

[TABLE]

By a direct computation, one can show that both (23) and (24) are equivalent to

[TABLE]

which clearly holds, and thus $\max\left(|r_{+}|,|r_{-}|\right)\leq M<1$ and the algorithm is convergent.

If $k=\frac{\omega}{C_{s}}$ then $\lambda_{1}=0$ and $\lambda_{2}=\frac{\omega\sqrt{C_{p}^{2}-C_{s}^{2}}}{C_{s}C_{p}}>0$. In this case the coefficients of the matrix $B$ are given by

[TABLE]

and the eigenvalues $r_{\pm}$ are

[TABLE]

and the algorithm therefore stagnates for $k=\frac{\omega}{C_{s}}$.

If $k\in\left(\frac{\omega}{Cs},\infty\right)$ then $\lambda_{1,2}\in\mathbb{R^{*}_{+}}$ and (20) gives $r_{\pm}=\dfrac{1}{D}(R\pm\mathrm{i}I)$ with

[TABLE]

A direct computation shows that for $C(\omega,k,C_{p},C_{s})\in\mathbb{R}^{*}$ a constant

[TABLE]

since $K^{2}\lambda_{1}\lambda_{2}-4Z_{1}Z_{2}+4\omega^{3}\lambda_{1}\lambda_{2}\frac{C_{p}}{C_{s}}=8\lambda_{1}\lambda_{2}C_{s}^{3}\left(4C_{s}^{3}k^{2}+C_{p}\omega^{2}\right)\left(k^{2}+\lambda_{1}^{2}\right)(k^{2}-\lambda_{1}^{2}-\frac{\omega^{2}}{C_{s}^{2}})$, and $k^{2}-\lambda_{1}^{2}-\frac{\omega^{2}}{C_{s}^{2}}=0$, and hence $|r_{\pm}|=1$ and the algorithm stagnates. ∎

The non-overlapping Schwarz algorithm with Taylor transmission conditions thus leads to good convergence for low frequencies, but stagnates for high frequencies. We now investigate if the combination of overlap and TTC can lead to a convergent optimized Schwarz algorithm. A first result for strictly positive overlap $\delta>0$ is the following, see also Figure 2 for an illustration:

Theorem 6** (Convergence of the overlapping Schwarz algorithm with TTC.).**

For $\delta>0$ small, the new overlapping Schwarz method (12) with Taylor transmission conditions (18) converges for

[TABLE]

diverges for $k\in(\frac{\omega}{C_{s}},k^{*})$, and stagnates for $k\in\{\frac{\omega}{C_{p}},\frac{\omega}{C_{s}},k^{*}\}$.

Proof.

Again the proof is quite technical: the eigenvalues of the iteration matrix are

[TABLE]

where the elements of the matrix $B$ are

[TABLE]

and $Z_{1}$, $Z_{2}$, $K$ and $D$ are given by

[TABLE]

We define $\bar{\lambda}_{j}\in\mathbb{R_{+}}$, $j=1,2$, as in (7) when ${\lambda}_{1}$ and/or ${\lambda}_{2}\in i\mathbb{R}$. When the overlap $\delta$ is small, a series expansion of the eigenvalues gives

[TABLE]

and the modulus of the eigenvalues becomes

[TABLE]

Again we need to distinguish several cases: if $k\in(0,\frac{\omega}{C_{p}})$ then $\lambda_{1,2}\in\mathrm{i}\mathbb{R_{+}}$ and $I_{1\pm}=R_{2\pm}=0$ for both eigenvalues. Therefore the series expansion (27) becomes

[TABLE]

where

[TABLE]

After simplifications, this gives exactly the same convergence factor as in the non-overlapping case for which we have proved in Theorem 5 that it is less than one. Therefore the algorithm is convergent in this case for $\delta>0$ small enough.

If $k=\frac{\omega}{C_{p}}$ then $\lambda_{1}=\mathrm{i}\frac{\omega\sqrt{C_{p}^{2}-C_{s}^{2}}}{C_{s}C_{p}}$ and $\lambda_{2}=0$. In this case the elements of the matrix $B$ are

[TABLE]

and the eigenvalues $r_{\pm}$ are

[TABLE]

Since $C_{p}^{3}-4C_{p}C_{s}^{2}+4C_{s}^{3}>0$, we have $|r_{-}|<1$, and thus $\rho_{T_{0}}=1$ which means the algorithm stagnates in this case.

If $k\in(\frac{\omega}{C_{p}},\frac{\omega}{C_{s}})$, then $\lambda_{1}\in\mathrm{i}\mathbb{R_{+}}$ and $\lambda_{2}\in\mathbb{R_{+}}$. The series expansion (27) becomes

[TABLE]

and the terms $(R_{1\pm}+\mathrm{i}I_{1\pm})$ are the same as in the non-overlapping case, and we already know from the proof of Theorem 5 that $\left(R_{1\pm}^{2}+I_{1\pm}^{2}\right)<1$. Therefore the algorithm is convergent in this case for $\delta>0$ small enough333From Figure 2 we see that actually the overlap makes the algorithm faster in this interval, and even slightly faster also in the first interval..

If $k=\frac{\omega}{C_{s}}$ then $\lambda_{1}=0,\,\lambda_{2}=\frac{\omega\sqrt{C_{p}^{2}-C_{s}^{2}}}{C_{s}C_{p}}>0.$ In this case the elements in the matrix $B$ are

[TABLE]

and the eigenvalues $r_{\pm}$ of the iteration matrix are given by

[TABLE]

which shows that the algorithm stagnates.

Finally, if $k\in(\frac{\omega}{C_{s}},\infty)$, then $\lambda_{1,2}\in\mathbb{R^{*}_{+}}$ and the eigenvalues are given by (26). We then use for $\delta>0$ small the series expansion (27) for $r_{\pm}$ and obtain

[TABLE]

where the values ($R,I,D$) are given in (25) and (22). Hence $R_{1\pm}^{2}+I_{1\pm}^{2}=1$ and

[TABLE]

since $(\lambda_{1}-\lambda_{2})<0<(Z_{1}+Z_{2})$. As the first eigenvalue is less than one ,

[TABLE]

we will focus now on $r_{-}\sim 1+R_{1-}R_{2-}+I_{1-}I_{2-}=:F(k)$, with

[TABLE]

Note that $\sqrt{(Z_{1}+Z_{2})^{2}-{\lambda}_{2}\lambda_{1}K^{2}}\in\mathbb{R}$ since we have

[TABLE]

which holds because the discriminant $K^{2}-4(4C_{s}^{3}k^{2}+C_{p}\omega^{2})Z_{1}=-\frac{4C_{p}\omega^{6}}{C_{s}}<0$. So we do not have real solutions and the dominant term being positive, we conclude this inequality holds for all $k>\frac{\omega}{C_{s}}$. We can conclude that $g(k)\in\mathbb{R}$ as we have seen previously. We now need to investigate under which conditions $g(k)<0$ which is equivalent to $r_{-}>1$. By a direct calculation, we obtain

[TABLE]

We next study the sign of $\tilde{g}$ in a neighborhood of $\frac{\omega}{C_{s}}$: we set $k=\frac{\omega}{C_{s}}+\varepsilon$, and then expand $\tilde{g}$ in a series for $\varepsilon$ small, which leads to

[TABLE]

For sufficiently small values of $\varepsilon$, that is for $k$ close to $\frac{\omega}{C_{s}}$, the leading term of this series being negative, we have $r_{-}>1$ for $\delta>0$ small enough and the algorithm diverges. On the other hand, because of the overlap, $\lim\limits_{k\to\infty}\rho_{T_{0}}(k,\omega,C_{p},C_{s},\delta)=0$ and by continuity there exist two values $k^{*}>\bar{k}>\frac{\omega}{C_{s}}$ such that for all $k>k^{*}$ we have ${\rho_{T_{0}}}(k)<1$, at $k=k^{*}$ we have $\rho_{T_{0}}(k^{*},\omega,C_{p},C_{s},\delta)=1$, and $|\rho_{T_{0}}(\bar{k})|>1$ with $\bar{k}=\mbox{argmax}_{k>\frac{\omega}{C_{s}}}|\rho_{T_{0}}|\in(\frac{\omega}{C_{s}},k^{*})$, which concludes the proof for small overlap $\delta>0$. ∎

It is possible to obtain an asymptotic estimate for $k^{*}$ and also the rate at which the method diverges for the frequency $k=\bar{k}$, see the PhD thesis [2, pp. 45 ff]. We focus however next on how to obtain a convergent algorithm. The results of Theorem 6 hold for overlap $\delta>0$ small enough: if the overlap is bigger, it is possible to obtain a convergent optimized Schwarz method except for the two isolated frequencies $k=\frac{\omega}{C_{p}}$ and $k=\frac{\omega}{C_{s}}$, as indicated in Figure 2 for $\delta=0.8$, where the bump in the convergence factor making it larger than one has disappeared. In the Helmholtz case, there is also one isolated frequency which is not convergent when using an optimized Schwarz method [24, 23, 25], and such isolated cases can be left to Krylov acceleration. We are therefore interested in estimating the value $\delta^{*}(C_{p},C_{s},\omega)$ for which the optimized Schwarz method with Taylor transmission conditions converges as soon as the overlap $\delta>\delta^{*}(C_{p},C_{s},\omega)$ like illustrated in Figure 3,

where we see with a zoom that $\delta=0.8$ is not quite enough for convergence, but $\delta=1$ is.

Theorem 7**.**

The new overlapping Schwarz algorithm (12) with Taylor transmission conditions (18) converges for $k\in\mathbb{R}_{+}\>\backslash\>\left\{\frac{\omega}{C_{p}},\frac{\omega}{C_{s}}\right\}$ if the overlap $\delta$ is bigger than

[TABLE]

where $\alpha$ is the positive root of

[TABLE]

Proof.

As illustrated in Figure 3, we need to investigate how the convergent algorithm turns into a divergent one when $\delta$ is decreased. For $k\in(0,\frac{\omega}{C_{s}})\backslash\{\frac{\omega}{C_{p}}\}$, the Schwarz algorithm with absorbing boundary conditions converges both without overlap (see Theorem 5) and with a small overlap (see Theorem 6), and a bigger overlap only improves the behavior, so divergence does not happen for those values of $k$. If $k\in\{\frac{\omega}{C_{p}},\frac{\omega}{C_{s}}\}$, we know that the convergence factor is independent of the size $\delta$ of the overlap and always equals $1$, so the algorithm stagnates there. Only if $k\in(\frac{\omega}{C_{s}},\infty)$, the algorithm could diverge, and we thus need to study the slope of the eigenvalues of the iteration matrix at $\frac{\omega}{Cs}$ coming from the right , see Figure 3. To do so, we set $k:=\frac{\omega}{C_{s}}+\varepsilon$ for $\varepsilon$ a parameter and expand $r_{\pm}$ in a series as in (27) for $\varepsilon$ small, with $\left(R_{j\pm},I_{j\pm}\right)\in\mathbb{R}$, and then obtain for the modulus of the eigenvalues

[TABLE]

For $r_{+}$, we obtain for the first term that

[TABLE]

and similarly for $r_{-}$ we get $R_{1-}+\mathrm{i}I_{1-}=1\;\Longrightarrow\;R_{1-}^{2}+I_{1-}^{2}=1$. For the second term, we get

[TABLE]

from which we can conclude that $|{r_{+}}|^{2}<1$. For the second eigenvalue, we get however

[TABLE]

We therefore need to study the function $f$ to investigate for which values of $\delta$ it is becoming negative, which means the algorithm will diverge. Computing the derivative, we obtain

[TABLE]

so the sign of $f^{\prime}$ is the opposite sign of $g$ given by

[TABLE]

Computing the derivative of $g$, we find

[TABLE]

and we have $g(0)=-8C_{s}(C_{p}+C_{s})(C_{p}^{2}-2C_{s}^{2})$. This shows that $g^{\prime}(\delta)>0$ for all $\delta>0$ and $g(0)<0$, since $C_{p}$ and $C_{s}$ are positive and $C_{p}^{2}>2C_{s}^{2}$, see (5). Now $\cosh$ is a strictly increasing function for positive arguments, and in our case all the parameters are real and positive, and for $\delta=0$ we have $\cosh(0)=1$. We therefore have

[TABLE]

and can thus estimate $g$ from below,

[TABLE]

Let $\bar{\delta}\in\mathbb{R}_{+}^{*}$ be the unique value of $\delta$ such that $\cosh\left(\frac{\omega}{C_{s}C_{p}}\sqrt{C_{p}^{2}-C_{s}^{2}}\bar{\delta}\right)=3$; then we obtain from (31) the lower bound

[TABLE]

Since $g(0)<0$ there exists by continuity a $\hat{\delta}\in(0,\bar{\delta})$ s.t. $g(\hat{\delta})=0$ and we know that $g$ is an increasing function. This implies, because $f(0)=0$ that $f$ is a strictly increasing function for $\delta<\hat{\delta}$, and a strictly decreasing function for $\delta>\hat{\delta}$, and by a direct calculation, we find for the second derivative

[TABLE]

therefore $\hat{\delta}$ is the absolute maximum for $f$. Since $\lim_{\delta\to\infty}f(\delta)=-\infty$, its graph will cut the x-axis only once. By solving the equation $f(\delta)=0$ w.r.t. $\delta$ we find

[TABLE]

where $\alpha$ the positive root of

[TABLE]

Note that $\alpha=0$ is also a solution but since ${\delta}^{*}>0$ we must have $\alpha>0$. ∎

4 Numerical results

We discretize the Navier equations by a finite element method using a triangulation $\mathcal{T}_{h}$ of the computational domain $\Omega$, and obtain a linear system $A\boldsymbol{u}=\boldsymbol{b}$ to solve. To present the discretized Schwarz method s, let $\{\mathcal{T}_{h,i}\}_{i=1}^{N}$ be a non-overlapping partition of the triangulation $\mathcal{T}_{h}$, obtained by using a mesh partitioner like METIS [28]. The overlapping partition needed in the Schwarz methods is defined as follows : for an integer value $l\geq 0$, we build the decomposition $\{\mathcal{T}_{h,i}^{l}\}_{i=1}^{N}$ such that $\mathcal{T}_{h,i}^{l}$ is the set of all triangles from $\mathcal{T}_{h,i}^{l-1}$ and all triangles from $\mathcal{T}_{h}\setminus\mathcal{T}_{h,i}^{l-1}$ that have non-empty intersection with $\mathcal{T}_{h,i}^{l-1}$, and $\mathcal{T}_{h,i}^{0}=\mathcal{T}_{h,i}$. With this definition the width of the overlap is $2l$ mesh layers. We denote by $W_{h}$ the finite element space associated with $\mathcal{T}_{h}$, and by $W_{h,i}^{l}$ the local finite element spaces on $\mathcal{T}_{h,i}^{l}$, which form a triangulation of $\Omega_{i}$. Let $\mathcal{N}$ be the set of indices of degrees of freedom of the global finite element space $W_{h}$ and $\mathcal{N}_{i}^{l}$ the set of indices of degrees of freedom of the local finite element spaces $W_{h,i}^{l}$ for $l\geq 0$. We define the restriction operators from the global set of degrees of freedom to the local one by ${R}_{i}:W_{h}\rightarrow W_{h,i}^{l}$. At the discrete level this is a rectangular matrix $|\mathcal{N}_{i}^{l}|\times|\mathcal{N}|$ containing zeros and ones such that if $\boldsymbol{v}$ is the vector of degrees of freedom of $v_{h}\in W_{h}$, then ${R}_{i}\boldsymbol{v}$ is the vector of degrees of freedom of $W_{h}$ in $\Omega_{i}$. The extension operator from $W_{h,i}^{l}$ to $W_{h}$ and its associated matrix are then given by $R_{i}^{T}$. In addition we introduce a partition of unity $D_{i}$ as a diagonal matrix $|\mathcal{N}_{i}^{l}|\times|\mathcal{N}_{i}^{l}|$ such that

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}I}\in\mathbb{R}^{|\mathcal{N}|\times|\mathcal{N}|} is the identity matrix. With these ingredients we can now present the Restricted Additive Schwarz (RAS) preconditioner as described in [12, Chapter 1.4] ,

[TABLE]

In our experiments we will also use the Optimized RAS (ORAS) preconditioner which is based on local boundary value problem s with absorbing boundary conditions. In this case, let ${B}_{i}$ be the matrix associated to a discretization of the corresponding local problems on the subdomains $\Omega_{i}$ with absorbing boundary conditions on $\partial\Omega_{i}\cap\partial\Omega_{j}$. The definition of the preconditioner is then very similar to (33) except that $R_{i}AR_{i}^{T}$ is replaced by ${B_{i}}$ ,

[TABLE]

It has been shown in [21] that the discretized parallel Schwarz algorithm is equivalent to the stationary iteration

[TABLE]

where the preconditioner $M^{-1}$ can either be $M^{-1}_{RAS}$ from (33) or $M^{-1}_{ORAS}$ from (34) ; see [31] for the precise result for the latter which contains an algebraic condition. For more information on the influence of the partition of unity, see [22].

The stationary iteration (35) can be accelerated using a Krylov method, which is equivalent to solving the preconditioned system

[TABLE]

using the Krylov method , see e.g. [12, Chapter 3]. We test our new Schwarz methods both as stationary iterations and as preconditioners for a Krylov method. In all the following test cases, we use as stopping criteri on the relative $L^{2}$ norm of the error,

[TABLE]

where $\boldsymbol{u}$ is the mono-domain solution and $\boldsymbol{u}_{m}$ denotes the approximation of $\boldsymbol{u}$ at the $m$-th iteration of the iterative solver. Note that when using Krylov acceleration, we can also use the relative residual to stop the iteration, which is also available when the solution $\boldsymbol{u}$ is not known.

We use a zero initial guess444 When studying optimized parameters, starting with a zero initial guess is not advisable, see [21, end of subsection 5.1]. in all our tests, and we vary the size of the overlap and the type of the decomposition (uniform or using METIS). Numerical simulations were done using the open source software Freefem++ [27], which is a high level language for the variational discretization of partial differential equations.

4.1 Two-subdomain case: optimized Schwarz with TTC

We first illustrate Theorem 6 which states that the optimized Schwarz algorithm with Taylor transmission conditions can have converge problems for frequencies $k$ slightly bigger than $\frac{\omega}{C_{s}}$ if the overlap is not big enough. We use the parameters $C_{p}=1$, $C_{s}=\frac{1}{2}$, $\rho=1$, the domain $\Omega=(-1,1)\times(0,1)$ with Dirichlet conditions on top and bottom, and absorbing boundary conditions on the left and right, and the two subdomains \Omega_{1}=(-1,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\delta})\times(0,1) and \Omega_{2}=(-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\delta},1)\times(0,1). We discretize the time-harmonic Navier equations using uniform P1 finite elements with mesh size $h=\frac{1}{80}$. We show in Figure 4

on the left the error in modulus at iteration $60$ of the optimized Schwarz method with Taylor transmission conditions for $\omega=5$ and overlap parameter $\delta=2h$. We see that the optimized Schwarz method stops converging: the interval for convergence problems predicted by Theorem 6 is $\left[\frac{\omega}{C_{s}},k^{\star}\right]=[10,k^{\star}]$ , and we observe that the error on the left in Figure 4 has $5$ bumps along the interface which corresponds well to the mode $|\sin(ky)|$ along the interface for k=5\pi\approx 15{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}>\frac{\omega}{C_{s}}=10}.

If we increase the overlap , $\delta=6h$, we see in Figure 5 on the right

that the optimized Schwarz method is now converging. The most slowly converging mode is shown on the left in Figure 5, and it also corresponds to a mode $|\sin(ky)|$ along the interface with $k=4\pi\approx 12>\frac{\omega}{C_{s}}=10$, so our Fourier analysis captures accurately the convergence behavior of the optimized Schwarz method.

4.2 Comparing Schwarz as solver and preconditioner

We next compare the performance of the Schwarz methods as solvers and preconditioners. We simulate the wave propagation through a computational domain given by the unit square $\Omega:=(0,1)^{2}$ with absorbing boundary conditions $\left({\cal T}^{({\bf n})}-\rm i\mathbf{\sigma}_{{\bf n}}\right){\bf u}=\mathbf{g}$ , where in the two-dimensional case considered here

[TABLE]

The source term $\mathbf{g}$ is chosen such that the exact solution is a plane wave $\mathbf{u}^{inc}$ consisting of both P- and S-waves, \mathbf{u}^{inc}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}:}={\bf d}\,e^{i\kappa_{p}{\bf x}\cdot{\bf d}}+{\bf d}^{\perp}\,e^{i\kappa_{s}{\bf x}\cdot{\bf d}},\,{\bf d}=\left(\cos\left(\frac{\pi}{3}\right),\cos\left(\frac{\pi}{3}\right)\right)^{T}. We choose the physical parameters $C_{p}=1$, $C_{s}=0.5$, $\rho=1$, $\lambda=\rho(C_{p}^{2}-2C_{s}^{2})$, $\mu=\rho C_{s}^{2}$, and $\omega=5$. We decompose the square domain $\Omega$ into $4\times 4$ equal subdomains $\Omega_{i}$ having each $40\times 40$ discretization points for a total number of 6400 degrees of freedom per subdomain. The convergence of the Schwarz algorithms as solvers and preconditioners for GMRES for different values of the overlap is shown in Figure 6.

As expected, the optimized Schwarz algorithm as solver converges, and the classical Schwarz algorithm diverges, for any size of the overlap. By increasing the overlap, as predicted by our two subdomain analyses in Theorem 6 and 2, the optimized Schwarz algorithm is getting better, whereas classical Schwarz is getting worse. With GMRES acceleration, overlap also helps the classical Schwarz algorithm, but it still takes substantially more iterations to converge than the optimized one.

4.3 Solving a circular transmission problem

We finally test our Schwarz methods for the Navier equations on a transmission problem formed by a circular inner part with radius 0.5 that has different material characteristics from the surrounding outer part, truncated with absorbing boundary conditions at the radius 1. The heterogeneous physical parameters are given in Table 1.

We use METIS to partition the unit disk $\{(x,y)|x^{2}+y^{2}\leq 1\}$ into $4$ subdomains as shown in Figure 7 on the left.

The solution of the transmission problem we compute is shown in Figure 7 on the right. We test the different Schwarz methods again both as solvers and as preconditioners for GMRES; the corresponding results are shown in Figure 8.

We see again that only the optimized Schwarz method with TTC converges when used as an iterative solver, the classical one diverges. This leads then naturally to a much better preconditioner for GMRES in the optimized Schwarz case for solving the transmission problem.

5 Conclusions

We presented a first study of the applicability of Schwarz methods for the solution of time-harmonic elastic waves modeled by the Navier equations. We showed by a detailed and technical analysis for two subdomains that the classical Schwarz method can not converge when applied to the Navier equations. We then introduced more physical transmission conditions and showed that optimal transmission conditions exist which make the algorithm converge in two steps. Since these optimal transmission conditions involve non-local operators, we also introduced a local, low-frequency approximation, and proved that the new, optimized Schwarz method is then convergent, provided the overlap is large enough. We then tested the Schwarz methods both for the two subdomain case, and also for many subdomains, including a heterogeneous transmission problem, and we observed numerically that the new, optimized Schwarz method can indeed be used as an iterative solver, while the classical one can not, since it is divergent. The new transmission conditions lead also to a much better Schwarz preconditioner for GMRES than the classical ones. Our analysis opens the path to further development, namely transmission conditions which do not only improve the low frequency behavior, but improve the convergence over the entire spectrum of the iteration operator, a topic which we are currently investigating.