# Optimized Schwarz domain decomposition methods for scalar and vector Helmholtz equations

Xavier Antoine, Christophe Geuzaine

► To cite this version:

Xavier Antoine, Christophe Geuzaine. Optimized Schwarz domain decomposition methods for scalar and vector Helmholtz equations. D. Lahaye, J. Tang and K. Vuik. Modern Solvers for Helmholtz Problems, Springer International Publishing, pp.189-213, 2016, Lecture notes in geosystems mathematics and computing, 10.1007/978-3-319-28832-1_8. hal-01276808

HAL Id: hal-01276808

https://hal.science/hal-01276808

Submitted on 23 Feb 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Optimized Schwarz Domain Decomposition Methods
for Scalar and Vector Helmholtz Equations

X. Antoine<fnref n="1" /> and C. Geuzaine<fnref n="2" />

February 20, 2016

Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Scalar Helmholtz Equation: Acoustic Waves</td><td>2</td></tr><tr><td>2.1</td><td>Domain Decomposition and Transmission Operators</td><td>3</td></tr><tr><td>2.2</td><td>Convergence Analysis on a Model Problem</td><td>7</td></tr><tr><td>2.3</td><td>Weak Formulations</td><td>9</td></tr><tr><td>3</td><td>Vector Helmholtz Equation: Electromagnetic Waves</td><td>11</td></tr><tr><td>3.1</td><td>Domain Decomposition and Transmission operators</td><td>12</td></tr><tr><td>3.2</td><td>Convergence Analysis on a Model Problem</td><td>14</td></tr><tr><td>3.3</td><td>Weak Formulations</td><td>16</td></tr><tr><td>4</td><td>Numerical Implementation</td><td>17</td></tr></table>

Abstract

In this chapter we review Schwarz domain decomposition methods for scalar and vector Helmholtz equations, with a focus on the choice of the associated transmission conditions between the subdomains. The methods are analyzed in both acoustic and electromagnetic settings, and generic weak formulations directly amenable to finite element discretization are presented. An open source solver along with ready-to-use examples is freely available online for further testing.

# 1 Introduction

Solving high-frequency time-harmonic wave problems is a very challenging problem, en-
countered in many physical applications, from acoustic noise propagation to seismology
and geophysical exploration to electromagnetic radiation. Among the various approaches

<footnote n="1">¹Institut Elie Cartan de Lorraine, Université de Lorraine, Inria Nancy-Grand Est EPI SPHINX, F-54506 Vandoeuvre-lès-Nancy Cedex, France, xavier.antoine@univ-lorraine.fr.</footnote>

<footnote n="2">²Université de Liège, Institut Montefiore B28, B-4000 Liège, Belgium, cgeuzaine@ulg.ac.be.</footnote>
for numerical simulation, the Finite Element Method (FEM) with an Absorbing Boundary Condition (ABC) or a Perfectly Matched Layer (PML) is well suited for tackling complex geometrical configurations and heterogeneous media. The brute-force application of the FEM in the high-frequency regime however requires the solution of extremely large, complex-valued and possibly indefinite linear systems [39]. Direct sparse solvers do not scale well for such large-size problems, and Krylov subspace iterative solvers exhibit slow convergence or diverge, while efficiently preconditioning proves difficult [24]. Domain decomposition methods provide an alternative, iterating between subproblems of smaller sizes, amenable to sparse direct solvers [49].

In [36], Lions introduced a converging Schwarz domain decomposition method without overlap for the Laplace equation by using Fourier-Robin boundary conditions on the interfaces instead of the standard Dirichlet or Neumann continuity conditions. For scalar or vector Helmholtz equations, these methods need to be adapted to lead to converging iterative algorithms. The first developments in this direction were introduced by Després [13, 14], who used simple impedance boundary conditions on the interfaces. A great variety of more general impedance conditions has been proposed since these early works, leading to so-called optimized Schwarz domain decomposition methods for time-harmonic wave problems [1, 9, 10, 11, 14, 15, 16, 19, 20, 25, 27, 43, 44, 45]. These methods can be used with or without overlap between the subdomains, and their convergence rate strongly depends on the transmission condition. Optimal convergence is obtained by using as transmission condition on each interface the non-local Dirichlet-to-Neumann (DtN) map [42] related to the complementary of the subdomain of interest [40, 41]. For acoustic waves, this DtN map links the normal derivative and the trace of the acoustic pressure on the interface. For electromagnetic waves, it links the magnetic and the electric surface currents (and is referred to in this case as the Magnetic-to-Electric, or MtE, map) [19]. However, using the DtN leads to a very expensive numerical procedure in practice, as this operator is non-local. Practical algorithms are thus based on local approximations of these operators, both for the acoustic case [13, 9, 10, 11, 27] and the electromagnetic one [1, 14, 15, 16, 20, 21, 43, 44, 45]. Recently, PMLs have also been used for this same purpose [23, 47, 51, 52].

In this chapter we provide a concise review of the most common transmission operators for optimized Schwarz methods applied to time-harmonic acoustic and electromagnetic wave problems, with the corresponding mathematical background. We analyze the behavior of these transmission operators on a model problem and derive generic weak formulations in view of their implementation in finite element codes. All the formulations are readily available for testing on several acoustic and electromagnetic cases using the open source GetDDM environment (http://onelab.info/wiki/GetDDM) [33, 48], based on the finite element solver GetDP (http://getdp.info) [17, 18, 28] and the mesh generator Gmsh (http://gmsh.info) [31, 32].

## 2 Scalar Helmholtz Equation: Acoustic Waves

Let $\Omega^-$ be an open subset of $\mathbb{R}^d$ ($d = 1, 2, 3$) with boundary $\Gamma := \partial\Omega^-$. The exterior domain of propagation is the complementary connected set defined by $\Omega^+ = \mathbb{R}^d \setminus \bar{\Omega}^-$. When considering a time-harmonic incident wave $u^{\text{inc}}$, the obstacle $\Omega^-$ creates a complex-
valued scattered field $u$ which is solution of the following problem

$$ \left\{ \begin{array}{ll} (\Delta + k^2)u = 0 & \text{in } \Omega^+, \\ u = -u^{\text{inc}} & \text{on } \Gamma, \\ u \text{ outgoing,} & \end{array} \right. \qquad (1) $$

fixing the time dependence under the form $e^{-i\omega t}$. The Laplacian operator is $\Delta = \sum_{i=1}^d \partial_{x_i}^2$ and the real-valued strictly positive wavenumber is given by $k = \omega/c$ (where $c = c(\mathbf{x})$ is the local speed of sound in the propagation medium). We denote by $\mathbf{a} \cdot \bar{\mathbf{b}}$ the inner product between two complex-valued vectors $\mathbf{a}$ and $\mathbf{b}$ in $\mathbb{C}^3$. We designate by $\bar{z}$ the complex conjugate of $z \in \mathbb{C}$ and the associated norm is $||\mathbf{a}|| := \sqrt{\mathbf{a} \cdot \bar{\mathbf{a}}}$. In this chapter, we fix a Dirichlet boundary condition on $\Gamma$ corresponding to the sound-soft obstacle case. Nevertheless, any other condition can be studied like e.g. for a Neumann, Fourier or even for a penetrable obstacle. The outgoing condition at infinity, better known as Sommerfeld radiation condition ($i$ being the square root of $-1$), is added

$$ \lim_{\|\mathbf{x}\| \to \infty} \|\mathbf{x}\|^{\frac{d-1}{2}} \left( \nabla u \cdot \frac{\mathbf{x}}{\|\mathbf{x}\|} - i\kappa u \right) = 0. $$

This allows to prove that the solution to (1) is unique. In addition, this translates the property that the scattered field $u$ is directed from $\Omega^-$ to infinity.

To numerically compute the solution to problem (1) by using e.g. the finite element method, $\Omega^+$ has to be truncated. This can be realized for example by introducing a Perfectly Matched Layer (PML) [7, 12] or a fictitious boundary $\Gamma^\infty$ with an Absorbing Boundary Condition (ABC) [6, 22] (see e.g. [4] for a review). If we consider an ABC on a fictitious boundary, we have to compute a field $\hat{u}$ approximating $u$ on the finite domain $\Omega$ with boundary $\Gamma^\infty \cup \Gamma$. After merging the notations $\hat{u}$ and $u$ for simplicity, the problem to be solved is

$$ \left\{ \begin{array}{ll} (\Delta + k^2)u = 0 & \text{in } \Omega, \\ u = -u^{\text{inc}} & \text{on } \Gamma, \\ \partial_{\mathbf{n}}u + \mathcal{B}u = 0 & \text{on } \Gamma^{\infty}, \end{array} \right. \qquad (2) $$

where the unit normal vector $\mathbf{n}$ is directed outside $\Omega$ (and thus inside $\Omega^-$ on $\Gamma$). The simplest local ABC, i.e., the Sommerfeld radiation condition at finite distance (zeroth-order condition), is obtained by setting

$$ \mathcal{B}u = -i\kappa u. \qquad (3) $$

The extension to more accurate ABCs or PMLs is standard [34].

## 2.1 Domain Decomposition and Transmission Operators

Let us consider now that $\Omega$ is decomposed into $N_{\text{dom}}$ disjoint subdomains $\Omega_i$ (the sub-structures) without overlap. For every $i = 0, \dots, N_{\text{dom}} - 1$, we set $\Gamma_i = \Gamma \cap \partial\Omega_i$, $\Gamma_i^\infty = \Gamma^\infty \cap \partial\Omega_i$, and, for $j = 0, \dots, N_{\text{dom}} - 1, j \neq i$, we introduce the transmission boundary $\Sigma_{ij} = \Sigma_{ji} = \bar{\partial}\Omega_i \cap \partial\Omega_j$. To simplify, let $D := \{0, \dots, N_{\text{dom}} - 1\}$ be the set of indices of the subdomains, and for $i \in D$, let $D_i := \{j \in D \text{ such that } j \neq i \text{ and } \Sigma_{ij} \neq \emptyset\}$ be the set of indices of the subdomains sharing at least a point with $\Omega_i$ (such a domain
is said to be connected to $\Omega_i$). Finally, for all $i \in D$, the unit normal $n_i$ is directed into the exterior of $\Omega_i$ and thus inside the obstacle $\Omega^-$ (if $\Gamma_i \neq \emptyset$).

Then the additive Schwarz domain decomposition method follows the steps at iteration $n+1$

1. For all $i \in D$, compute $u_i^{n+1}$ solution to the boundary-value problem

$$ egin{cases} (\Delta + k^2)u_i^{n+1} = 0 & \text{in } \Omega_i, \\ u_i^{n+1} = -u_i^{\text{inc}} & \text{on } \Gamma_i, \\ \partial_{n_i} u_i^{n+1} + \mathcal{B} u_i^{n+1} = 0 & \text{on } \Gamma_i^\infty, \\ \partial_{n_i} u_i^{n+1} + \mathcal{S} u_i^{n+1} = g_{ij}^n & \text{on } \Sigma_{ij}, \quad \forall j \in D_i. \end{cases} \quad (4) $$

2. For all $i \in D$ and $j \in D_i$, update the interface unknowns with respect to the relation

$$ g_{ji}^{n+1} = -\partial_{n_i} u_i^{n+1} + \mathcal{S} u_i^{n+1} = -g_{ij}^n + 2\mathcal{S} u_i^{n+1}, \quad \text{on } \Sigma_{ij}. \quad (5) $$

The operator $\mathcal{S}$ is a transmission operator that will be described later. A more compact writing of the $(n+1)^{th}$ iteration is

1. For all $i \in D$, compute the volume solution $u_i^{n+1}$ of problem (4), which is written here as $u_i^{n+1} = \mathcal{V}_i(u_i^{\text{inc}}, g^n)$, where $g^n = (g_{ji}^n)_{i \in D, j \in D_i}$ is the vector that collects all the contributions related to the interface unknowns.

2. For all $i \in D$ and $j \in D_i$, update the surface fields $g_{ji}^{n+1}$ following relation (5). This is written as $g_{ji}^{n+1} = \mathcal{T}_{ji}(g_{ij}^n, u_i^{n+1})$ in the sequel of the paper.

In the boundary-value problem (4), only the case of Dirichlet sources is considered; however, any kinds such as volume sources could be handled similarly in the algorithm. These sources are called physical sources in contrast with the artificial sources $g_{ij}^n$ related to the transmission boundaries.

The algorithm described by (4) and (5) can be understood as a Jacobi iteration for a linear operator equation. For every $n \in \mathbb{N}$, the field $u_i^{n+1}$ can be decomposed by linearity as $u_i^{n+1} = v_i^{n+1} + \tilde{u}_i^{n+1}$, with

$$ v_i^{n+1} = \mathcal{V}_i(u_i^{\text{inc}}, 0) \quad \text{and} \quad \tilde{u}_i^{n+1} = \mathcal{V}_i(0, g^n). \quad (6) $$

The function $v_i^{n+1}$ does not depend on the iteration $n$ and can be written as $v_i := v_i^n$, $\forall n \in \mathbb{N}, \forall i \in D$. Therefore, Equation (5) can be written

$$ g_{ji}^{n+1} = \mathcal{T}_{ji}(g_{ij}^n, u_i^{n+1}) = \mathcal{T}_{ji}(g_{ij}^n, \tilde{u}_i^{n+1}) + 2\mathcal{S}v_i, \quad \text{on } \Sigma_{ij}. \quad (7) $$

Let us define the vector $b = (b_{ji})_{i \in D, j \in D_i}$, with $b_{ji} = 2(\mathcal{S}v_i)|_{\Sigma_{ij}}$, and $\mathcal{A}: g^n \mapsto \mathcal{A}g^n$ as the operator such that

$$ \forall i \in D \quad \left\{ \begin{array}{l} \tilde{u}_i^{n+1} = \mathcal{V}_i(0, g^n), \\ (\mathcal{A}g^n)_{ji} = \mathcal{T}_{ji}(g_{ij}^n, \tilde{u}_i^{n+1}), \quad \forall j \in D_i. \end{array} \right. \quad (8) $$

One iteration of the domain decomposition method writes

$$ g^{n+1} = \mathcal{A}g^n + b. \quad (9) $$
This can be interpreted as an iteration of the Jacobi method for solving the system

$$ (\mathcal{I} - \mathcal{A})g = b, \qquad (10) $$

where the identity operator is $\mathcal{I}$. An interesting consequence of (10) is that any iterative linear solver can be used for solving the equation. For example, Krylov subspace methods can be applied such as GMRES [46]. When a Krylov subspace solver is used, the resulting method is called a substructured preconditioner [26].

An important remark is that the iteration unknowns in (9), (10) are the surface quantities $g$ and not the volume unknowns $u$. To get the volume quantities from the surface unknowns, $u_i = \mathcal{V}_i(u^{\text{inc}}, g)$ needs to be solved on every subdomain $\Omega_i$. Algorithm 1 summarizes the Schwarz method with Krylov solver.

**Algorithm 1** Schwarz algorithm with Krylov solver.

1. Compute the right-hand side $b$

$$ \left\{ \begin{array}{ll} \forall i \in D, & v_i = \mathcal{V}_i(u^{\text{inc}}, 0), \\ \forall i \in D, \forall j \in D_i, & b_{ji} = \mathcal{T}_{ji}(0, v_i). \end{array} \right. $$

2. Solve the following system $(\mathcal{I} - \mathcal{A})g = b$ iteratively by using a Krylov subspace solver, where the operator $\mathcal{A}$ is given by (8).

3. At convergence, compute the solution: $\forall i \in D$, $u_i = \mathcal{V}_i(u^{\text{inc}}, g)$.

The convergence rate of the iterative solver is strongly related to the choice of the transmission operator $\mathcal{S}$ [10]. The so-called Dirichlet-to-Neumann (DtN) map for the complement of each subdomain [40, 41] appears as being optimal. Unfortunately, this operator is nonlocal and consequently costly to use in an iterative solver. An alternative approach consists in using local approximations based on polynomial or rational approximations of the total symbol of the surface DtN operator in the free-space, or a volume representation through PMLs. We detail below four specific examples which are also implemented in GetDDM for a generic transmission boundary $\Sigma$

* *Evanесcent Modes Damping Algorithm* [9, 11]:

$$ S_{\text{IBC}(\chi)}u = (-\imath k + \chi)u, $$

where $\chi$ is a real-valued constant. This zeroth-order polynomial approximation is a generalization of the well-known Després condition [14], which corresponds to $\chi = 0$. We will denote this family of impedance transmission conditions as IBC($\chi$) in the sequel of the paper.

* *Optimized second-order transmission condition* [27]:

$$ S_{\text{GIBC}(a,b)}u = au + b\Delta_{\Sigma}u, \qquad (11) $$

where $\Delta_{\Sigma}$ designates the Laplace-Beltrami operator on $\Sigma$, and $a$ and $b$ are two complex-valued numbers computed by solving a min-max optimization problem
involving the rate of convergence (spectral radius) of the iteration operator. At the symbol level, this condition yields a second-order polynomial approximation of the DtN symbol. In the following, this family of generalized impedance transmission conditions is denoted by GIBC(a, b). A zeroth-order optimized condition can be built similarly.

• Padé-localized square-root transmission condition [10]:

$$S_{\text{GIBC}(N_p, \alpha, \varepsilon)} u = -\iota k C_0 u - \iota k \sum_{\ell=1}^{N_p} A_{\ell} \text{div}_{\Sigma} \left( \frac{1}{k_{\varepsilon}^2} \nabla_{\Sigma} \right) \left( \mathcal{I} + B_{\ell} \text{div}_{\Sigma} \left( \frac{1}{k_{\varepsilon}^2} \nabla_{\Sigma} \right) \right)^{-1} u, \quad (12)$$

setting

$$k_{\varepsilon} = k + i\varepsilon. \quad (13)$$

The complex-valued coefficients $C_0$, $A_\ell$ and $B_\ell$ are

$$C_0 = e^{i\alpha/2} R_{N_p} (e^{-i\alpha} - 1), \quad A_{\ell} = \frac{e^{-\frac{i\alpha}{2}} a_{\ell}}{(1 + b_{\ell}(e^{-i\alpha} - 1))^2}, \quad B_{\ell} = \frac{e^{-i\alpha} b_{\ell}}{1 + b_{\ell}(e^{-i\alpha} - 1)}. \quad (14)$$

The parameter $\alpha$ is a rotation angle in the complex plane (usually taken as $\pi/4$) and $R_{N_p}$ are the standard real-valued Padé approximation of order $N_p$ of $\sqrt{1+z}$

$$R_{N_p}(z) = 1 + \sum_{\ell=1}^{N_p} \frac{a_{\ell}z}{1 + b_{\ell}z},$$

with

$$a_{\ell} = \frac{2}{2N_p + 1} \sin^2 \left( \frac{\ell \pi}{2N_p + 1} \right) \quad \text{and} \quad b_{\ell} = \cos^2 \left( \frac{\ell \pi}{2N_p + 1} \right). \quad (15)$$

This transmission condition is a complex-valued rational approximation [37] of the nonlocal pseudodifferential operator

$$S_{\text{GIBC}(\text{sq}, \varepsilon)} u = -\iota k \sqrt{1 + \text{div}_{\Sigma} \left( \frac{1}{k_{\varepsilon}^2} \nabla_{\Sigma} \right)} u.$$

Fixing $\varepsilon = 0$ leads to the principal symbol of the exact DtN operator for the half-space. The introduction of the parameter $\varepsilon$ regularizes this operator to model glancing rays at the surface of a curved interface. An optimal choice of $\varepsilon$ is explained below in section 2.2. In the sequel of the paper, we denote this family of generalized impedance transmission conditions as GIBC($N_p, \alpha, \varepsilon$) and GIBC(sq, $\varepsilon$), respectively.

• PML transmission condition [23, 47, 51, 52]: The operator $S_{\text{PML}(\sigma)}$ is constructed by appending a layer $\Omega^{\text{PML}}$ to the transmission interface, in which a PML transformation with absorption profile $\sigma$ is applied. For example, in cartesian coordinates, the singular profile

$$\sigma(x_{\text{PML}}) = \frac{1}{k(x_{\text{PML}} - \delta)}$$

can be used, where $\delta$ corresponds to the thickness of the PML layer and $x_{\text{PML}}$ is the local coordinate inside the PML [8, 38].
All these methods are referred to as optimized Schwarz domain decomposition methods. Note that GIBC($N_p$, $\alpha$, $\varepsilon$) and PML($\sigma$) have in common that they introduce additional surface/volume unknowns, whereas the other two transmission conditions do not. Also, the first three transmission conditions can be formulated explicitly through sparse surface equations (see e.g. the weak formulations (20)–(25) below), while a sparse formulation of the PML transmission condition requires a volume representation (see e.g. (26)–(27)), a surface representation being dense [50].

## 2.2 Convergence Analysis on a Model Problem

To study the impact of the various transmission conditions on the convergence of DDM, we analyze the model problem depicted in Figure 1 which couples two subdomains: a disk-shaped *bounded* subdomain $\Omega_1$ of radius $R_0$ and an *unbounded* domain $\Omega_0 = \mathbb{R}^2 \setminus \Omega_1$:

$$ \Omega_0 := \{\mathbf{x} \in \mathbb{R}^2, |\mathbf{x}| > R_0\}, \quad \Omega_1 := \{\mathbf{x} \in \mathbb{R}^2, |\mathbf{x}| < R_0\}, \qquad (16) $$

with $\partial\Omega_0 = \partial\Omega_1 = \Sigma$. We analyze the spectral properties of the iteration operator $\mathcal{A}$ obtained from the domain decomposition algorithm coupling these two subdomains. Understanding the coupling of a curved bounded and unbounded subdomains allows us to clarify the main properties that one could not be analyzed by considering two bounded (e.g. a square domain divided in two) or two unbounded (e.g. two half-planes) subdomains. The considered model problem essentially contains the main features arising when solving exterior scattering problems in homogeneous media. It is thus not directly applicable to the PML-based transmission conditions, which introduce a fictitious heterogeneous medium, even for a radial profile.

<figure><img src="image_2.png" /><figcaption>Figure 1: Model problem with two subdomains and a circular interface.</figcaption></figure>

For this problem, the iteration operator $\mathcal{A}$ can be expanded as $\mathcal{A} = \sum_{m=-\infty}^{+\infty} \mathcal{A}_m e^{im\theta}$. We report in Figure 2 the modal spectral radius $\rho(\mathcal{A}_m)$ with respect to the Fourier mode $m$ for the transmitting boundary conditions IBC(0), IBC($k/2$), GIBC($a, b$) and GIBC(sq, 0). We fix $k = 6\pi$, $R_0 = 1$ and the maximal number of modes is set to $m^{\max} = [10kR_0]$ (where $[10kR_0]$ denotes the integer part of $10kR_0$). Clearly, IBC(0) leads to a spectral radius equal to 1 for the evanescent modes, which is improved by IBC($k/2$)—for which the radius of convergence is always strictly less than one. Using GIBC($a, b$) further improves
<figure><img src="image_3.png" /><figcaption>Figure 2: Spectral radius of the modal iteration operator $\mathcal{A}_m$ vs. the Fourier mode $m$.</figcaption></figure>

over IBC($k/2$), particularly for large spatial modes $m$. We recall here that the GIBC($a, b$) method is based on optimizing the coefficients $a$ and $b$ in relation (11) according to a min-max problem posed in the Fourier space [10, 27]. For the square-root transmission condition with $\varepsilon = 0$ (GIBC(sq,0)), we clearly observe an optimal convergence rate in the evanescent part of the spectrum. We also see a significant improvement over the IBC(0), IBC($\chi$) and GIBC($a, b$) algorithms on the propagating modes. The damping parameter $\varepsilon$ can be optimized to further improve the spectrum of the iteration operator corresponding to the modes in the transition zone. The optimization problem can be formulated as a min-max problem: find $\varepsilon_{opt} > 0$ such that it minimizes the spectral radius $\rho(\mathcal{A}_m)$ of the iteration operator (associated with GIBC(sq, $\varepsilon$)) for the mode $m \in \mathbb{Z}$ where it is maximal. Mathematically, this leads to solving the problem

$$ \rho^{\text{sq}, \varepsilon_{\text{opt}}} = \min_{\varepsilon \in \mathbb{R}^+} \left( \max_{m \in \mathbb{Z}} |\rho(\mathcal{A}_m)| \right), \qquad (17) $$

resulting in the estimate $\varepsilon_{opt} = 0.4k^{1/3}\mathcal{H}^{2/3}$ [10] of the optimal value of the damping parameter, where $\mathcal{H}$ is the mean curvature on $\Sigma$. We see in Figure 2 that the spectral radius of the iteration operator is indeed locally minimized for $\varepsilon_{opt}$.

Fast convergence of the GMRES solver is known to be strongly linked to the existence of eigenvalues clustering of the operator to solve, i.e. $(I - \mathcal{A})$ in our case. We report in Figure 3 (left) the spectrum of the iteration operator for IBC(0), IBC($k/2$), GIBC($a, b$) and GIBC(sq, $\varepsilon_{opt}$) (again for $kR_0 = 6\pi$ and $m^{\max} = [10kR_0]$). For all transmission operators, the spectrum lies in the right half-plane, which makes the GMRES converging. Nevertheless, many eigenvalues spread out in the complex plane for IBC(0). A slightly better clustering occurs for IBC($k/2$) and GIBC($a, b$), while there is an excellent clustering of the eigenvalues for GIBC(sq, $\varepsilon_{opt}$). Most particularly, only a few eigenvalues associated with the propagating modes do not cluster but are very close to (1, 0). In addition, the eigenvalues linked to the evanescent modes seem to cluster at (1, 0). The eigenvalues clustering for the evanescent modes can be shown in numerical experiments
<figure><img src="image_4.png" /><figcaption>Figure 3: Left: Eigenvalues distribution in the complex plane for ($I - A$) and different transmission operators. Right: Eigenvalue distribution in the complex plane for the exact and Padé-localized square-root transmission operator of order 4.</figcaption></figure>

to lead to a quasi-optimal GMRES convergence rate that is independent of the density of discretization points per wavelength $n_λ$ [10].

As said before, the square-root operator (2.1) is a first-order nonlocal pseudodifferential operator. Therefore, it is impractical in a finite element setting since it would lead to consider full complex-valued matrices at the transmission interfaces. Fortunately, a localization process of this operator can be efficiently realized and based on partial differential (local) operators to have a sparse matrix representation. In [3, 5, 37], this is done by using a rotating branch-cut approximation of the square-root and next applying complex Padé approximants of order $N_p$, leading to the transmission operator (12). We report in Figure 3 (right) the spectrum of the modal iteration operators GIBC(sq, $ε_{opt}$) and GIBC(4, $π/4$, $ε_{opt}$). As already noticed, there is an almost perfect clustering of the eigenvalues for GIBC(sq, $ε_{opt}$). As expected, the larger $N_p$, the better the approximation of the spectrum of the square-root. Moreover, $N_p$ allows to adjust the spectrum accuracy for large modes $m$ (evanescent modes which numerically correspond to mesh refinement in a finite element context). Numerical simulations show that in practice relatively small values of $N_p$ ($N_p = 2, 4, 8$) give optimal convergence results.

## 2.3 Weak Formulations

For the finite element approximation, we consider some variational formulations. Two kinds of PDEs are involved when using optimized Schwarz methods: firstly, a volume system (in the present case, the scalar Helmholtz equation) given by $\mathcal{Y}_i$, and, secondly, a surface system on the transmission interfaces, fixed by $\mathcal{T}_{ji}$. The variational formulations are first provided for a general transmission operator $S$. To simplify the presentation, we consider the situation where no contribution comes from $\partial\Sigma_{ij}$ through an integration by parts. However, in some cases (e.g. when $\Sigma_{ij} \cap \Gamma^\infty \neq \emptyset$), a special attention must be
directed towards the inclusion of these terms into the variational formulations.

Without loss of generality, we only detail the case of a particular subdomain $\Omega_i$, for $i \in D$, without incident wave contribution (i.e. homogeneous Dirichlet boundary condition). We consider the general setting where PML layers $\Omega_i^{\text{PML}} = \cup_{j \in D_i} \Omega_{ij}^{\text{PML}}$ are potentially appended to the artificial interfaces $\Sigma_{ij}$, and define $\Omega_i^* := \Omega_i \cup \Omega_i^{\text{PML}}$. In what follows, the space $H^1(\Omega_i^*) := \{\tilde{u}_i \in L^2(\Omega_i^*) \text{ such that } \nabla \tilde{u}_i \in (L^2(\Omega_i^*))^3\}$ is the classical Sobolev space and $H_0^1(\Omega_i^*)$ is the space of functions $\tilde{u}_i \in H^1(\Omega_i^*)$ such that $\tilde{u}_i|_{\Gamma_i} = 0$, which slightly differs from its usual definition (the Dirichlet condition is here set only on part of $\partial \Omega_i^*$). Then,

• the volume PDE $\tilde{u}_i^{n+1} = \mathcal{V}_i(0, g^n)$ has the following weak formulation

$$ \left\{ \begin{array}{l} \text{Find } \tilde{u}_i^{n+1} \text{ in } H_0^1(\Omega_i^*) \text{ such that, for every } \tilde{u}_i' \in H_0^1(\Omega_i^*): \\ \displaystyle \int_{\Omega_i} \nabla \tilde{u}_i^{n+1} \cdot \nabla \tilde{u}_i' \, d\Omega_i - \int_{\Omega_i} k^2 \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Omega_i + \int_{\Gamma_i^\infty} \mathcal{B} \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Gamma_i^\infty \\ \displaystyle \qquad + \sum_{j \in D_i} \int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij} = \sum_{j \in D_i} \int_{\Sigma_{ij}} g_{ij}^n \tilde{u}_i' \, d\Sigma_{ij}, \end{array} \right. \quad (18) $$

• and the surface PDE $g_{ji}^{n+1} = \mathcal{T}_{ji}(g_{ij}^n, \tilde{u}_i^{n+1})$ has the following one:

$$ \left\{ \begin{array}{l} \text{Find } g_{ji}^{n+1} \text{ in } H^1(\Sigma_{ij}) \text{ such that, for every } g_{ji}' \in H^1(\Sigma_{ij}): \\ \displaystyle \int_{\Sigma_{ij}} g_{ji}^{n+1} g_{ji}' \, d\Sigma_{ij} = - \int_{\Sigma_{ij}} g_{ij}^n g_{ji}' \, d\Sigma_{ij} + 2 \int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} g_{ji}' \, d\Sigma_{ij}. \end{array} \right. \quad (19) $$

Depending on the choice of the transmission operator $\mathcal{S}$, the quantities $\int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij}$ and $\int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} g_{ji}' \, d\Sigma_{ij}$ write as follows:

• IBC($\chi$):

$$ \int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} (-\iota k + \chi) \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij}; \quad (20) $$

$$ \int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} g_{ji}' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} (-\iota k + \chi) \tilde{u}_i^{n+1} g_{ji}' \, d\Sigma_{ij}. \quad (21) $$

• GIBC($a, b$):

$$ \int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} a \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij} - \int_{\Sigma_{ij}} b \nabla \tilde{u}_i^{n+1} \cdot \nabla \tilde{u}_i' \, d\Sigma_{ij}; \quad (22) $$

$$ \int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} g_{ji}' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} a \tilde{u}_i^{n+1} g_{ji}' \, d\Sigma_{ij} - \int_{\Sigma_{ij}} b \nabla \tilde{u}_i^{n+1} \cdot \nabla g_{ji}' \, d\Sigma_{ij}. \quad (23) $$

• GIBC($N_p, \alpha, \varepsilon$):

$$ \int_{\Sigma_{ij}} \mathcal{S} \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij} := -\iota k C_0 \int_{\Sigma_{ij}} \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij} + \iota k \sum_{\ell=1}^{N_p} A_\ell \int_{\Sigma_{ij}} \frac{1}{k_\varepsilon^2} \nabla_{\Sigma_{ij}} \varphi_\ell \cdot \nabla_{\Sigma_{ij}} \tilde{u}_i' \, d\Sigma_{ij}, \quad (24) $$
where, for every $\ell = 1, \dots, N_p$, the function $\varphi_\ell$ is obtained through the resolution of

$$ egin{cases} 
\text{Find } \varphi_\ell \text{ in } H^1(\Sigma_{ij}) \text{ such that, for every } \varphi'_\ell \in H^1(\Sigma_{ij}): \\ 
- \int_{\Sigma_{ij}} \tilde{u}_i^{n+1} \varphi'_\ell \, d\Sigma_{ij} - B_\ell \int_{\Sigma_{ij}} \frac{1}{k_\varepsilon^2} \nabla_{\Sigma_{ij}} \varphi_\ell \cdot \nabla_{\Sigma_{ij}} \varphi'_\ell \, d\Sigma_{ij} + \int_{\Sigma_{ij}} \varphi_\ell \cdot \varphi'_\ell \, d\Sigma_{ij} = 0; \\ 
\int_{\Sigma_{ij}} S \tilde{u}_i^{n+1} g'_{ji} \, d\Sigma_{ij} := -\iota k C_0 \int_{\Sigma_{ij}} \tilde{u}_i^{n+1} g'_{ji} \, d\Sigma_{ij} - \iota k \sum_{\ell=1}^{N_p} \frac{A_\ell}{B_\ell} \int_{\Sigma_{ij}} (\tilde{u}_i^{n+1} - \varphi_\ell) g'_{ji} \, d\Sigma_{ij}. 
\end{cases} 
\quad (25) $$

• PML($\sigma$):

$$ \int_{\Sigma_{ij}} S \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Sigma_{ij} := \int_{\Omega_{ij}^{\text{PML}}} D \nabla \tilde{u}_i^{n+1} \cdot \nabla \tilde{u}_i' \, d\Omega_{ij}^{\text{PML}} - \int_{\Omega_{ij}^{\text{PML}}} k^2 E \tilde{u}_i^{n+1} \tilde{u}_i' \, d\Omega_{ij}^{\text{PML}}, \quad (26) $$

$$ \int_{\Sigma_{ij}} S \tilde{u}_i^{n+1} g'_{ji} \, d\Sigma_{ij} := \int_{\Omega_{ij}^{\text{PML}}} D \nabla \tilde{u}_i^{n+1} \cdot \nabla g'_{ji} \, d\Omega_{ij}^{\text{PML}} - \int_{\Omega_{ij}^{\text{PML}}} k^2 E \tilde{u}_i^{n+1} g'_{ji} \, d\Omega_{ij}^{\text{PML}}, \quad (27) $$

where $D = \text{diag}(\frac{1}{\gamma_x}, \gamma_x, \gamma_x)$ and $E = \gamma_x$, with $\gamma_x(x_{\text{PML}}) = 1 + \frac{i}{\omega}\sigma_x(x_{\text{PML}})$, that is, we consider a 1D PML with an absorption function that grows only in the direction normal to the interface. In (27) the domain of definition of the test functions $g'_{ji}$ on $\Sigma_{ij}$ is extended to the neighboring PML layer $\Omega_{ij}^{\text{PML}}$, effectively resulting at the discrete level in the integration of the functions associated with the nodes of the interface in the layer of volume elements connected to the interface.

## 3 Vector Helmholtz Equation: Electromagnetic Waves

We now consider the case of an incident electromagnetic wave $E^{\text{inc}}$ illuminating a perfectly conducting obstacle $\Omega^{-}$ with boundary $\Gamma$, in a three dimensional medium. The scattered electric field $E$ is solution to the following exterior electromagnetic scattering problem:

$$ \left\{ \begin{array}{l} \text{curl curl } E - k^2 E = 0, \quad \text{in } \Omega^+, \\ \gamma^T(E) = -\gamma^T(E), \quad \text{on } \Gamma, \\ \lim_{\|\mathbf{x}\| \to \infty} \|\mathbf{x}\| \left( \frac{\mathbf{x}}{\|\mathbf{x}\|} \times \text{curl } E + \iota k E \right) = 0, \end{array} \right. \quad (28) $$

where $k := 2\pi/\lambda$ is again the wavenumber and $\lambda$ the wavelength, $\mathbf{n}$ is the outward unit normal to $\Omega^+$ (thus, inward to the obstacle) and $\gamma^T$ is the tangential component trace operator

$$ \gamma^T : \mathbf{v} \mapsto \mathbf{n} \times (\mathbf{v} \times \mathbf{n}). $$

The curl operator is defined by $\text{curl } \mathbf{a} := \nabla \times \mathbf{a}$, for a complex-valued vector field $\mathbf{a} \in \mathbb{C}^3$, and the notation $\mathbf{a} \times \mathbf{b}$ designates the cross product between two complex-valued vectors $\mathbf{a}$
and **b**. The last equation of system (28) is the so-called Silver-Müller radiation condition at infinity, which provides the uniqueness of the solution to the scattering boundary-value problem (28).

As in the acoustic case, solving (28) numerically with a volume discretization method requires the truncation of the exterior propagation domain with a PML or with an ABC on a fictitious boundary $\Gamma^\infty$ surrounding $\Omega^-$. For an ABC the problem to be solved is then defined on the bounded domain $\Omega$, with boundaries $\Gamma$ and $\Gamma^\infty$:

$$ \left\{ \begin{array}{ll} \text{curl } \text{curl } \mathbf{E} - k^2 \mathbf{E} = 0, & \text{in } \Omega, \\ \gamma^T(\mathbf{E}) = -\gamma^T(\mathbf{E}), & \text{on } \Gamma, \\ \gamma^t(\text{curl } \mathbf{E}) + \mathcal{B}(\gamma^T(\mathbf{E})) = 0, & \text{on } \Gamma^\infty, \end{array} \right. \qquad (29) $$

with $\gamma^t$ the tangential trace operator:

$$ \gamma^t : \mathbf{v} \mapsto \mathbf{n} \times \mathbf{v}. $$

As above, the unit normal **n** is outwardly directed to $\Omega$ and, to simplify, the solution of the above problem is still designated by **E**. The operator $\mathcal{B}$ is an approximation of the Magnetic-to-Electric (MtE) operator. The well-known Silver-Müller ABC at finite distance is obtained with $\mathcal{B} = \iota k$, similar to (3) for acoustics modulo the sign (due to the trace operator definitions). The extension to more accurate ABCs or PMLs is standard.

## 3.1 Domain Decomposition and Transmission operators

The optimized Schwarz domain decomposition without overlap for the Maxwell problem (29) can be set up in exactly the same way as for the scalar Helmholtz equation. The domain $\Omega$ is decomposed as described in Section 2.1, and the same notations are used. The iterative Jacobi algorithm for the computation of the electric fields $(\mathbf{E}_i^{n+1})_{i \in D}$ at iteration $n+1$ involves, first, the solution of the $N_{dom}$ following problems

$$ \left\{ \begin{array}{ll} \text{curl } \text{curl } \mathbf{E}_i^{n+1} - k^2 \mathbf{E}_i^{n+1} = \mathbf{0}, & \text{in } \Omega_i, \\ \gamma_i^T(\mathbf{E}_i^{n+1}) = -\gamma_i^T(\mathbf{E}^{\text{inc}}), & \text{on } \Gamma_i, \\ \gamma_i^t(\text{curl } \mathbf{E}_i^{n+1}) + \mathcal{B}(\gamma_i^T(\mathbf{E}_i^{n+1})) = \mathbf{0}, & \text{on } \Gamma_i^\infty, \\ \gamma_i^t(\text{curl } \mathbf{E}_i^{n+1}) + \mathcal{S}(\gamma_i^T(\mathbf{E}_i^{n+1})) = \mathbf{g}_{ij}^n, & \text{on } \Sigma_{ij}, \forall j \in D_i, \end{array} \right. \qquad (30) $$

and then forming the quantities $\mathbf{g}_{ji}^{n+1}$ through

$$ \mathbf{g}_{ji}^{n+1} = \gamma_i^t(\text{curl } \mathbf{E}_i^{n+1}) + \mathcal{S}(\gamma_i^T(\mathbf{E}_i^{n+1})) = -\mathbf{g}_{ij}^n + 2\mathcal{S}(\gamma_i^T(\mathbf{E}_i^{n+1})), \quad \text{on } \Sigma_{ij}, \qquad (31) $$

where, for $i \in D$, $\mathbf{E}_i = \mathbf{E}_{|\Omega_i}$, $\mathcal{S}$ is a transmission operator through the interfaces $\Sigma_{ij}$ and $\gamma_i^t$ and $\gamma_i^T$ are the local tangential trace and tangential component trace operators:

$$ \gamma_i^t : \mathbf{v}_i \mapsto \mathbf{n}_i \times \mathbf{v}_{i|\partial\Omega_i} \quad \text{and} \quad \gamma_i^T : \mathbf{v}_i \mapsto \mathbf{n}_i \times (\mathbf{v}_{i|\partial\Omega_i} \times \mathbf{n}_i), $$

with $\mathbf{n}_i$ the outward-pointing unit normal to $\Omega_i$.

Following the same procedure as in section 2.1, we introduce the two families of operators $(\mathcal{V}_i)_{i \in D}$ and $(\mathcal{T}_{ji})_{i \in D, j \in D_i}$ as:
1. $\mathbf{E}_i^{n+1} = \mathcal{V}_i(\mathbf{E}^{\text{inc}}, \mathbf{g}^n) \iff \mathbf{E}_i^{n+1}$ is solution of problem (30), where $\mathbf{g}^n = (\mathbf{g}_{ji}^n)_{i \in D, j \in D_i}$ collects all the unknowns at iteration $n$;

2. $\mathbf{g}_{ji}^{n+1} = \mathcal{T}_{ji}(\mathbf{g}_{ij}^n, \mathbf{E}_i^{n+1}) \iff \mathbf{g}_{ji}^{n+1}$ is solution of problem (31).

By linearity, we decompose the field $\mathbf{E}_i^{n+1}$ as $\mathbf{E}_i^{n+1} = \mathbf{F}_i^{n+1} + \tilde{\mathbf{E}}_i^{n+1}$, where

$$ \mathbf{F}_i^{n+1} = \mathcal{V}_i(\mathbf{E}^{\text{inc}}, 0) \quad \text{and} \quad \tilde{\mathbf{E}}_i^{n+1} = \mathcal{V}_i(0, \mathbf{g}^n). \qquad (32) $$

The quantity $\mathbf{F}_i^{n+1}$ is independent of the iteration number $n$ and can hence be written as $\mathbf{F}_i := \mathbf{F}_i^n, \forall n \in \mathbb{N}, \forall i \in D$. The whole algorithm can then be recast into a linear system:

$$ (\mathcal{I} - \mathcal{A}) \mathbf{g} = \mathbf{b}, \qquad (33) $$

that can be solved by a Krylov subspace solver.

As in the scalar case, for a vector $\mathbf{g}^n$, the quantity $\mathcal{A}\mathbf{g}^n$ is given by, for $i \in D$ and $j \in D_i$, $(\mathcal{A}\mathbf{g}^n)_{ji} = \mathcal{T}_{ji}(\mathbf{g}_{ij}^n, \tilde{\mathbf{E}}_i^{n+1})$. The information about the incident wave is contained in the right-hand side: $\mathbf{b}_{ji} = \mathcal{T}_{ji}(0, \mathbf{F}_i)$. The domain decomposition algorithm for the Maxwell system is then exactly the same as the one described in Algorithm 1 for the scalar Helmholtz equation, by formally replacing $v_i, u^{\text{inc}}, g$ and $u_i$ by $\mathbf{F}_i, \mathbf{E}^{\text{inc}}, \mathbf{g}$ and $\mathbf{E}_i$, respectively.

Similarly to the acoustic case, optimal convergence of the domain decomposition algorithm would be achieved by using the (nonlocal) MtE operator as transmission condition. Local approximations based on polynomial or rational approximations of the total symbol of the surface free-space MtE have been proposed, as well as volume representations through Perfectly Matched Layers. We detail four of those approximations below, for a generic transmission boundary $\Sigma$:

* Zeroth-order transmission condition [14]:

$$ \mathcal{S}_{\text{IBC}(0)}(\gamma^T(\mathbf{E})) = \iota k \gamma^T(\mathbf{E}). \qquad (34) $$

* Optimized second-order transmission condition [45]:

$$ \mathcal{S}_{\text{GIBC}(a,b)}(\gamma^T(\mathbf{E})) = \iota k \left( \mathcal{I} + \frac{a}{k^2} \nabla_{\Sigma} \text{div}_{\Sigma} \right)^{-1} \left( \mathcal{I} - \frac{b}{k^2} \text{curl}_{\Sigma} \text{curl}_{\Sigma} \right) (\gamma^T(\mathbf{E})), \quad (35) $$

where the curl operator is the dual operator of **curl** and where $a$ and $b$ are chosen so that an optimal convergence rate is obtained for the (TE) and (TM) modes; see [45] for the expression of $a$ and $b$ in the half-plane case. An optimized transmission condition using a single second-order operator was proposed in [1]:

$$ \mathcal{S}_{\text{GIBC}(a)}(\gamma^T(\mathbf{E})) = \iota k a \left( \mathcal{I} - \frac{1}{k^2} \text{curl}_{\Sigma} \text{curl}_{\Sigma} \right) (\gamma^T(\mathbf{E})). \qquad (36) $$

* Padé-localized square-root transmission condition [19, 21]:

$$ \mathcal{S}_{\text{GIBC}(N_p, \alpha, \varepsilon)}(\gamma^T(\mathbf{E})) = \iota k \left( C_0 + \sum_{\ell=1}^{N_p} A_{\ell} X (\mathcal{I} + B_{\ell} X)^{-1} \right)^{-1} \left( \mathcal{I} - \text{curl}_{\Sigma} \frac{1}{k_{\varepsilon}^2} \text{curl}_{\Sigma} \right) (\gamma^T(\mathbf{E})), \qquad (37) $$
with $X := \nabla_{\Sigma} \frac{1}{k_{\varepsilon}^2} \mathrm{div}_{\Sigma} - \mathbf{curl}_{\Sigma} \frac{1}{k_{\varepsilon}^2} \mathrm{curl}_{\Sigma}$, and where $k_{\varepsilon}$, $C_0$, $A_{\ell}$ and $B_{\ell}$ are defined by (13) and (14). This transmission condition corresponds to a rational approximation of the nonlocal operator

$$ \mathcal{S}_{\text{GIBC}(\text{sq}, \varepsilon)}(\gamma^T(\mathbf{E})) = \iota k (\mathcal{I} + X)^{-1/2} \left( \mathcal{I} - \mathbf{curl}_{\Sigma} \frac{1}{k_{\varepsilon}^2} \mathbf{curl}_{\Sigma} \right) (\gamma^T(\mathbf{E})), $$

which for $\varepsilon = 0$ is the principal symbol of the exact MtE operator for the half-space. As in the scalar Helmholtz case, the parameter $\varepsilon$ is introduced to regularize this operator for grazing rays on curved interfaces, and the rational approximation generalizes the polynomial approximations underlying (34), (36) and (35).

* PML transmission condition [51, 52]: The operator $\mathcal{S}_{\text{PML}(\sigma)}$ is constructed by appending a layer $\Omega^{\text{PML}}$ to the transmission interface, into which a PML transformation with absorption profile $\sigma$ is applied in the same way as for the acoustic case.

## 3.2 Convergence Analysis on a Model Problem

In order to study the convergence rate and spectral properties of the DDM algorithm we consider a similar setting as for the scalar case, but in three dimensions: the whole domain $\Omega = \mathbb{R}^3$ is separated in two curved subdomains $\Omega_1$ and $\Omega_2$ by a spherical boundary of radius $R_0$

$$ \Omega_0 := \{\mathbf{x} \in \mathbb{R}^3, |\mathbf{x}| > R_0\}, \quad \Omega_1 := \{\mathbf{x} \in \mathbb{R}^3, |\mathbf{x}| < R_0\}, \quad (38) $$

with $\partial\Omega_0 = \partial\Omega_1 := \Sigma$. Again, in this homogeneous medium setting we only consider the transmission operators that lead to a sparse surface representation. Using the same strategy as in Section 2.2, we fix $R_0 = 1$ and $k = 6\pi$, and consider a maximal number of modes $m^{\max} = [10kR]$. We report on Figure 4 the modal spectral radius $\rho(\mathcal{A}_m)$ for the transmission conditions IBC(0), GIBC(a), GIBC(a, b) and GIBC(sq, $\varepsilon$). For GIBC(a) and GIBC(a, b), the optimal parameters a and b are numerically computed by solving the min-max problem

$$ \min_{(a,b) \in \mathbb{C}^2} \max_{m \ge 1} \rho(\mathcal{A}_m) \quad (39) $$

with the Matlab function `fminsearch`. Analytical solutions of (39) for the half-space case are provided in [1] for GIBC(a) and in [45] for GIBC(a, b). Contrary the scalar Helmholtz case and to the half-space case [16], where IBC(0) leads to a convergence factor that is exactly 1 for the evanescent modes, in this model problem IBC(0) leads to $\rho(\mathcal{A}_m) < 1$ in the whole spectrum, although $\rho(\mathcal{A}_m)$ is very close to 1 for the evanescent modes, which results in a globally slowly converging DDM. For GIBC(a), we see that $\rho(\mathcal{A}_m) < 1$, for all $m$, which is improved further for GIBC(a, b). GIBC(sq, 0) leads to a better convergence rate still, which can furthermore be optimized in the transition zone by using GIBC(sq, $\varepsilon$) (with the value parameter $\varepsilon = 0.4k^{1/3}R^{-2/3}$). Finally, a numerical study using the exact series solution shows that GIBC(a) can lead to a spectral radius larger than one if the parameter a is chosen as in the half-plane case, which highlights the need for careful geometry-dependent optimization of the parameters.
<figure><img src="image_5.png" /><figcaption>Figure 4: Spectral radius of the modal iteration operator $\mathcal{A}_m$ vs. the Fourier mode $m$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Left: Residual history of GMRES vs. #iter for the various transmission conditions. Right: Eigenvalues distribution of the operator ($\mathcal{I} - \mathcal{A}$) for the different transmission conditions.</figcaption></figure>

The history of the GMRES residual with respect to the number of iterations #iter is displayed on Figure 5 (left) for the various transmission conditions. As we can observe, there is a hierarchy in the convergence curves that is directly connected to the increasing order of the GIBCs, the best convergence being obtained for $GIBC(sq, \varepsilon)$. Note that when using $GIBC(a, b)$ with the optimal parameters for the half-plane, the number of iterations is about the same as for $GIBC(a)$. Also, numerical tests show that using the Jacobi method instead of GMRES can lead to a convergence failure for $IBC(0)$, $GIBC(a)$ and $GIBC(a, b)$. The eigenvalues distribution of the operator ($\mathcal{I} - \mathcal{A}$) is displayed on Figure 5 (right). As in the scalar Helmholtz case, the improvement in the clustering of the eigenvalues around (1,0) is again observed when improving the approximation of the MtE. Finally, it shown in [21] that the localization of $GIBC(sq, \varepsilon)$ using Padé approximants behaves very similarly to the scalar Helmholtz case.
## 3.3 Weak Formulations

Without loss of generality, only the case of a particular subdomain $\Omega_i$, for $i \in D$, with no incident wave (homogeneous Dirichlet boundary condition) is detailed. We consider the same general setting as in the scalar Helmholtz case, i.e., where PML layers $\Omega_i^{PML} = \cup_{j \in D_i} \Omega_{ij}^{PML}$ are potentially appended to the artificial interfaces $\Sigma_{ij}$, and define $\Omega_i^* := \Omega_i \cup \Omega_i^{PML}$. The space of complex-valued curl-conforming vector fields on $\Omega_i^*$ is denoted by $\mathbf{H}(\text{curl}, \Omega_i^*) := \{\mathbf{W} \in (L^2(\Omega_i^*))^3 \text{ such that } \text{curl}(\mathbf{W}) \in (L^2(\Omega_i^*))^3\}$. The functional space $\mathbf{H}_0(\text{curl}, \Omega_i^*)$ is the space of functions $\mathbf{W}_i$ in $\mathbf{H}(\text{curl}, \Omega_i^*)$ such that $\gamma_i^T(\mathbf{W}_i) = 0$ on $\Gamma_i = 0$ (the boundary condition is only imposed on a part $\partial \Omega_i^*$).

* The volume PDE $\tilde{\mathbf{E}}_i^{n+1} = \mathcal{V}_i(0, \mathbf{g}^n)$ has the following weak formulation:

$$ \left\{ \begin{array}{l} \text{Find } \tilde{\mathbf{E}}_i^{n+1} \in \mathbf{H}_0(\text{curl}, \Omega_i) \text{ such that, for every } \tilde{\mathbf{E}}_i' \in \mathbf{H}_0(\text{curl}, \Omega_i): \\ \displaystyle \int_{\Omega_i} \text{curl } \tilde{\mathbf{E}}_i^{n+1} \cdot \text{curl } \tilde{\mathbf{E}}_i' \, d\Omega_i - \int_{\Omega_i} k^2 \tilde{\mathbf{E}}_i^{n+1} \cdot \tilde{\mathbf{E}}_i' \, d\Omega_i - \int_{\Gamma_i^\infty} \mathcal{B}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \tilde{\mathbf{E}}_i' \, d\Gamma_i^\infty \\ \displaystyle - \sum_{j \in D_i} \int_{\Sigma_{ij}} \mathcal{S}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij} = - \sum_{j \in D_i} \int_{\Sigma_{ij}} \mathbf{g}_{ij}^n \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij}. \end{array} \right. \qquad (40) $$

* The surface PDE $\mathbf{g}_{ji}^{n+1} = \mathcal{T}_{ji}(\mathbf{g}_{ij}^n, \tilde{\mathbf{E}}_i^{n+1})$ has the following one:

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{g}_{ji}^{n+1} \text{ in } \mathbf{H}(\text{curl}, \Sigma_{ij}) \text{ such that, for every } \mathbf{g}_{ji}' \in \mathbf{H}(\text{curl}, \Sigma_{ij}): \\ \displaystyle \int_{\Sigma_{ij}} \mathbf{g}_{ji}^{n+1} \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij} = - \int_{\Sigma_{ij}} \mathbf{g}_{ij}^n \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij} + 2 \int_{\Sigma_{ij}} \mathcal{S}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij}. \end{array} \right. $$

On the transmission boundaries, we have:

* IBC(0):

$$ \int_{\Sigma_{ij}} \mathcal{S}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} \imath k(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij}; \qquad (41) $$

$$ \int_{\Sigma_{ij}} \mathcal{S}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} \imath k(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij}. \qquad (42) $$

* GIBC(a, b):

$$ \int_{\Sigma_{ij}} \mathcal{S}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} \imath k \mathbf{r} \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij}, \qquad (43) $$

where the function $\mathbf{r} \in \mathbf{H}(\text{curl}, \Sigma_{ij})$ is obtained through the solution of

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{r} \text{ in } \mathbf{H}(\text{curl}, \Sigma_{ij}) \text{ and } \rho \text{ in } H^1(\Sigma_{ij}) \text{ such that } \forall \mathbf{r}' \in \mathbf{H}(\text{curl}, \Sigma_{ij}) \\ \text{and } \forall \rho' \in H^1(\Sigma_{ij}): \\ \displaystyle - \int_{\Sigma_{ij}} \frac{a}{k^2} \nabla_{\Sigma_{ij}} \rho \cdot \mathbf{r}' \, d\Sigma_{ij} - \int_{\Sigma_{ij}} \mathbf{r} \cdot \mathbf{r}' \, d\Sigma_{ij} + \int_{\Sigma_{ij}} \gamma_i^T(\tilde{\mathbf{E}}_i^{n+1}) \cdot \mathbf{r}' \, d\Sigma_{ij} \\ \displaystyle - \int_{\Sigma_{ij}} \frac{b}{k^2} \text{curl}_{\Sigma_{ij}}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \text{curl}_{\Sigma_{ij}} \mathbf{r}' \, d\Sigma_{ij} = 0, \\ \displaystyle \int_{\Sigma_{ij}} \rho \rho' \, d\Sigma_{ij} + \int_{\Sigma_{ij}} \mathbf{r} \cdot \nabla_{\Sigma_{ij}} \rho' \, d\Sigma_{ij} = 0; \end{array} \right. \qquad (44) $$
$$ \int_{\Sigma_{ij}} S(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} \iota k \mathbf{r} \cdot \tilde{\mathbf{g}}_{ji}' \, d\Sigma_{ij}. \quad (45) $$

• GIBC($N_p$, $\alpha$, $\varepsilon$):

$$ \int_{\Sigma_{ij}} S(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} \iota k \mathbf{r} \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij}, \quad (46) $$

where the function $\mathbf{r} \in \mathbf{H}(\text{curl}, \Sigma_{ij})$ is obtained through the solution of

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{r} \text{ in } \mathbf{H}(\text{curl}, \Sigma_{ij}), \text{ and for } \ell = 1, \dots, N_p, \varphi_\ell \text{ in } \mathbf{H}(\text{curl}, \Sigma_{ij}) \text{ and } \rho_\ell \text{ in } H^1(\Sigma_{ij}) \\ \text{such that } \forall \mathbf{r}' \in \mathbf{H}(\text{curl}, \Sigma_{ij}), \forall \varphi'_\ell \in \mathbf{H}(\text{curl}, \Sigma_{ij}) \text{ and } \forall \rho'_\ell \in H^1(\Sigma_{ij}): \\ \int_{\Sigma_{ij}} C_0 \mathbf{r} \cdot \mathbf{r}' \, d\Sigma_{ij} - \int_{\Sigma_{ij}} \gamma_i^T(\tilde{\mathbf{E}}_i^{n+1}) \cdot \mathbf{r}' \, d\Sigma_{ij} + \int_{\Sigma_{ij}} \frac{1}{k_\varepsilon^2} \text{curl}_{\Sigma_{ij}}(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \text{curl}_{\Sigma_{ij}} \mathbf{r}' \, d\Sigma_{ij} \\ + \sum_{\ell=1}^{N_p} A_\ell \left[ \int_{\Sigma_{ij}} \nabla_{\Sigma_{ij}} \rho_\ell \cdot \mathbf{r}' \, d\Sigma_{ij} - \int_{\Sigma_{ij}} \frac{1}{k_\varepsilon^2} \text{curl}_{\Sigma_{ij}} \varphi_\ell \cdot \text{curl}_{\Sigma_{ij}} \mathbf{r}' \, d\Sigma_{ij} \right] = 0, \\ \int_{\Sigma_{ij}} \varphi_\ell \cdot \varphi_\ell' \, d\Sigma_{ij} + B_\ell \left[ \int_{\Sigma_{ij}} \nabla_{\Sigma_{ij}} \rho_\ell \cdot \varphi_\ell' \, d\Sigma_{ij} - \int_{\Sigma_{ij}} \frac{1}{k_\varepsilon^2} \text{curl}_{\Sigma_{ij}} \varphi_\ell \cdot \text{curl}_{\Sigma_{ij}} \varphi_\ell' \, d\Sigma_{ij} \right] \\ - \int_{\Sigma_{ij}} \mathbf{r} \cdot \varphi_\ell' \, d\Sigma_{ij} = 0, \quad \ell = 1, \dots, N_p, \\ \int_{\Sigma_{ij}} \rho_\ell \rho_\ell' \, d\Sigma_{ij} + \int_{\Sigma_{ij}} \frac{1}{k_\varepsilon^2} \varphi_\ell \cdot \nabla_{\Sigma_{ij}} \rho_\ell' \, d\Sigma_{ij} = 0, \quad \ell = 1, \dots, N_p; \end{array} \right. \quad (47) $$

$$ \int_{\Sigma_{ij}} S(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij} := \int_{\Sigma_{ij}} \iota k \mathbf{r} \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij}. \quad (48) $$

• PML($\sigma$):

$$ \int_{\Sigma_{ij}} S(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \tilde{\mathbf{E}}_i' \, d\Sigma_{ij} := \int_{\Omega_{ij}^{\text{PML}}} D^{-1} \text{curl} \tilde{\mathbf{E}}_i^{n+1} \cdot \text{curl} \tilde{\mathbf{E}}_i' \, d\Omega_{ij}^{\text{PML}} \\ - \int_{\Omega_{ij}^{\text{PML}}} D k^2 \tilde{\mathbf{E}}_i^{n+1} \cdot \tilde{\mathbf{E}}_i' \, d\Omega_{ij}^{\text{PML}}; \quad (49) $$

$$ \int_{\Sigma_{ij}} S(\gamma_i^T(\tilde{\mathbf{E}}_i^{n+1})) \cdot \mathbf{g}_{ji}' \, d\Sigma_{ij} := \int_{\Omega_{ij}^{\text{PML}}} D^{-1} \text{curl} \tilde{\mathbf{E}}_i^{n+1} \cdot \text{curl} \mathbf{g}_{ji}' \, d\Omega_{ij}^{\text{PML}} \\ - \int_{\Omega_{ij}^{\text{PML}}} D k^2 \tilde{\mathbf{E}}_i^{n+1} \cdot \mathbf{g}_{ji}' \, d\Omega_{ij}^{\text{PML}}, \quad (50) $$

where the tensor $D$ is defined as for the acoustic case and the test functions $\mathbf{g}_{ji}'$ are again extended to the volume of the PML layers.

# 4 Numerical Implementation

The domain decomposition methods analyzed above are all readily available for testing using finite element methods in the open source GetDDM software environment [33, 48],
<figure><img src="image_7.png" /><figcaption>Figure 6: Sample models available online at http://onelab.info/wiki/GetDDM. a., b., c., d.: acoustic or electromagnetic (c. and d. only) scattering by cylindrical or spherical obstacles, with concentric or radial subdomains [10, 21]. e., f.: guided acoustic or electromagnetic waves in rectangular waveguides [51]. g.: guided acoustic or electromagnetic waves in the COBRA benchmark defined by the JINA98 workgroup [52]. h.: acoustic waves in the underground Marmousi model [47].</figcaption></figure>

available online on the web site of the ONELAB projet [29, 30]: http://onelab.info/wiki/GetDDM. GetDDM is based on the open source finite element solver GetDP (http://getdp.info) [17, 18, 28] and the open source mesh generator Gmsh (http://gmsh.info) [31, 32]. Various 2D and 3D test-cases are provided online (see Figure 6) for both acoustic and electromagnetic wave problems, as well as detailed instructions on how to build the software for parallel computer architectures. Pre-compiled, serial versions of the software for Windows, MacOS and Linux are also available for development and testing.

While GetDDM is written in C++, all the problem-specific data (geometry description, finite element formulation with appropriate transmission condition, domain decomposition algorithm) are directly written in input ASCII text files, using the code's built-in language. This general implementation allows to solve a wide variety of problems with the same software, without recompilation, and hides all the complexities of the finite element implementation from the end-user (in particular the MPI-based parallelization). Moreover, the software is designed to work both on small- and medium-scale problems (on a workstation, a laptop, a tablet or even a mobile phone) and on large-scale problems on high-performance computing clusters, without changing the input files.

One of the main features of the environment is the closeness between the input data files and the symbolic mathematical expressions of the problems. In particular, the weak formulations presented in Sections 2.3 and 3.3 are directly transcribed symbolically in the input files. For example, the relevant terms of the finite element formulation for the Maxwell problem using IBC(0) as transmission condition are directly written as follows
in the input file:

Galerkin { [ Dof{Curl E^{i}}, {Curl E^{i}} ];
         In Omega^{i}; Integration I; Jacobian V; }
Galerkin { [ -k[^2 * Dof{E^{i}}, {E^{i}} ];
         In Omega^{i}; Integration I; Jacobian V; }
Galerkin { [ -I[] * k[] * N[] /\ (Dof{E^{i}} /\ N[]), {E^{i}} ];
         In GammaInf^{i}; Integration I; Jacobian S; }
Galerkin { [ g^{i}[], {E^{i}} ];
         In Sigma^{i}; Integration I; Jacobian S; }
Galerkin { [ -I[] * k[] * N[] /\ (Dof{E^{i}} /\ N[]), {E^{i}} ];
         In Sigma^{i}; Integration I; Jacobian S; }

where $Dof{E^{i}}$ corresponds to the discrete unknown in the $i^{th}$ subdomain $Omega^{i}$ and $[.,]$ denotes the inner product. Other transmission conditions are implemented in a similar way, as is the update relation. The parallel implementation of the iterative algorithm uses the built-in function `IterativeLinearSolver`, which takes as argument the operations that implement the matrix-vector product required by Krylov subspace solvers, and is based on PETSc [5] and MUMPS [2] for the parallel (MPI-based) implementation of the linear algebra routines.

For illustration purposes, Figure 7 presents some other cases that have been solved using GetDDM. Published references are provided, which contain further information about the specific test cases, mathematical models and numerical results.

## Acknowledgements

This work was supported in part by the Wallonia-Brussels Federation (ARC grant for Concerted Research Actions ARC WAVES 15/19-03), the Belgian Science Policy (PAI grant P7/02), the Walloon Region (WIST3 grants ONELAB and ALIZEES), the French ANR (grant MicroWave NT09 460489 "Programme Blanc") and the "EADS Foundation" (High-BRID project, grant 089-1009-1006). Computational resources have been provided by CÉCI, funded by F.R.S.-FNRS (Fonds de la Recherche Scientifique) under grant n°2.5020.11, and the Tier-1 supercomputer of the Fédération Wallonie-Bruxelles, funded by the Walloon Region under grant n°1117545.

## References

[1] A. Alonso-Rodriguez and L. Gerardo-Giorda. New nonoverlapping domain decomposition methods for the harmonic Maxwell system. *SIAM J. Sci. Comput.*, 28(1):102–122, 2006.

[2] P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L'Excellent. A fully asynchronous multifrontal solver using distributed dynamic scheduling. *SIAM Journal on Matrix Analysis and Applications*, 23(1):15–41, 2001.
<figure><img src="image_8.png" /><figcaption>Figure 7: Sample models solved with GetDDM. Top: acoustic waves around a submarine (image reproduced from [10]). Bottom: electromagnetic waves around a Falcon aircraft (images reproduced from [21]).</figcaption></figure>
[3] X. Antoine, M. Darbas, and Y.Y. Lu. An improved surface radiation condition for high-frequency acoustic scattering problems. *Comput. Methods Appl. Mech. Engrg.*, 195(33-36):4060-4074, 2006.

[4] X. Antoine, C. Geuzaine, and K. Ramdani. *Wave Propagation in Periodic Media - Analysis, Numerical Techniques and Practical Applications*, volume 1, chapter Computational Methods for Multiple Scattering at High Frequency with Applications to Periodic Structures Calculations, pages 73–107. Progress in Computational Physics, 2010.

[5] S. Balay, M. F. Adams, J. Brown, P. Brune, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. Curfman McInnes, K. Rupp, B. F. Smith, and H. Zhang. PETSc Web page. http://www.mcs.anl.gov/petsc, 2015.

[6] A. Bayliss, M. Gunzburger, and E. Turkel. Boundary conditions for the numerical solution of elliptic equations in exterior regions. *SIAM J. Appl. Math.*, 42(2):430–451, 1982.

[7] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. *J. Comput. Phys.*, 114(2):185–200, 1994.

[8] A. Bermúdez, L. Hervella-Nieto, A. Prieto, and R. Rodríguez. An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems. *J. Comput. Phys.*, 223(2):469–488, 2007.

[9] Y. Boubendir. An analysis of the BEM-FEM non-overlapping domain decomposition method for a scattering problem. *J. Comput. Appl. Math.*, 204(2):282–291, 2007.

[10] Y. Boubendir, X. Antoine, and C. Geuzaine. A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation. *Journal of Computational Physics*, 231(2):262–280, 2012.

[11] Y. Boubendir, A. Bendali, and M. B. Fares. Coupling of a non-overlapping domain decomposition method for a nodal finite element method with a boundary element method. *Internat. J. Numer. Methods Engrg.*, 73(11):1624–1650, 2008.

[12] F. Collino and P. Monk. The perfectly matched layer in curvilinear coordinates. *SIAM J. Sci. Comput.*, 19(6):2061–2090 (electronic), 1998.

[13] B. Després. *Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique. Le théorème de Borg pour l'équation de Hill vectorielle*. PhD thesis, Rocquencourt, 1991. Thèse, Université de Paris IX (Dauphine), Paris, 1991.

[14] B. Després, P. Joly, and J. E. Roberts. A domain decomposition method for the harmonic Maxwell equations. In *Iterative methods in linear algebra (Brussels, 1991)*, pages 475–484, Amsterdam, 1992. North-Holland.

[15] V. Dolean, J. M. Gander, S. Lanteri, J.-F. Lee, and Z. Peng. Optimized Schwarz methods for curl-curl time-harmonic Maxwell's equations. 2013.
[16] V. Dolean, M. J. Gander, and L. Gerardo-Giorda. Optimized Schwarz methods for Maxwell's equations. *SIAM J. Sci. Comput.*, 31(3):2193–2213, 2009.

[17] P. Dular and C. Geuzaine. GetDP Web page, http://getdp.info, 2015. [online]. available: http://getdp.info.

[18] P. Dular, C. Geuzaine, F. Henrotte, and W. Legros. A general environment for the treatment of discrete problems and its application to the finite element method. *IEEE Transactions on Magnetics*, 34(5):3395–3398, September 1998.

[19] M. El Bouajaji, X Antoine, and C. Geuzaine. Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations. *Journal of Computational Physics*, 279(15):241–260, 2014.

[20] M. El Bouajaji, V. Dolean, M. Gander, and S. Lanteri. Optimized Schwarz methods for the time-harmonic Maxwell equations with damping. *SIAM Journal on Scientific Computing*, 34(4):A2048–A2071, 2012.

[21] M. El Bouajaji, B. Thierry, X. Antoine, and C. Geuzaine. A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations. *Journal of Computational Physics*, 294(1):38–57, 2015.

[22] B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. *Math. Comp.*, 31(139):629–651, 1977.

[23] B. Engquist and L. Ying. Sweeping preconditioner for the Helmholtz equation: moving perfectly matched layers. *Multiscale Model. Simul.*, 9(2):686–710, 2011.

[24] O.G. Ernst and M.J. Gander. Why it is difficult to solve Helmholtz problems with classical iterative methods. In Ivan G. Graham, Thomas Y. Hou, Omar Lakkis, and Robert Scheichl, editors, *Numerical Analysis of Multiscale Problems*, volume 83 of *Lecture Notes in Computational Science and Engineering*, pages 325–363. Springer Berlin Heidelberg, 2012.

[25] M. Gander. Optimized Schwarz methods. *SIAM Journal on Numerical Analysis*, 44(2):699–731, 2006.

[26] M. Gander and L. Halpern. Méthode de décomposition de domaine. Encyclopédie électronique pour les ingénieurs, 2012.

[27] M. J. Gander, F. Magoulès, and F. Nataf. Optimized Schwarz methods without overlap for the Helmholtz equation. *SIAM J. Sci. Comput.*, 24(1):38–60 (electronic), 2002.

[28] C. Geuzaine. GetDP: a general finite-element solver for the de Rham complex. In *PAMM Volume 7 Issue 1. Special Issue: Sixth International Congress on Industrial Applied Mathematics (ICIAM07) and GAMM Annual Meeting, Zürich 2007*, volume 7, pages 1010603–1010604. Wiley, 2008.
[29] C. Geuzaine, F. Henrotte, E. Marchandise, J.-F. Remacle, P. Dular, and R. Vazquez Sabariego. ONELAB: Open Numerical Engineering LABoratory. *Proceedings of the 7th European Conference on Numerical Methods in Electromagnetism (NUMELEC2012)*, 2012.

[30] C. Geuzaine, F. Henrotte, E. Marchandise, J.-F. Remacle, and R. Vazquez Sabariego. ONELAB Web page, http://onelab.info, 2015. [online]. available: http://onelab.info.

[31] C. Geuzaine and J.-F. Remacle. Gmsh Web page, http://gmsh.info, 2015. [online]. available: http://gmsh.info.

[32] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. *Internat. J. Numer. Methods Engrg.*, 79(11):1309–1331, 2009.

[33] C. Geuzaine, B. Thierry, N. Marsic, D. Colignon, A. Vion, S. Tournier, Y. Boubendir, M. El Bouajaji, and X. Antoine. An open source domain decomposition solver for time-harmonic electromagnetic wave problems. In *Proceedings of the 2014 IEEE Conference on Antenna and Measurements and Applications, CAMA 2014*, November 2014.

[34] D. Givoli. Computational absorbing boundaries. In Steffen Marburg and Bodo Nolte, editors, *Computational Acoustics of Noise Propagation in Fluids - Finite and Boundary Element Methods*, pages 145–166. Springer Berlin Heidelberg, 2008.

[35] R. Kerchroud, X. Antoine, and A. Soulaimani. Numerical accuracy of a Padé-type non-reflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency. *International Journal for Numerical Methods in Engineering*, 64(10):1275–1302, 2005.

[36] P.-L. Lions. On the Schwarz alternating method. III. A variant for nonoverlapping subdomains. In *Third International Symposium on Domain Decomposition Methods for Partial Differential Equations (Houston, TX, 1989)*, pages 202–223. SIAM, Philadelphia, PA, 1990.

[37] FA Milinazzo, CA Zala, and GH Brooke. Rational square-root approximations for parabolic equation algorithms. *Journal of the Acoustical Society of America*, 101(2):760–766, FEB 1997.

[38] A. Modave, E. Delhez, and C. Geuzaine. Optimizing perfectly matched layers in discrete contexts. *International Journal for Numerical Methods in Engineering*, 99(6):410–437, 2014.

[39] A. Moiola and E. A. Spence. Is the Helmholtz equation really sign-indefinite? *SIAM Rev.*, 56(2):274–312, 2014.

[40] F. Nataf. Interface connections in domain decomposition methods. *NATO Science Series II*, 75, 2001.
[41] F. Nataf and F. Nier. Convergence rate of some domain decomposition methods for overlapping and nonoverlapping subdomains. *Numer. Math.*, 75:357–377, 1997.

[42] J.-C. Nédélec. *Acoustic and electromagnetic equations*, volume 144 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 2001. Integral representations for harmonic problems.

[43] Z. Peng and J. Lee. A scalable nonoverlapping and nonconformal domain decomposition method for solving time-harmonic Maxwell equations in $\mathbb{R}^3$. *SIAM Journal on Scientific Computing*, 34(3):A1266–A1295, 2012.

[44] Z. Peng, V. Rawat, and J.-F. Lee. One way domain decomposition method with second order transmission conditions for solving electromagnetic wave problems. *Journal of Computational Physics*, 229(4):1181–1197, 2010.

[45] V. Rawat and J.-F. Lee. Nonoverlapping domain decomposition with second order transmission condition for the time-harmonic Maxwell’s equations. *SIAM J. Scientific Computing*, 32(6):3584–3603, 2010.

[46] Y. Saad and M. H. Schultz. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. *SIAM J. Sci. Statist. Comput.*, 7(3):856–869, 1986.

[47] C. Stolk. A rapidly converging domain decomposition method for the Helmholtz equation. *Journal of Computational Physics*, 241(0):240–252, 2013.

[48] B. Thierry, A. Vion, S. Tournier, M. El Bouajaji, D. Colignon, N. Marsic, X. Antoine, and C. Geuzaine. GetDDM: an open framework for testing optimized Schwarz methods for time-harmonic wave problems. *Submitted to Computer Physics Communications*, 2015.

[49] A. Toselli and O. Widlund. *Domain decomposition methods—algorithms and theory*, volume 34 of *Springer Series in Computational Mathematics*. Springer-Verlag, Berlin, 2005.

[50] A. Vion, R. Bélanger-Rioux, L. Demanet, and C. Geuzaine. A DDM double sweep preconditioner for the Helmholtz equation with matrix probing of the DtN map. In *Mathematical and Numerical Aspects of Wave Propagation WAVES 2013*, June 2013.

[51] A. Vion and C. Geuzaine. Double sweep preconditioner for optimized Schwarz methods applied to the Helmholtz problem. *Journal of Computational Physics*, 266(0):171–190, 2014.

[52] A. Vion and C. Geuzaine. Parallel double sweep preconditioner for the optimized Schwarz algorithm applied to high frequency Helmholtz and Maxwell equations. In *LNCSE, Proc. of DD22*, 2014.