M J Gander

A Tonnoir

ANALYSIS OF SCHWARZ METHODS FOR CONVECTED HELMHOLTZ LIKE EQUATIONS

Keywords: Complex advection, convected Helmholtz equation, Schwarz methods. MSC codes. 65M55, 65N55, 65F10

We present and analyze Schwarz domain decomposition methods for a general diffusion problem with complex advection. The complex advection term changes completely the nature of the solution and makes it more Helmholtz like. We analyze in detail the influence of the outer boundary conditions on the performance of the Schwarz algorithm, including PML conditions to emulate free space problems, and optimized transmission conditions, also for multiple subdomains. Our results show that the performance of Schwarz methods for such Helmholtz like problems is much better on free space configurations than in waveguides or closed cavities. Equations with complex advection appear in diverse applications, for example the convected Helmholtz equation, the Gross-Pitaevskii equation, Schrödinger equations, and also as important component in the wave-ray multigrid algorithm for Helmholtz problems. We show as an example the performance of our Schwarz methods for a potential flow around a schematic submarine.

1. Introduction. We are interested in solving numerically a partial differential equation (PDE) with a complex (!) advection term of the form (1.1) -div(A∇u) + ıa

• ∇u + µu = f in Ω, ı := √ -1,
where Ω is a subset of R 2 , A is a 2 × 2 positive definite matrix function, µ is a real function and a is a vector function in R 2 . We will assume that the source term f is compactly supported, and (1.1) must be equipped with appropriate boundary conditions that we will specify later. Equation (1.1) is very different from a classical advection diffusion equation with real advection term, and can have Helmholtz character even when µ has the good sign, i.e. µ ≥ 0. Equation (1.1) appears in various contexts:

• The convected Helmholtz equation: in this case, µ = -ω 2 with ω the pulsation of the wave, a = -2ωv with v the underlying flow (with convention e -ıωt for the time variable), and the solution u represents a pressure field. If the underlying flow is assumed to be incompressible, then we have

A = c 2 0 Id -v v T ,
with c 0 > 0 the sound speed, see e.g. [START_REF] Marchner | Stable Perfectly Matched Layers with Lorentz transformation for the convected Helmholtz equation[END_REF][START_REF] Barucq | HDG and HDG+ methods for harmonic wave problems with convection[END_REF][START_REF] Bécache | Perfectly matched layers for the convected Helmholtz equation[END_REF]. Note that to ensure that the matrix A is positive definite, the flow speed v must be small enough with respect to the sound speed c 0 (under mach 1).

• The Gross-Pitaevskii equation: equation (1.1) also appears as an intermediate problem for computing ground states of the Gross-Pitaevskii equation (which consists in solving a minimization problem), see [11, p.1107] or [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]. Solving equation (1.1) is an essential ingredient to compute the Sobolev gradient of the cost functional.

• The linearized Schrödinger equation: when looking for traveling wave solutions of the form ψ(t, x) = u(x-at) to the linearized Schrödinger equation ı∂ t ψ + 1 2 ∆ψ -V ψ = 0, see [START_REF] Bethuel | Travelling waves for the Gross-Pitaevskii equation I[END_REF] or [10, p.198], equation (1.1) appears with A = 1 2 Id and µ = V .

• The ray equation: equation (1.1) also appears as a fundamental ingredient in the wave-ray multigrid method for solving the Helmholtz equation [START_REF] Brandt | Wave-ray multigrid method for standing wave equations[END_REF][START_REF] Livshits | An algebraic multigrid wave-ray algorithm to solve eigenvalue problems for the Helmholtz operator[END_REF][START_REF] Livshits | Accuracy properties of the wave-ray multigrid algorithm for Helmholtz equations[END_REF][START_REF] Verburg | Multi-level Wave-Ray method for 2d Helmholtz equation[END_REF],

-∆ũ -ω 2 ũ = 0, when seeking the ray component of the form ũ(x) = e ık•x u, where k is a given direction in R 2 satisfying the dispersion relation k 2 2 = ω 2 . Then, u satisfies equation (1.1) with A = Id, a = -2k and µ = 0.

When equipped with classical Dirichlet, Neumann or Robin (impedance) boundary conditions (BCs), one can show that problem (1.1) is of Fredholm type, since the operator -div(A∇•) + • is coercive. We deduce then that the problem is of type coercive + compact, see [3, p.6], and therefore admits a unique solution, except for at most a countable set of parameters.

Depending on the situation of interest from the list above, we will consider Dirichlet, Neumann or Robin BCs. Furthermore, we will also consider the case where we have Perfectly Matched Layers (PMLs) surrounding the domain of interest, which is important for wave-like problems on unbounded domains. The derivation of the PML formulation is not straightforward for equations of the type (1.1), see for instance [START_REF] Bécache | Perfectly matched layers for the convected Helmholtz equation[END_REF][START_REF] Marchner | Stable Perfectly Matched Layers with Lorentz transformation for the convected Helmholtz equation[END_REF] for the convected Helmholtz equation, and we will briefly recall the PML construction hereafter.

Our goal is to analyze convergence properties of a Schwarz Domain Decomposition Method (DDM) with overlap using classical Fourier analysis, see [START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF][START_REF] Dolean | An introduction to domain decomposition methods: algorithms, theory, and parallel implementation[END_REF]. In particular, we wish to emphasize the impact of considering PML to truncate the computational domain. There is an important body of literature dedicated to the study of Schwarz methods for the Helmholtz equation, see [START_REF] Gander | Schwarz methods by domain truncation[END_REF][START_REF] Gong | Convergence of parallel overlapping domain decomposition methods for the Helmholtz equation[END_REF][START_REF] Graham | Domain decomposition with local impedance conditions for the Helmholtz equation with absorption[END_REF][START_REF] Gander | A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF][START_REF] Graham | Domain decomposition preconditioning for highfrequency Helmholtz problems with absorption[END_REF][START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF][START_REF] Gander | An optimized Schwarz method with twosided Robin transmission conditions for the Helmholtz equation[END_REF][START_REF] Gander | Optimized Schwarz methods without overlap for the Helmholtz equation[END_REF][START_REF] Després | Domain decomposition method and the Helmholtz problem[END_REF] and reference therein. However, only few results exist for the convected Helmholtz equation; an exception is the recent paper [START_REF] Lieu | A nonoverlapping Schwarz domain decomposition method with high-order finite elements for flow acoustics[END_REF], in which the authors study a nonoverlapping DDM for the convected Helmholtz equation. In fact, in the case of constant parameters A, a and µ, one can reformulate, as we explain hereafter, equation

(1.1) as a classical Helmholtz equation, using an appropriate change of variables. This

shows in particular that we will clearly face the same difficulties as in the Helmholtz case [START_REF] Ernst | Why it is difficult to solve Helmholtz problems with classical iterative methods[END_REF] for solving (1.1), but we can also benefit from the results known for the Helmholtz case.

The rest of our paper is organized as follows: First in section 2 we recall the link between the classical Helmholtz equation and equation (1.1), and explain how we can derive a stable PML formulation and first order Absorbing BCs (ABCs). Then, in section 3 we present a Fourier analysis of a Schwarz DDM considering vertical slicing and Robin transmission conditions. We study the impact of various outer PML truncations on the performance of the method, and explain how to properly take them into account in the implementation. Finally, in section 4 we give some concluding remarks.

Remark 1.1. Schwarz methods have been intensively studied for a formally similar equation, namely the advection-diffusion equation, see [START_REF] Achdou | A domain decomposition preconditioner for an advection-diffusion problem[END_REF][START_REF] Gerardo-Giorda | Optimized Schwarz methods for unsymmetric layered problems with strongly discontinuous and anisotropic coefficients[END_REF][START_REF] Dubois | Optimized Schwarz methods with Robin conditions for the advection-diffusion equation[END_REF], but the mathematical character of this equation with real advection term is very different from our equation (1.1). Also, the anisotropic aspect of diffusion was studied for Schwarz methods in [START_REF] Gander | Optimized Schwarz methods with general Ventcell transmission conditions for fully anisotropic diffusion with Discrete Duality Finite Volume discretizations[END_REF][START_REF] Gander | Discrete optimization of Robin transmission conditions for anisotropic diffusion with Discrete Duality Finite Volume methods[END_REF], but again without the fundamentally character changing term of the complex advection in (1.1).

2. Reformulation as a Helmholtz equation and related results. In this section, as well as for the analysis in the next section, we will suppose that A, a and µ are constant parameters. Note that for the construction of ABCs and for the PML formulation, a generalization to locally perturbed parameters is possible.

2.1. Link with the Helmholtz equation. For the convected Helmholtz equation, in the case of constant parameters, it is well-known that there exist coordinate transformations that map the convected Helmholtz equation into the Helmholtz equation, see [START_REF] Marchner | Stable Perfectly Matched Layers with Lorentz transformation for the convected Helmholtz equation[END_REF][START_REF] Hu | On the use of a Prandtl-Glauert-Lorentz transformation for acoustic scattering by rigid bodies with a uniform flow[END_REF]. Let us explain a similar idea for (1.1), namely to consider u(x) = v(x)e ık•x with k to be suitably chosen later. Then, we have

∇u = (∇v)e ık•x + ve ık•x ık, and
div(A∇u) = div(A∇v)e ık•x + 2ıA∇v • ke ık•x -v k 2 A e ık•x ,
where we denote by k 2 A := Ak • k. Introducing these results into (1.1), we get that v satisfies

(2.1) -div(A∇v) + ∇v • (ia -2ıAk) + (µ + k 2 A -a • k)v = e -ık•x f. If we choose now k := 1 2 A -1 a, equation (2.1) simplifies to (2.2) -div(A∇v) + (µ - a 2 A -1 4 )v = e -ıA -1 a 2 •x f. Remark 2.1. Note that (2.2) is a classical Helmholtz problem if µ - a 2
A -1

4

< 0, so even for µ ≥ 0, (1.1) has a Helmholtz character if

a 2 A -1 4
is large enough: the complex convection term is the reason for this Helmholtz character of equation (1.1). Now, from (2.2), using an appropriate linear coordinate transformation of the form

(2.3) x = Sx,
we can rewrite the operator div(A∇•) as a Laplace operator, because

(2.4) ∇• = S T ∇ • =⇒ div(A∇•) = div (SAS T ∇ •),
and since A is symmetric positive definite, we can use the Cholesky decomposition A = GG T and thus take S = G -1 to simplify (2.4). Equation (2.2) then simplifies to

(2.5) -∆ v + (µ - 1 4 a 2 A -1 )v = f , where v (x ) = v(x) and f (x ) = e -ıA -1 a 2 •S -1 x f (S -1 x ).
Remark 2.2. In our case, the change of variables (2.3) preserves the vertical boundaries, because the matrix S = G -1 is lower triangular, so that

{x = α} ⇐⇒ {x = s 11 α},
where s ij = [S] i,j . In contrast, horizontal boundaries are deformed into oblique boundaries (if a 12 = 0), in the same spirit as in [START_REF] Tonnoir | Dirichlet-to-Neumann operator for diffraction problems in stratified anisotropic acoustic waveguides[END_REF].

3
This manuscript is for review purposes only. for the convected Helmholtz equation, see [START_REF] Barucq | Prandtl-Glauert-Lorentz based absorbing boundary conditions for the convected Helmholtz equation[END_REF], or PMLs [START_REF] Bécache | Perfectly matched layers for the convected Helmholtz equation[END_REF][START_REF] Marchner | Stable Perfectly Matched Layers with Lorentz transformation for the convected Helmholtz equation[END_REF], albeit using a different change of variables (the choice is not unique). Suppose that 1 µ -1 4 a 2 A -1 < 0, then we can easily deduce the equivalent of the classical ABC for the Helmholtz equation (taking the convention e -ıωt for the time variable):

∇ v • n -ı ωv = 0 ⇐⇒ S -T ∇v • n -ı ωv = 0, ⇐⇒ A∇v • (G -T n ) =n/ G T n -ı ωv = 0, ⇐⇒ A∇u • n -ı 1 2 a • nu -ı ω n A u = 0, (2.6)
where ω := -µ + [START_REF] Barucq | Prandtl-Glauert-Lorentz based absorbing boundary conditions for the convected Helmholtz equation[END_REF][START_REF] Marchner | Construction and numerical assessment of local absorbing boundary conditions for heterogeneous time-harmonic acoustic problems[END_REF] for more details and higher order ABCs.

As a numerical illustration, we consider the convected Helmholtz problem

(2.7) -div (A∇u) -2ıωv • ∇u -ω 2 u = δ in Ω = (0, 4) 2 , A∇u • n + ıωv • nu -ı ω n A u = 0 on ∂Ω,
where δ is the Dirac source term, and A = Id -vv T . We show in Figure 1 the solution u we obtain using the ABC (2.6) for the model problem (2.7) with problem parameters

(2.8) v := Ma cos(θ) sin(θ) , θ = π 4 , Ma = 1 2
, ω = 10.

Derivation of a Cartesian PML formulation.

In this subsection, we will also assume that A is a diagonal matrix so that the coordinate transformation 1 In what follows, we will always assume that µ -1 4 a 2 A -1 < 0, since otherwise, the problem is coercive and has lost its difficult Helmholtz character.

(2.3) preserves both horizontal and vertical boundaries, since the matrix S is diagonal, which simplifies the construction of a Cartesian PML. We refer to [START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains[END_REF] for the construction of a PML on a polygonal domain. We emphasize that the assumptions that the parameters are constant and A is diagonal are necessary only in the PML region.

Under these hypotheses, it is well-known that the PML formulation for the Helmholtz equation in (x , y ) coordinates reads (2.9)

-div (D P M L ∇ v ) + s x s y (µ - 1 4 a 2 A -1 )v = f ,
where

D P M L = s y /s x 0 0 s x /s y .
The complex valued functions s x and s y are defined by (2.10)

s x (x ) := 1 if x ∈ (a + , b -), 1 + ıσ x otherwise, and (2.11) s y (y ) := 1 if y ∈ (c + , d -), 1 + ıσ y otherwise,
where σ x > 0 and σ y > 0 are the strength of the PML in each direction, and > 0 is the depth of the PML. In (x , y ) coordinates, the computational domain would then be the square (a , b ) × (c , d ). Returning to the (x, y) coordinates, we get

-div (D P M L ∇ v ) + s x s y (µ - 1 4 a 2 A -1 )v = f , ⇐⇒ -div (A P M L ∇v) + s x s y (µ - 1 4 a 2 A -1 )v = f , (2.12)
where

A PML = a 11 s y /s x 0 0 a 22 s x /s y ,
and s x (x) = s x (x ) and s y (y) = s y (y ). Recalling that v(x) = u e -i 1 2 A -1 a•x , we get

A PML ∇v = A PML ∇u -ıu 1 2 A PML A -1 a e -ı 1 2 A -1 a•x ,
and, observing that A PML depends on (x, y),

div(A PML ∇v) = div(A PML ∇u) -ı 1 2 a PML • ∇u -ı 1 2 div (ua PML ) -u 1 4 a 2 A -1 PML e -i 1 2 A -1 a•x ,
where Thus, inserting this expression into (2.12), we obtain the PML formulation for (1.1), namely

A -1 PML := A -1 A PML A -1 , a PML := A PML A -1 a.
-div(A PML ∇u) + ı 1 2 a PML • ∇u+ı 1 2 div(u a PML ) +s x s y (µ - 1 4 a 2 A -1 + 1 4s x s y a 2 A -1 PML )u = f.
(2.13)

Note that the PDE remains unchanged from the initial one in (1.1) in the physical region.

As a numerical illustration, let us consider once again problem (2.7), but this time with the problem parameters

(2.14) v := M cos(θ) sin(θ) , θ := 0, M := 4 5
, ω := 20.

In Figure 2, we show the real part of the solution computed using a naive classical PML (left), which is known to have instabilites in some configurations, as we can clearly see here, and the PML formulation (2.13) (right), which works perfectly.

3. Fourier analysis of a classical Schwarz algorithm. We now present and analyze Schwarz domain decomposition methods for (1.1) in a specific geometry: find

u ∈ H 1 (Ω) satisfying (3.1) -div(A∇u) + ıa • ∇u + µu = f in Ω = (a, b) × (c, d),
where A is a diagonal matrix, a < b and c < d with {a, b, c, d} ∈ R. For the boundary conditions, we will consider four configurations:

• Dirichlet-Dirichlet: we impose homogeneous Dirichlet boundary conditions on both vertical and horizontal boundaries,

• Dirichlet-PML: we impose also homogeneous Dirichlet boundary conditions on the left and right but a PML on the bottom and top boundaries (which terminates with a homogeneous Dirichlet boundary condition),

• PML-Dirichlet: the same idea but with PML on the vertical boundaries,

• PML-PML: imposing PML on all sides of the domain.

The first case models a bounded domain, the second and third cases a waveguide with different orientation, and the last case a free space problem. We decompose the domain Ω first into two overlapping subdomains This manuscript is for review purposes only. interface of Ω 1 (within Ω 2 ) and Γ 2,1 := {x = α}×(0, 1) the interface of Ω 2 (within Ω 1 ).

a b c d Γ 1,2 α Γ 2,1 β a b c d Γ 1,2 α Γ 2,1 β a b c d Γ 1,2 α Γ 2,1 β a b c d Γ 1,2 α Γ 2,1 β
Then, a general iterative Schwarz algorithm computes for iteration index n = 1, 2, . . .

the subdomain solutions (3.2) -div(A∇u n 1 ) + ıa • ∇u n 1 + µu n 1 = f 1 in Ω 1 , u n 1 = 0 on ∂Ω 1 \ Γ 1,2 , (a 11 ∂ x + p 1,2 -ı a1 2 )u n 1 = (a 11 ∂ x + p 1,2 -ı a1 2 )u n-1 2 on Γ 1,2 , -div(A∇u n 2 ) + ıa • ∇u n 2 + µu n 2 = f 2 in Ω 2 , u n 2 = 0 on ∂Ω 2 \ Γ 2,1 , (-a 11 ∂ x + p 2,1 + ı a1 2 )u n 2 = (-a 11 ∂ x + p 2,1 + ı a1 2 )u n 1 on Γ 2,1 ,
where f i is the restriction of f to Ω i , i ∈ {1, 2}, and p 1,2 , p 2,1 are complex constants.

We emphasize that when using PML, the algorithm should be written with complex stretched coordinates, or equivalently with the PML formulation as described in subsection 2.3. Also, note that we consider here a particular Robin type transmission condition to get a condition similar to the ABC (2.6) at the interfaces. Moreover, in the PML formulation it is interesting to note that the boundary term coming from the integration by parts of div(ua P M L ) is canceled by this choice of transmission condition, see Remark 3.3 for more details.

To study the convergence of the Schwarz algorithm (3.2) as n goes to infinity, we consider the error u -u n i | Ωi , i ∈ {1, 2}, which amounts to consider the algorithm (3.2) with zero source terms. Using the equivalence with the Helmholtz equation, the iterative algorithm (3.2) for the error becomes

(3.3) -∆ (v ) n 1 -ω 2 (v ) n 1 = 0 in Ω 1 , (v ) n 1 = 0 on ∂Ω 1 ∩ Ω , ∂ x + p 1,2 (v ) n 1 = ∂ x + p 1,2 (v ) n-1 2 on Γ 1,2 , -∆ (v ) n 2 -ω 2 (v ) n 2 = 0 in Ω 2 , (v ) n 2 = 0 on ∂Ω 2 ∩ Ω , -∂ x + p 2,1 (v ) n 2 = -∂ x + p 2,1 (v ) n 1 on Γ 2,1 , where Γ 1,2 := {x = β } × (c , d ) and Γ 2,1 := {x = α } × (c , d ), and (3.4) p 2,1 = p 2,1 g 11 and p 1,2 = p 1,2 g 11 ,
and we recall that (g ij ) ij are the coefficients of the lower triangular matrix G from the Cholesky decomposition A = GG T . As a consequence, to study the convergence of the Schwarz algorithm (3.2), we will study the convergence of the algorithm rewritten for the Helmholtz equation (3.3). A similar idea of using an equivalent algorithm to remove the anisotropy and advection term can be found in [START_REF] Gander | Optimized Schwarz methods for a diffusion problem with discontinuous coefficient[END_REF]. Remark 3.3. For the implementation, in the PML context, the transmission conditions on Γ 1,2 and Γ 2,1 should be equivalently rewritten as

s y s x a 11 ∂ x u n 1 + s y s x (p 1,2 -ı a 1 2 )u n 1 = s y s x a 11 ∂ x u n-1 2 + s y s x (p 1,2 -ı a 1 2 )u n-1 2 on Γ 1,2 ,
and

- s y s x a 11 ∂ x u n 2 + s y s x (p 2,1 + ı a 1 2 )u n 2 = - s y s x a 11 ∂ x u n 1 + s y s x (p 2,1 + ı a 1 2 )u n 1 on Γ 2,1 ,
to get natural variational conditions. This is different from implementing the Després like transmission conditions

s y s x a 11 ∂ x u n 1 + (p 1,2 -ı a 1 2 )u n 1 = s y s x a 11 ∂ x u n-1 2 + (p 1,2 -ı a 1 2 )u n-1 2 on Γ 1,2 ,
and

- s y s x a 11 ∂ x u n 2 + (p 2,1 + ı a 1 2 )u n 2 = - s y s x a 11 ∂ x u n 1 + (p 2,1 + ı a 1 2 )u n 1 on Γ 2,1 ,
which can lead to a divergent algorithm when algorithm (3.2) is convergent! 3.1. Computation of the convergence factor in the two subdomain case.

We show the computations for the PML-PML case, the other cases can be deduced by simply taking σ x = 0 or σ y = 0, see equation (2.10) for the definition of σ x and σ y .

We introduce the complex stretched coordinates in the modified coordinate system (x , y ),

(3.5) x (x ) =    x + ıσ(x -(a + )) if x ∈ (a , a + ), x if x ∈ (a + , b -), x + ıσ(x -(b -)) if (b -, b ), and (3.6) y (y ) =    y + ıσ y (y -(c + )) if y ∈ (c , c + ), y if y ∈ (c + , d -), y + ıσ y (y -(d -)) if y ∈ (d -, d ).
Due to the rectangular geometry of the domain Ω , and since A is assumed to be diagonal, we can use separation of variables to analytically obtain the errors in the Schwarz algorithm (3.3),

(3.7) (v ) n i = k∈N * ψ k (y ) A n i (k)e ıλ(ξ k ) x (x ) + B n i (k)e -ıλ(ξ k ) x (x ) , i ∈ {1, 2}, where λ(ξ k ) = ω 2 -ξ 2 k .
The functions ψ k and the complex numbers ξ k are the eigenfunctions and eigenvalues of the eigenvalue problem

(3.8) -∂ 2 y y ψ k = ξ 2 k ψ k for y ∈ (c , d ), ψ k = 0 on y ∈ {c , d }, 8
This manuscript is for review purposes only. c ,d )). This does not hold any more when considering PML (σ y > 0). In fact, although one can show that the family is a complete basis [START_REF] Knockaert | On the completeness of eigenmodes in a parallel plate waveguide with a perfectly matched layer termination[END_REF], it is neither an orthonormal basis nor a Riesz basis. A consequence of this result is that the decomposition (3.7) is still justified, but cannot be computed in practice given an arbitrary Robin data on Γ 1,2 or Γ 2,1 .

In the expressions (3.7), the amplitudes A n i (k) and B n i (k) should be chosen to satisfy the vertical BCs, namely

• on Γ 0 = {x = a } × (c , d ) and Γ 1,2 for i = 1,
• and on

Γ 2 = {x = b } × (c , d ) and Γ 2,1 for i = 2,
the horizontal BCs on (a , b ) × {c , d } being already satisfied. To ensure these BCs, we must impose

B n 1 (k) = -A n 1 (k)e 2ıλ(ξ k )(a -iσx ) and B n 2 (k) = -A n 2 (k)e 2iλ(ξ k )(b +iσx ) .
Now, the BC on Γ 1,2 (on x = β ) implies

A n 1 (k) ıλ(ξ k ) e ıλ(ξ k )β + e 2ıλ(ξ k )(a -ıσx ) e -ıλ(ξ k )β +p 1,2 e ıλ(ξ k )β -e 2ıλ(ξ k )(a -ıσx ) e -ıλ(ξ k )β = A n-1 2 (k) ıλ(ξ k ) e ıλ(ξ k )β +e 2ıλ(ξ k )(b +ıσx ) e -ıλ(ξ k )β +p 1,2 e ıλ(ξ k )β -e 2ıλ(ξ k )(b +ıσx ) e -ıλ(ξ k )β , so that A n 1 (k) = ρ 1 (k)A n-1 2 (k),
with the first convergence factor component

(3.11) ρ 1 (k) = ıλ(ξ k ) + p 1,2 + e 2ıλ(ξ k )(b +ıσx -2β ) ıλ(ξ k ) -p 1,2 ıλ(ξ k ) + p 1,2 + e 2ıλ(ξ k )(a -iσx -2β ) ıλ(ξ k ) -p 1,2 .
In the same way, we get using the BC on Γ 2,1 (on x = α ) that

A n 2 (k) -ıλ(ξ) e iλ(ξ k )α + e 2ıλ(ξ k )(b +ıσx ) e -ıλ(ξ k )α +p 2,1 e ıλ(ξ k )α -e 2ıλ(ξ k )(b +ıσx ) e -iλ(ξ k )α = A n 1 (k) -ıλ(ξ k ) e ıλ(ξ k )α +e 2ıλ(ξ k )(a -ıσx ) e -ıλ(ξ k )α +p 2,1 e ıλ(ξ k )α -e 2ıλ(ξ k )(a -ıσx ) e -ıλ(ξ k )α , 9
This manuscript is for review purposes only. so that

A n 2 (k) = ρ 2 (k)A n 1 (k),
with the second convergence factor component

(3.12) ρ 2 (k) = -ıλ(ξ k ) + p 2,1 -e 2ıλ(ξ k )(a -ıσx -α ) ıλ(ξ k ) + p 2,1 -ıλ(ξ k ) + p 2,1 -e 2ıλ(ξ k )(b +ıσx -α ) ıλ(ξ k ) + p 2,1 .
Thus, the convergence factor of the Schwarz method is ρ = 0.1. We see that the convergence factor is highly influenced by using a PML on the outer boundary. In particular, the more we open the domain by adding PMLs, the better the convergence factor becomes. This can be understood for wave like problems in the sense that when the domain is open, error components can leave the domain to infinity, or equivalently they are damped by the PML which emulates the unbounded domain. Other boundary conditions reflect these error components and inject them back into the iteration, leading to worse convergence, or even divergence.

(k) = ρ 1 (k)ρ 2 (k).
A second remark we can make for waveguide problems, corresponding to the D-PML and PML-D cases, is that cutting the waveguide in the infinite direction or in the transverse direction is very different. Indeed, in the PML-D case, the convergence factor is good for small k, whereas in the the D-PML case, the convergence factor is better for larger k, A last remark is the fact that computing the convergence factor in a vertical waveguide Ω = (a , b ) × R using a Fourier transform in the y -direction would lead exactly to the same convergence factor as in the D-D case, since the only change is the continuous summation with eigenfunctions e iξy which replaces the discrete summation, but the computed solution and performance of the Schwarz method is very different in an open wave guide or a closed cavity. In contrast, using a horizontal PML as in the D-PML case leads to a very different convergence factor, whereas the computed Schwarz iterates in the physical domain correspond to the solution in the unbounded domain! This shows that the two-subdomain analysis is very different if we consider the PML or not. Remark 3.5. Note that if the PML parameters σ x and σ y are too large, then the convergence factor deteriorates. In particular, if σ x = σ y = σ → +∞, we do not recover the convergence factor one would get in the full space R 2 , as in [START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF]. However,

11
This manuscript is for review purposes only.

we recover this convergence factor if the length of the PML tends to +∞.

3.2.

Generalization to more subdomains. The Fourier analysis above can be generalized to more subdomains if we still consider vertical slicing of the domain to allow us to use separation of variables. Let us consider N s subdomains Ω

i = (α i , β i ) × (c, d), i ∈ {1, • • • , N S }, where a = α 1 < α 2 < β 1 < α 2 < • • • < α Ns < β Ns-1 < β Ns = b.
In that case, for simplicity we will consider the parallel version of algorithm (3.2), which gives for its equivalent Helmholtz formulation the error equations

(3.15) -∆ (v ) n i -ω 2 (v ) n i = 0 in Ω i , (v ) n i = 0 on ∂Ω i ∩ Ω , (∂ x + p i,i+1 )(v ) n i = (∂ x + p i,i+1 )(v ) n-1 i+1 on Γ i,i+1 , (-∂ x + p i,i-1 )(v ) n i = (-∂ x + p i,i-1 )(v ) n-1 i-1 on Γ i,i-1 , where Γ i,i-1 := {x = α i } × (c , d ) and Γ i,i+1 := {x = β i } × (c , d ).
Then, with separation of variables, we still have as in (3.7), for all i ∈ {1,

• • • , N s } (3.16) (v ) n i = k∈N * ψ k (y ) A n i (k)e ıλ(ξ k ) x (x ) + B n i (k)e -ıλ(ξ k ) x (x ) .
Here again, the horizontal BCs are satisfied by definition of ψ k . For each mode, the

BC on {x = a } × (c , d ) imposes that (3.17) B n 1 (k) = -A n 1 (k)e 2ıλ(ξ k )(a -iσx ) , whereas the BC on {x = b } × (c , d ) imposes that (3.18) B Ns,n (k) = -A Ns,n (k)e 2ıλ(ξ k )(b +iσx ) .
The transmission conditions on Γ i,i+1 give

A n i (k) ıλ(ξ k ) + p i,i+1 e ıλ(ξ k ) x (βi) + B n i (k) -ıλ(ξ k ) + p i,i+1 e -ıλ(ξ k ) x (βi) = A n-1 i+1 (k) ıλ(ξ k ) + p i,i+1 e ıλ(ξ k ) x (βi) + B n-1 i+1 (k) -ıλ(ξ k ) + p i,i+1 e -ıλ(ξ k ) x (βi) ,
and similarly the transmission conditions on Γ i,i-1 give

A n i (k) -ıλ(ξ k ) + p i,i-1 e ıλ(ξ k ) x (αi) + B n i (k) ıλ(ξ k ) + p i,i-1 e -ıλ(ξ k ) x (αi) = A n-1 i-1 (k) -ıλ(ξ k ) + p i,i-1 e ıλ(ξ k ) x (αi) + B n-1 i-1 (k) ıλ(ξ k ) + p i,i-1 e -ıλ(ξ k ) x (αi) .
Combining these relations, we get for the mode k the iteration relation

c n (k) = I(k) c n-1 (k), where I(k) = D -1 (k)K(k), 12 where c n (k) := A n 1 (k) B n 1 (k) • • • A Ns,n (k) B Ns,n (k) T . The matrix D is block diagonal D =          D 1 0 • • • • • • 0 0 D 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • • • • 0 D Ns          , where the matrices D i are 2 × 2 matrices s.t. for all i ∈ {2, • • • , N s -1} D i = (-ıλ(ξ k ) + p i,i-1 )e 2ıλ(ξ k ) x (αi) ıλ(ξ k ) + p i,i-1 (ıλ(ξ k ) + p i,i+1 )e 2ıλ(ξ k ) x (βi) -ıλ(ξ k ) + p i,i+1 ,
and, to take into account the Dirichlet BCs (3.17) and (3.18), we have

D 1 = e 2ıλ(ξ k ) x (α1) 1 (ıλ(ξ k ) + p 1,2 )e 2ıλ(ξ k ) x (β1) -ıλ(ξ k ) + p 1,2
and

D Ns = (-ıλ(ξ k ) + p Ns,Ns-1 )e 2ıλ(ξ k ) x (α Ns ) ıλ(ξ k ) + p Ns,Ns-1 e 2ıλ(ξ k ) x (β Ns ) 1 .
Similarly, the matrix K is given by

K =          0 K 1,2 0 • • • 0 K 2,1 0 K 2,3 . . . 0 . . . . . . . . . . . . . . . . . . . . . K Ns-1,Ns 0 • • • • • • K Ns,Ns-1 0          , where K i,i-1 = (-ıλ(ξ k ) + p i,i-1 )e 2ıλ(ξ k ) x (αi) ıλ(ξ k ) + p i,i-1 0 0 and K i,i+1 = 0 0 (ıλ(ξ k ) + p i,i+1 )e 2ıλ(ξ k ) x (βi) -ıλ(ξ k ) + p i,i+1
.

Thus the convergence factor is the spectral radius of the matrix I(k),

(3.19) ρ(k) = ρ(I(k)).
Remark 3.6. In the two subdomain case, N s = 2, if we eliminate B n 1 (k) and B n 2 (k) using the outer Dirichlet BCs, the iteration matrix becomes

0 ρ 1 (k) ρ 2 (k) 0 13
This manuscript is for review purposes only. where ρ 1 (k) and ρ 2 (k) are defined in (3.11) and (3.12). In particular, the convergence factor is in that case the square root of the convergence factor defined in the previous section. Note that this is simply linked to the fact that in this section, we have considered the parallel version of the Schwarz algorithm, whereas before we studied the alternating version for two subdomains.

Remark 3.7. Let us also note that the matrix D(k) is not invertible if (and only if) p i,i-1 = -p i,i+1 = ±ıλ(ξ k ), which corresponds to the case where the subproblem in Ω i is not well-posed: the mode k is a non zero solution of the homogeneous problem. Moreover, if p i,i-1 = p i,i+1 = -ıλ(ξ k ), then D(k) is diagonal and one can show that ρ(k) → 0 as the length of the PML tends to +∞.

As a numerical illustration, let us consider again example (3.13) of the previous section with the same parameters, except that this time the domain is split into 5 subdomains. The subdomains are defined by (3.20)

α 1 = 0, α i = 0.2(i -1) -0.05 for i ∈ {2, • • • , 5}, β 5 = 1, β i = 0.2i + 0.05 for i ∈ {1, • • • , 4}.
For the parameters p i,i+1 and p i,i-1 , we chose them s.t. p i,i+1 = p i,i-1 = -ı ω, as before. In Figure 6, we show the convergence factor ρ(k) for the D-D, D-PML, PML-D and PML-PML configurations. As one could expect, the convergence factor is less good than for the two subdomain case, but the remarks for the two subdomain case still hold. In particular, the more we open up the domain with PML outer boundary conditions, the better the convergence becomes.

3.3. Optimized transmission conditions. Now that we have obtained the convergence factor, we can look for optimized parameters p i,i+1 and p i,i-1 , and deduce In particular, the optimized parameters will be different depending on the outer BCs.

Also, in contrast to the usual convergence factor in free space, see [START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF] for instance, here the convergence factor can be smaller than 1 at the cut-off frequency ξ k = w. This is due to the fact that we consider a bounded domain for the analysis. As a consequence, the optimization can be done for all k, as in [START_REF] Bouajaji | Optimized Schwarz methods for the time-harmonic maxwell equations with damping[END_REF], where the equation contained damping.

Finally, note that solving analytically this min-max problem is difficult, so we use a simple optimization process and the function fmin of scipy.optimize. For first results on a many subdomain optimization for a diffusive problem, see [START_REF] Dolean | Closed form optimized transmission conditions for complex diffusion with many subdomains[END_REF]. As a first numerical example, we consider the case of the Helmholtz equation, so

that u = v , (3.22) -∆u -ω 2 u = 0 in Ω = (0, 1) 2 .
We split the domain into 5 subdomains defined as before in equation (3.20). Taking ω = 50, we show in Figure 7 the convergence factor with p i,i+1 = p i,i-1 = -ıω and with optimized parameters in the different cases (for the PML cases, we take = 0.02 and σ x = σ y = 10). In particular for the PML-D case, one can see that optimized parameters allow us to get a convergent algorithm.

In Figure 8, we show the error evolution versus the iterations using either the Schwarz algorithm as iterative solver, or as preconditioner for GMRES. As we can see, in each case the optimized parameters improve the convergence for the Schwarz algorithm. Yet, this is no more true for GMRES. This can be explained since the optimization problem (3.21) optimizes the convergence of the iterative Schwarz algorithm. Therefore, when using it as a preconditioner for GMRES, a priori, we are not ensured that the optimized parameters are optimized parameters for GMRES. Remark 3.8. Let us emphasize that for the mesh discretization, one must consider a sufficiently fine mesh to get accurate results that match the theoretical convergence properties. In particular, if the mesh in the PML is too coarse, then the Schwarz algorithm can be divergent even if the continuous convergence factor is less than one.

As a second, more realistic example, let us consider the case of the convected Helmholtz equation

(3.23) -div (A∇u) -2ıωv • ∇u -ω 2 u = δ in Ω = (0, 1) 2 \ O, A∇u • n = 0 on ∂Ω,
where the obstacle O has the rough shape of a submarine, see Figure 9. We consider a potential flow v = ∇ϕ coming from the left, which we compute solving the Laplace problem

-∆ϕ = 0 in Ω, ∇ϕ • n = 0 on ∂O ∪ (0, 1) × {0, 1}, ∇ϕ • n = -1 on {0} × (0, 1), ∇ϕ • n = 1 on {1} × (0, 1),
with the same mesh. Note that to get a well-posed problem, we simply impose a value of ϕ inside Ω. We compute the gradient of ϕ inside each cell of the mesh to get v. Then, the velocity is normalized

v = Ma v v ∞ , where v ∞ = sup (x,y)∈Ω v(x, y) 2 ,
where Ma is the mach number, i.e. the ratio between the velocity of the fluid and the sound speed in the medium. Thus, the matrix A in (3.23) is given by

A = Id -v v T .
For this example, we took2 Ma = 0. which corresponds to an overlap of size 0.03. The mesh we use is unstructured, so that the interfaces between the subdomains are not perfectly straight any more. Moreover, since a is no more constant in the physical domain, a natural generalization of the transmission conditions, similar to the ABC (2.6), is

A P M L ∇u • n -ı 1 2 a P M L • n + p i,i±1 u.
In particular, we compare in this example the following choices of the parameters p i,i±1 :

• First, a classical ABC condition as in (2.6), (ABC) p i,i±1 = -ı ω n A where ω = -µ + 1 4 a 2 A -1 .

Note that w and n A are variable. where we recall that s y and s x are the PML parameters.

• Third, a condition that takes the PML into account, (OPT) p i,i±1 = s y s x q i,i±1 , where q i,i±1 are (numerical) solutions of the min-max problem (3.21) considering the medium with no obstacle and with a = -2ω( v ext , 0) T constant.

In Figure 10, we show the evolution of the residual considering these three transmission conditions. The best results are obtained with the condition ABC PML and OPT.

We see that clearly, taking into account the PML coefficient in the parameter is very important, as already mentioned in Remark 3.3.

4. Concluding remarks. We studied Schwarz domain decomposition methods for a general diffusion problem with complex advection, which appears in several important applications. The complex advection term changes fundamentally the nature of the diffusion problem and makes it Helmholtz like. We have shown that for such problems the outer boundary conditions imposed on the global domain have a strong influence on the convergence of the Schwarz method, and on how one should choose optimized parameters. Not taking into account the PML coefficients in the transmission conditions deteriorates the convergence of the Schwarz algorithm, both when used as iterative solver and as preconditioner for GMRES. Our analysis covers both two subdomain and many subdomain situations for decompositions into strips, and allowed us to formulate the min-max problem one has to solve to compute optimized parameters, which turns out to be difficult to treat theoretically. Furthermore, computing optimized parameters for GMRES is currently out of reach, for a special case in a splitting method, see [START_REF] Benzi | Optimization of the Hermitian and skew-Hermitian splitting iteration for saddle-point problems[END_REF].

Fig. 1 .

 1 Fig. 1. Example of using the ABC (2.6) with a point source in the center of the domain.

2. 2 .

 2 Derivation of a simple Absorbing Boundary Condition. The reformulation (2.5) as a Helmholtz equation is used in the literature for constructing ABCs

Fig. 2 .

 2 Fig. 2. Solution obtained using (on the left) the classical PML and (on the right) the modified PML.

  Ω 1 := (a, β) × (c, d) and Ω 2 := (α, b) × (c, d) with α ≤ β, see Figure 3. We denote by Γ 1,2 := {x = β} × (0, 1) the 6

Fig. 3 .

 3 Fig. 3. Domain decomposition for the model problem, from left to right: D-D, D-PML, PML-D and PML-PML. In blue we show the domain Ω 1 , in orange the domain Ω 2 and in gray the PML region. The overlapping area is delimited by the boundaries Γ 1,2 and Γ 2,1 .

Remark 3 . 1 .

 31 The choice of a diagonal matrix A ensures that in the reformulation (3.3) the domain Ω is still a square, which is important for the analytical solution we use below.Remark 3.2. From the reformulation (3.3) of the Schwarz algorithm (3.2), we can obtain optimized transmission conditions, using the optimized parameters p 1,2 and p 2,1 from[START_REF] Gander | Optimized Schwarz methods with overlap for the Helmholtz equation[END_REF] and relation(3.4).

  and we have, up to normalization (3.9) ψ k ∝ sin(ξ k (y -c )) and ξ k = kπ d -c if σ y = 0, and (3.10) ψ k ∝ sin(ξ k ( y (y ) -y (0))) and ξ k = kπ d -c + 2i σ y if σ y > 0. Remark 3.4. If there is no horizontal PML (σ y = 0), then the family (ψ k ) k is an orthonormal basis of L 2 ((

Fig. 4 .

 4 Fig. 4. Convergence factor ρ DD (on the left) with ρ DD (1) 1.11 (in blue), ρ DD (2) 0.463 (in orange) and ρ DD (3) 0.202 (in green). Error evolution versus iterations (on the right) initializing the error system with mode 1, 2 and 3 (curves in blue, orange and green respectively).

  Let us emphasize once again that the convergence factor depends on the case considered (D-D, D-PML, PML-D or PML-PML) through σ x but also through the eigenvalue problem (3.8). As a numerical illustration, let us consider the wave-ray equation, (3.13) -∆u + ıa • ∇u = 0 in Ω = (0, 1) 2 , u = 0 on ∂Ω, with a := (10, 0). To solve this problem, we implemented algorithm (3.2) (not its equivalent Helmholtz formulation (3.3)) with the parameters (3.14) p 12 = -ı ωg 11 and p 21 = -ı ωg 11 , which corresponds to p 12 = p 21 = -ı ω (and to the ABC (2.6)). We took α = 0.45 and β = 0.55 so that the overlap is of size 0.1. In the D-D case, the convergence factor is shown in Figure 4 (where the variable k is "continuified" in the abscissa). Since we have ρ(-k) = ρ(k), we show the convergence factor only for k ≥ 0. The vertical dotted lines correspond to integer values of k. We see that the algorithm is not convergent, since for the first (non-zero) mode, we have ρ DD (1) 1.115. In the same Figure, we also show the error evolution versus the iterations initializing the error equations with mode k = 1, 2 and 3. We see that for modes 2 and 3, the algorithm is convergent (as expected from the convergence factor) up to a point where the round-off error makes the first mode appear, like in power iterations. Computing the slopes of the three lines on the right gives the convergence factors ρ DD (1) 1.115, ρ DD (2) 0.463 and ρ DD (3) 0.202, matching well the theoretical prediction on the left. For the same example, we show the convergence factors ρ DD , ρ DP M L , ρ P M LD and ρ P M LP M L in Figure 5. For the PML parameters, we took σ x = σ y = 10 and 10This manuscript is for review purposes only.

Fig. 5 .

 5 Fig. 5. Convergence factor ρ DD (top left), ρ DP M L (top right), ρ P M LD (bottom left) and ρ P M LP M L (bottom right).

Fig. 6 .

 6 Fig. 6. Convergence factor ρ DD (top left), ρ DP M L (top right), ρ P M LD (bottom left) and ρ P M LP M L (bottom right) in the case of 5 subdomains.

Fig. 7 .

 7 Fig. 7. Convergence factor ρ DD (top left), ρ DP M L (top right), ρ P M LD (bottom left) and ρ P M LP M L (bottom right) in the case of 5 subdomains for the Helmholtz equation. Comparison between the classical transmission conditions -ıω (in blue) and optimized parameters (in red).

  p i,i+1 and p i+1,i from equation (3.4). More precisely, we have to solve the classical min-max problem (3.21) min (p 1,2 ,p 2,1 ,p 2,3 ,••• ,p Ns ,Ns -1 )∈C 2Ns max k∈N + |ρ(k)|.

Fig. 8 .

 8 Fig. 8. Relative residual versus iterations in the D-D case (top left), D-PML case (top right), PML-D (bottom left) and in the PML-PML case (bottom right).

Fig. 9 .

 9 Fig. 9. Potential flow around a submarine on the left (the background color corresponds to the norm of v/ v ∞). Diffracted field from a Dirac point source on the right.

  [START_REF] Bethuel | Travelling waves for the Gross-Pitaevskii equation I[END_REF] and ω = 200. For the PML, we took = 0.02 and σ x = σ y = 15. Also, we assume the flow to be constant and horizontal in thePML region v = ( v ext , 0) T ,although this is not exact: the flow is almost horizontal and constant. The domain is decomposed into 5 subdomains defined by (3.24) α 1 = 0, α i = 0.2(i -1) -0.015 for i ∈ {2, • • • , 5}, β 5 = 1, β i = 0.2i + 0.015 for i ∈ {1, • • • , 4},

Fig. 10 .a 2 A - 1 ,

 1021 Fig. 10. GMRES residual versus iterations for the three types of transmission conditions (ABC), (ABC PML) and (OPT).

This manuscript is for review purposes only.

This value is in fact not realistic since in water the sound speed is much larger than in air. As far as we know, the fastest speed one can reach in water is around mach 0.075 with supercavitation. Nevertheless, this very high speed is only obtained close to corners of the submarine, where the velocity increases a lot. In the rest of the domain, the speed is more realistic.

Acknowledgments. This research was supported by the Swiss National Science Foundation.