# Several ways to achieve robustness when solving wave propagation problems

Niall Bootland, Victorita Dolean, Pierre Jolivet, Frédéric Nataf, S. Operto, Pierre-Henri Tournier

## ► To cite this version:

Niall Bootland, Victorita Dolean, Pierre Jolivet, Frédéric Nataf, S. Operto, et al. Several ways to achieve robustness when solving wave propagation problems. 2022. hal-03837294

HAL Id: hal-03837294

https://cnrs.hal.science/hal-03837294

Preprint submitted on 3 Nov 2022

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
# Several ways to achieve robustness when solving wave propagation problems

Niall Bootland and Victorita Dolean and Pierre Jolivet and Frédéric Nataf and Stéphane Operto and Pierre-Henri Tournier

**Abstract** Wave propagation problems are notoriously difficult to solve. Time-harmonic problems are especially challenging in mid and high frequency regimes. The main reason is the oscillatory nature of solutions, meaning that the number of degrees of freedom after discretisation increases drastically with the wave number, giving rise to large complex-valued problems to solve. Additional difficulties occur when the problem is defined in a highly heterogeneous medium, as is often the case in realistic physical applications. For time-discretised problems of Maxwell type, the main challenge remains the significant kernel in curl-conforming spaces, an issue that impacts on the design of robust preconditioners. This has already been addressed theoretically for a homogeneous medium but not yet in the presence of heterogeneities. In this review we provide a big-picture view of the main difficulties encountered when solving wave propagation problems, from the first step of their discretisation through to their parallel solution using two-level methods, by showing their limitations on a few realistic examples. We also propose a new preconditioner inspired by the idea of subspace decomposition, but based on spectral coarse spaces, for curl-conforming discretisations of Maxwell's equations in heterogeneous media.

<footnote>Niall Bootland
University of Strathclyde, Dept. of Maths and Stats, e-mail: niall.bootland@strath.ac.uk</footnote>

<footnote>Victorita Dolean
University of Strathclyde, Dept. of Maths and Stats and University Côte d'Azur, CNRS, LJAD
e-mail: work@victoritadolean.com</footnote>

<footnote>Pierre Jolivet
University of Toulouse, CNRS, IRIT, e-mail: pierre.jolivet@enseeiht.fr</footnote>

<footnote>Frédéric Nataf
Sorbonne Université, CNRS, LJLL, e-mail: frederic.nataf@sorbonne-universite.fr</footnote>

<footnote>Stéphane Operto
University Côte d'Azur, CNRS, Géoazur, e-mail: stephane.operto@geoazur.unice.fr</footnote>

<footnote>Pierre-Henri Tournier
Sorbonne Université, CNRS, LJLL e-mail: tournier@ann.jussieu.fr</footnote>
# 1 Motivation and challenges

Why do we need robust solution methods for wave propagation problems? Very often in applications, as for example in seismic inversion, we need to reconstruct the a priori unknown physical properties of an environment from given measurements. From a mathematical point of view, this means solving inverse problems by applying an optimisation algorithm to a misfit functional between the computation and the data. At each iteration of this algorithm we need to solve a boundary value problem involving the Helmholtz equation

$$-\Delta u - \frac{\omega^2}{c^2} u = f, \quad (1)$$

where $c = \sqrt{\rho c_P^2}$, $\rho$ is the density of the medium and $c_P$ is the speed of longitudinal waves. Here, $\omega$ is usually given as being the frequency of a localised source and we wish to reconstruct $c = \frac{1}{n}$ from the measured data (here, $n$ is also called the refraction index).

The Helmholtz equation is also known as the reduced wave equation or time-harmonic wave equation. Solving this equation is mathematically difficult, especially for high wave number $k = \frac{\omega}{c}$, as the solution is oscillatory and becomes more so with increasing $k$. Note that the notion of a high frequency problem is to be understood relative to the size of the computational domain: meaning how many wavelengths are present in the latter. In geophysics, the typically large size of the domain, and therefore the presence of hundreds of wavelengths, makes the problem difficult.

## 1.1 Why the time-harmonic problem in mid and high frequency is hard

What happens if one wants to approximate this problem with a numerical method? A simple computation in the one-dimensional case shows that the numerical refraction index is different from the physical one and the error depends on the product between the spacing of the grid $h$ and the frequency $\omega$, in other words numerical waves travel at a different speed to physical waves and this is also reflected in the size of error. This is also called the pollution effect and was first highlighted in the seminal paper [3]. For quasi optimality in the finite element sense we require that $h^p \omega^{p+1}$ be bounded, where $p$ is the order or the precision of the method, as shown in [21]. To summarise, the high-frequency solution $u$ oscillates at a scale $\frac{1}{\omega}$, therefore the mesh size should be chosen as at least $h \sim \frac{1}{\omega}$ leading to a large number of degrees of freedom. The pollution effect requires $h \ll \frac{1}{\omega}$, namely $h \sim \omega^{-1-\frac{1}{p}}$, therefore in practice one needs an even larger number of degrees of freedom. Note that in order to get a bounded finite element error the constraint is weaker, being $h \sim \omega^{-1-\frac{1}{2p}}$, as shown in [14]. A trade-off should be found between the number of points per wavelength (ppwl) $G = \frac{d}{h} = \frac{2\pi}{\omega h}$ and the polynomial degree $p$ in order to minimise pollution and this
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Dispersion curves for finite elements and spectral elements of order 2 and 3: normalised phase velocity (the ratio between the numerical wave speed and the physical one) as a function of the reciprocal number of points per wavelength $\frac{1}{G} = \frac{\omega h}{2\pi}$ for different discretisations. Notice that the use of higher order elements minimises dispersion even for a low value of G.</figcaption></figure>

is usually the object of dispersion analysis [1]. This is illustrated in Figure 1, where we see that the best dispersion properties are achieved when we increase the order of the discretisation or we increase G.

Suppose now that we have discretised the equation following the previous rules. We end up with a huge linear system (for a typical application we should expect millions of unknowns) whose size increases with $\omega$ very quickly, especially with more spatial dimensions. The matrix is symmetric and non-Hermitian which makes this system difficult to solve by standard iterative methods, as shown in review paper by Ernst and Gander [15] or the most recent one by Gander and Zhang [16]. Our aim should be to find the solution in optimal time for large frequencies and our algorithms should not only have good parallel properties but they should also be robust with respect to heterogeneities.

It is well-known that direct solvers, while being robust, have two main drawbacks: their high memory storage and poor parallel properties. On the other hand, iterative methods are not robust but very easy to parallelise. For this reason we consider hybrid methods, such as the naturally parallel compromise of domain decomposition methods, to obtain the best of both worlds. However, how large is truly large? In real applications, problems do not need to be over-resolved (for example, 4 ppwl are enough to perform Full Waveform Inversion with a finite-difference scheme that is specifically tuned to minimise numerical dispersion for this discretisation rule [2]) and time-harmonic Helmholtz equations with 50 million degrees of freedom were solved by a parallel direct method [20]. On the other side, when we consider much larger domains (for example via the use of a separate network of nodes rather than with cables) and that the number of nodes is limited, we must switch to iterative or hybrid methods of domain decomposition type. The methods we develop are not only motivated by the current trend in seismic imaging, meaning the development of sparse node devices (OBN) for data acquisition in the oil industry [4], but in the last decades, since the seminal work of Desprès [9], they have become the method of choice when solving the discretised Helmholtz equations.
## 2 What is the best coarse space for Helmholtz?

Consider the decomposition of the computational domain $\Omega$ into $N$ overlapping subdomains $\Omega_j$. The construction of these domains is explained later in Section 2.1 and illustrated in Figure 3. We usually solve the system $Au = b$ stemming from the finite element discretisation of (1) by a preconditioned GMRES method, e.g., in the form $M^{-1}Au = M^{-1}b$ with

$$M^{-1} = \sum_{j=1}^{N} R_j^T D_j B_j^{-1} R_j, \qquad (2)$$

where $R_j: \Omega \to \Omega_j$ is the restriction operator, $R_j^T: \Omega_j \to \Omega$ the prolongation operator and $D_j$ corresponds to the partition of unity, i.e., it is chosen such that $\sum_{j=1}^N R_j^T D_j R_j = I$. Note also that local matrices $B_j$ are stiffness matrices of local Robin boundary problems

$$(-\Delta - k^2)(u_j) = f \qquad \text{in } \Omega_j, \\ \left(\frac{\partial}{\partial n_j} + ik\right)(u_j) = 0 \qquad \text{on } \partial\Omega_j \setminus \partial\Omega.$$

We call (2) the one-level preconditioner, in particular it is the ORAS preconditioner.

Conventional wisdom in domain decomposition, backed by the definitions of strong and weak scaling, says that one-level preconditioners are not scalable (i.e., their behaviour deteriorates with the number of subdomains $N$). The crucial idea is to add a second level: that is, coarse information that is cheap to compute and immediately available to all subdomains/processors. Suppose that the coarse space is spanned by a matrix $Z$, then $E = Z^*AZ$ is the coarse matrix and $H = ZE^{-1}Z^*$ is the coarse space correction. This coarse space correction can be combined with the one-level preconditioner in an additive or hybrid manner via projectors $P$ and $Q$ ($P = Q = I$ for additive while $P = I - AH, Q = I - HA$ provides a hybrid variant)

$$M_2^{-1} = QM^{-1}P + H.$$

This coarse correction can be understood as a solution of a coarser problem on a geometrical grid with a larger spacing for example. For time-harmonic wave propagation problems, the size of the coarse grid is, however, constrained by the wave number. The theory of the grid CS (coarse space) has been introduced by Graham et al. [17] for a two-level approach to the Helmholtz problem using an equivalent problem with absorption; it has since been extended to the time-harmonic Maxwell equations [5]. This preconditioner is based on local Dirichlet boundary value problems within the one-level method. An extension to Robin transmission conditions was recently provided in [18].

The questions we would like to answer are the following: Is the grid coarse space the best choice for heterogeneous problems? Note also that the definition of the coarse space does not have to be geometrical, we can build more sophisticated coarse
spaces based on solving eigenvalue problems. Can we further improve performance by extending the idea of spectral coarse spaces to Helmholtz problems? And, if yes, what kind of modes should be included in the coarse space?

## 2.1 Spectral coarse spaces for Helmholtz

There are already now a few spectral versions of two-level preconditioners and these are DtN, H-GenEO and Δ-GenEO. For the first two there is no theory available and, while a theory has been developed for the latter, this preconditioner works mainly for low frequency and mildly non-symmetric problems.

The idea of the DtN coarse space was first introduced in [22] for elliptic problems, further analysed in [13], and extended to the Helmholtz equation in [8]. Let $D \subset \Omega$ with internal boundary $\Gamma_D = \partial D \setminus \partial \Omega$ and $v_{\Gamma_D}: \Gamma_D \to \mathbb{C}$. Then the DtN operator is defined as $\text{DtN}_D(v_{\Gamma_D}) = \frac{\partial v}{\partial n}|_{\Gamma_D}$ where $v: D \to \mathbb{C}$ is the Helmholtz extension of $v_{\Gamma_D}$ (the solution to a local boundary value problem with Dirichlet value $v_{\Gamma_D}$ on $\Gamma_D$). The DtN coarse space (introduced in [13]) is based on eigenvalue problems of the DtN operator local to each subdomain: find $(u_{\Gamma_j}, \lambda) \in V(\Gamma_j) \times \mathbb{C}$ such that

$$ \text{DtN}_{\Omega_j}(u_{\Gamma_j}) = \lambda u_{\Gamma_j}. $$

To provide the modes in the coarse space we use the Helmholtz extension $v$. We choose only eigenfunctions with $\lambda$ such that $\text{Re}(\lambda) < k_j$ where $k_j = \max_{x \in \Omega_j} k(x)$. Note that this criterion depends on the local heterogeneity in the problem and is purely heuristic (as explained in detail in [8]). In practice, finding the coarse space vectors amounts to solving local problems depending on Schur complements and mass matrices on the interfaces. By a local Helmholtz extension, we obtain vectors that, after multiplication by the partition of unity and extension by zero, form the matrix $Z$. We do this in each subdomain and combine to give the global coarse space.

The GenEO (Generalised Eigenproblems in the Overlap) coarse space was first developed in [24] for SPD problems with heterogeneous coefficients, where the heterogeneities do not align with the subdomain decomposition. More precisely, in each $\Omega_j$ we solve discrete eigenproblems with local Dirichlet matrices $A_j$ weighted by the partition of unity on one side and the local Neumann matrix $\tilde{A}_j$ on the other:

$$ D_j A_j D_j u = \lambda \tilde{A}_j u. \qquad (3) $$

We then choose only eigenfunctions with eigenvalue $\lambda$ such that $\lambda > \lambda_{\min}$. Note that if we try to replicate this exactly for Helmholtz, the method will fail. For this reason, we need to make some adaptations. The first idea is to use a nearby positive problem to build the coarse space and then use these modes for Helmholtz. This approach is called Δ-GenEO and it is amenable to theory. The second idea is more Helmholtz related in the sense that we only modify the right-hand side of the generalised eigenvalue problem (3) and thus the wave number $k$ is included in the eigenproblem:
<figure><img src="image_3.png" /><figcaption>Fig. 2 The real part of the solution to the Marmousi problem at 20 Hz. The domain is 9.2 km×3 km.</figcaption></figure>

$$D_j L_j D_j u = \lambda \tilde{A}_j u,$$

where $L_j$ corresponds to the Laplacian part of the problem and $\tilde{A}_j$ is the Neumann matrix for the Helmholtz operator. Eigenvectors associated to the eigenvalues with $Re(\lambda) > \lambda_{min}$ are now those put into the coarse space. We call this method H-GenEO.

## 2.2 Comparison of coarse spaces

We show how these three two-level methods (grid CS, DtN and H-GenEO) compare on the Marmousi<fnref n="1" /> problem [25] (see also Figure 2.2), which is a 2D geophysical benchmark problem consisting of propagation of seismic waves in a heterogeneous medium from a point source situated towards the surface. This problem is high frequency because of the large number of wavelengths in the domain. For more extensive results and comparative performance tests with these methods on other benchmark problems, see [6].

From the practical point of view, a coarse mesh is generated (from which we build the grid coarse space) and this coarse mesh is refined to give the fine mesh; see Figure 3. Alternatively, we can refine on the underlying non-overlapping decomposition and then take minimum overlap. For the discretisation by finite elements (here P2 Lagrange finite elements), we have used FreeFEM. For the domain decomposition and solver we use the FreeFEM library ffddm along with HPDDM and PETSc.<fnref n="2" />

Note that the grid CS is applied naturally to the configuration (c) (with coarse overlap), whereas for the one-level and spectral methods we can choose between minimum and coarse overlap. In the following tables we report the best results for each method in the most favourable configuration (overlap and number of modes for the spectral coarse spaces). These are the iteration counts of the GMRES method applied to the preconditioned problem with the two-level domain decomposition preconditioner in order to achieve a relative residual tolerance of $10^{-6}$. We consider

<footnote n="1">¹ https://reproducibility.org/RSF/book/data/marmousi/paper_html/node2.html</footnote>

<footnote n="2">² Software available at FreeFem-sources/examples/ffddm (ffddm) within FreeFEM, https://github.com/hpddm/hpddm (HPDDM), and https://www.mcs.anl.gov/petsc (PETSc).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Left to right: (a) the coarse mesh, (b) use of minimum overlap (refine the non-overlapping decomposition) and (c) use of coarse overlap (refine the coarse overlapping mesh directly).</figcaption></figure>

<table><caption>**Table 1** Results using the one-level and coarse grid methods for the Marmousi problem when using 5 points per wavelength, varying the frequency *f* and the number of subdomains *N*.</caption><thead><tr><th rowspan="2">f</th><th rowspan="2">\ #dofs</th><th colspan="4">One-level</th><th colspan="4">Coarse grid</th><th colspan="4">H-Geneo</th><th colspan="4">DtN</th></tr><tr><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th></tr></thead><tbody><tr><td>1</td><td>4 × 10<sup>3</sup></td><td>26</td><td>39</td><td>47</td><td>64</td><td>—</td><td>15</td><td>18</td><td>19</td><td>20</td><td>—</td><td>9</td><td>11</td><td>17</td><td>21</td><td>—</td><td>6</td><td>7</td><td>9</td><td>6</td><td>—</td></tr><tr><td>5</td><td>1 × 10<sup>5</sup></td><td>53</td><td>76</td><td>105</td><td>154</td><td>213</td><td>26</td><td>29</td><td>28</td><td>29</td><td>31</td><td>15</td><td>17</td><td>26</td><td>37</td><td>56</td><td>7</td><td>19</td><td>10</td><td>8</td><td>19</td></tr><tr><td>10</td><td>5 × 10<sup>5</sup></td><td>68</td><td>102</td><td>158</td><td>212</td><td>302</td><td>32</td><td>35</td><td>41</td><td>40</td><td>42</td><td>33</td><td>40</td><td>45</td><td>56</td><td>73</td><td>18</td><td>19</td><td>21</td><td>48</td><td>29</td></tr><tr><td>20</td><td>2 × 10<sup>6</sup></td><td>82</td><td>125</td><td>178</td><td>248</td><td>347</td><td>34</td><td>35</td><td>42</td><td>43</td><td>44</td><td>64</td><td>83</td><td>121</td><td>134</td><td>157</td><td>43</td><td>75</td><td>77</td><td>61</td><td>35</td></tr></tbody></table>

<table><caption>**Table 2** Results using the one-level and coarse grid methods for the Marmousi problem when using 10 points per wavelength, varying the frequency *f* and the number of subdomains *N*.</caption><thead><tr><th rowspan="2">f</th><th rowspan="2">\ #dofs</th><th colspan="4">One-level</th><th colspan="4">Coarse grid</th><th colspan="4">H-Geneo</th><th colspan="4">DtN</th></tr><tr><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th><th>10</th><th>20</th><th>40</th><th>80</th><th>160</th></tr></thead><tbody><tr><td>1</td><td>2 × 10<sup>4</sup></td><td>30</td><td>43</td><td>63</td><td>97</td><td>—</td><td>16</td><td>18</td><td>19</td><td>21</td><td>—</td><td>7</td><td>8</td><td>8</td><td>13</td><td>—</td><td>4</td><td>7</td><td>5</td><td>6</td><td>—</td></tr><tr><td>5</td><td>5 × 10<sup>5</sup></td><td>58</td><td>87</td><td>126</td><td>175</td><td>246</td><td>29</td><td>29</td><td>34</td><td>34</td><td>36</td><td>10</td><td>9</td><td>10</td><td>10</td><td>12</td><td>10</td><td>11</td><td>12</td><td>17</td><td>24</td></tr><tr><td>10</td><td>2 × 10<sup>6</sup></td><td>78</td><td>124</td><td>172</td><td>251</td><td>346</td><td>35</td><td>41</td><td>43</td><td>46</td><td>45</td><td>20</td><td>16</td><td>14</td><td>13</td><td>13</td><td>19</td><td>23</td><td>25</td><td>25</td><td>24</td></tr><tr><td>20</td><td>8 × 10<sup>6</sup></td><td>92</td><td>142</td><td>198</td><td>272</td><td>389</td><td>39</td><td>47</td><td>48</td><td>49</td><td>49</td><td>45</td><td>40</td><td>34</td><td>25</td><td>19</td><td>35</td><td>46</td><td>48</td><td>56</td><td>59</td></tr></tbody></table>

two scenarios: the *under-resolved* case with a discretisation of 5 points per wavelength (Table 1) and the *over-resolved* case with a discretisation of 10 points per wavelength (Table 2) and vary the frequency and the number of subdomains. Low resolution is motivated by applications where high precision is not needed, especially when solving inverse problems by FWI (Full Waveform Inversion). In this case, since the test cases are large, one needs to find a good trade-off between precision and the size of the system to be solved. We refer the reader to the references [11, 12] where a more extensive numerical study was performed. We notice that in the first scenario the grid CS outperforms the spectral methods (with a slight advantage over the DtN method) whereas in the second scenario the H-GenEO method displays the best performance.

We conclude this comparison by noting that there is no clear advantage in one method over another, all depends on the frequency and precision desired. We have not sought an optimal implementation and the grid CS is the finest possible (which is in principle very expensive), in this sense the timings are not relevant, even if the
cost per iteration might be different. In the case of multiple right-hand sides, spectral coarse spaces may have an advantage, although we have not studied this aspect here.

For large-scale geophysical example problems, we have explored extensively the performance of the grid coarse space in [10, 11]. A few conclusions are stated below:

* The use of higher order finite elements allow us to minimise dispersion with a minimum number of ppwl, as shown in Figure 1. A good compromise is the choice of P3 finite elements for which, with 5 ppwl on unstructured meshes, we note a reduction by a factor 2 in the number of degrees of freedom with respect to a finite difference discretisation on uniform meshes.

* Local solves in domain decomposition methods are usually done by direct methods such as Cholesky factorisation, which is part of the setup phase ahead of the application of the GMRES method. We can already improve performance by replacing the Cholesky method with incomplete Cholesky factorisation.

* Precision is also important in the parsimony of the computation and the use of single precision highly decreases both the setup and solution times.

### 3 Can we improve on the auxiliary subspace preconditioner?

Let us consider the positive (or time-discretised) Maxwell equations

$$ 
\begin{array}{ll}
\nabla \times (\mu_r^{-1} \nabla \times \mathbf{u}) + \alpha \varepsilon_r \mathbf{u} = \mathbf{f} & \text{in } \Omega, \\
\mathbf{u} \times \mathbf{n} = 0 & \text{on } \partial \Omega.
\end{array}
$$

Here $\mathbf{u}$ is the vector-valued electric field, $\mathbf{f}$ is the source term, $\alpha > 0$ is a constant (e.g., stemming from the time discretisation), and $\mu_r$ and $\varepsilon_r$ are electromagnetic parameters which are uniformly bounded and strictly positive but which we allow to be heterogeneous. We suppose $\Omega$ is a polyhedral computational domain and $\mathbf{n}$ is the outward normal to $\partial \Omega$. After discretisation by Nédélec elements we obtain

$$ A\mathbf{U} := (K + \alpha M)\mathbf{U} = \mathbf{b}, \quad (4) $$

where $K \in \mathbb{R}^{n \times n}$ represents the discretisation of the curl-curl operator $\nabla \times (\mu_r^{-1} \nabla \times)$ and $M \in \mathbb{R}^{n \times n}$ is the $\varepsilon_r$-weighted mass matrix computed in the edge element space. Note that matrix $K$ has a huge kernel (all the gradients of $H^1$ functions are part of the kernel of the curl operator) so designing efficient preconditioners for this problem can be challenging.

There is a well-established preconditioner in the literature known as the (nodal) auxiliary space preconditioner (ASP) [19] which is based on a splitting of the space, here $\mathbf{H}(\text{curl})$, by isolating the kernel. The auxiliary space then uses a nodal (Lagrangian) discretisation. The preconditioner is given as

$$ M_{ASP}^{-1} = \text{diag}(A)^{-1} + P(\tilde{L} + \alpha \tilde{Q})^{-1} P^T + \alpha^{-1} C L^{-1} C^T, $$
where $\tilde{L}$ is the nodal discretisation of the $\mu_r^{-1}$-weighted vector Laplacian operator, $\tilde{Q}$ is the nodal $\varepsilon_r$-weighted vector mass matrix, $P$ is the matrix form of the nodal interpolation operator between the Nédélec space and nodal element space, and $C$ is the “gradient matrix”, which is exactly the null space matrix of $A$ here.

The spectral condition number $\kappa_2(M_{ASP}^{-1}A)$ of the preconditioned problem is independent of the mesh size but might depend on any heterogeneities present. The natural question is then whether we can improve upon this preconditioner in the case of heterogeneous Maxwell problems?

In order to do this, we use extensively the fictitious space lemma (FSL) of Nepomnyaschickh, which can be considered the Lax–Milgram theorem of domain decomposition [23].

**Lemma 1 (Nepomnyaschickh, 1991)** Consider two Hilbert spaces $H$ and $H_D$ along with positive symmetric bilinear forms $a: H \times H \to \mathbb{R}$ and $b: H_D \times H_D \to \mathbb{R}$. The operators $A$ and $B$ are defined as follows

* $A: H \to H$ such that $(Au, v) = a(u, v)$ for all $u, v \in H$;

* $B: H_D \to H_D$ such that $(Bu_D, v_D)_D = b(u_D, v_D)$ for all $u_D, v_D \in H_D$.

Suppose we have a linear surjective operator $\mathcal{R}: H_D \to H$ verifying the properties

* Continuity: $\exists c_R > 0$ such that $\forall u_D \in H_D$ we have

$$a(\mathcal{R}u_D, \mathcal{R}u_D) \le c_R b(u_D, u_D).$$

* Stable decomposition: $\exists c_T > 0$ such that $\forall u \in H \exists u_D \in H_D$ with $\mathcal{R}u_D = u$ and

$$c_T b(u_D, u_D) \le a(\mathcal{R}u_D, \mathcal{R}u_D) = a(u, u).$$

Consider the adjoint operator $\mathcal{R}^*$: $H \to H_D$ given by $(\mathcal{R}u_D, u) = (u_D, \mathcal{R}^*u)_D$ for all $u_D \in H_D$ and $u \in H$. Then for all $u \in H$ we have the spectral estimate

$$c_T a(u, u) \le a(\mathcal{R}B^{-1}\mathcal{R}^*Au, u) \le c_R a(u, u).$$

Thus, the eigenvalues of the preconditioned operator $\mathcal{R}B^{-1}\mathcal{R}^*A$ are bounded from below by $c_T$ and from above by $c_R$.

In this lemma we have a few ingredients: two Hilbert spaces with the associated scalar products (that are linked by the surjective operator $\mathcal{R}$) and two symmetric positive bilinear forms. The first of each comes from our problem while the second is for the preconditioner. Under the assumptions of continuity and stable decomposition, the spectral estimate tells us that the spectral condition number of the preconditioned problem is bounded solely in terms of the constants $c_R$ and $c_T$.

Discretised problems which are perturbations of a singular operator, such as the Maxwell problem in (4) when $\alpha$ is small, have a huge near-kernel $G \subset \mathbb{R}^n$ of $A$, given by the gradient of all $H^1(\Omega)$ functions for example. This near-kernel will be within a space $V_G \subset \mathbb{R}^n$, which is the vector space spanned by the sequence $(R_i^T D_i R_i G)_{1 \le i \le N}$ so that $G \subset V_G$. These spaces may not be equal due to the fact
that not all the elements of $R_i^T D_i R_i G$ are in $G$, for example, corresponding to the degrees of freedom for which $D_i$ is not locally constant. Nevertheless, since the $D_i$ are related to a partition of unity, we guarantee the inclusion. The space $V_G$ can now serve as a “free” coarse space. We denote the coarse space $V_0 := V_G$ and let $Z \in \mathbb{R}^{n_0 \times n}$ be a rectangular matrix whose columns are a basis of $V_0$. The coarse space matrix is then defined in the usual way by $E = Z^T AZ$.

We now need to define all the other ingredients in the FSL. The second Hilbert space is the product space of vectors stemming, for example, from the $n_i$ degrees of freedom on the local subdomains $\Omega_i$ and the $n_0$ coarse space vectors

$$H_D := \mathbb{R}^{n_0} \times \prod_{i=1}^{N} \mathbb{R}^{n_i}.$$

The bilinear form $b$ for the preconditioner is given by the sum of local bilinear forms $b_i$ and the coarse space contribution

$$b(\mathcal{U}, \mathcal{V}) := (E\mathbf{U}_0, \mathbf{V}_0) + \sum_{i=1}^{N} b_i(\mathbf{U}_i, \mathbf{V}_i), \quad b_i(\mathbf{U}_i, \mathbf{V}_i) := (R_i A R_i^T \mathbf{U}_i, \mathbf{V}_i),$$

for $\mathcal{U} = (\mathbf{U}_0, (\mathbf{U}_i)_{1 \le i \le N}) \in H_D$, $\mathcal{V} = (\mathbf{V}_0, (\mathbf{V}_i)_{1 \le i \le N}) \in H_D$. Finally, the surjective operator $\mathcal{R}_{AS}: H_D \to H$ corresponding the additive Schwarz method is given by

$$\mathcal{R}_{AS}(\mathcal{U}) := Z\mathbf{U}_0 + (I - P_0) \sum_{i=1}^{N} R_i^T \mathbf{U}_i,$$

where $P_0$ is the A-orthogonal projection on the coarse space $V_0$. By applying the FSL we obtain a spectral condition number estimate $\kappa(M_{AS}^{-1}A) \le C$, with a bound $C$ that can be large due to heterogeneities in the problem.

How can we improve this preconditioner in this case in order to be robust? We simply build a GenEO space from local generalised eigenproblems in the orthogonal complement of the “free” coarse space: find $(\mathbf{V}_{jk}, \lambda_{jk}) \in \mathbb{R}^{n_j} \setminus \{0\} \times \mathbb{R}$ such that

$$(I - \xi_{0j}^T) D_j R_j A R_j^T D_j (I - \xi_{0j}) \mathbf{V}_{jk} = \lambda_{jk} \tilde{A}_j \mathbf{V}_{jk},$$

where $\xi_{0j}$ denotes the $b_j$-orthogonal projection from $\mathbb{R}^{n_j}$ on $G_j = R_j G$ and $\tilde{A}_j$ is the local Neumann matrix for the problem. We define $V_{j,geneo}^T \subset \mathbb{R}^n$ to be the vector space spanned by the family of vectors $(R_j^T D_j (I - \xi_{0j}) \mathbf{V}_{jk})_{\lambda_{jk}>\tau}$ corresponding to eigenvalues larger than a chosen threshold parameter $\tau$. Now, collecting over all subdomains $j$, we let $V_{geneo}^T$ be the span of all $(V_{j,geneo}^T)_{1 \le j \le N}$, which will lead to a new coarse space

$$V_0 := V_G + V_{geneo}^T.$$

Applying the FSL now yields a spectral condition number estimate of the resulting two-level Schwarz method which is independent of the heterogeneity in the problem.
Several other variants of this approach can be formulated, including with the use of inexact coarse solves in order to more efficiently handle the large coarse space; these theoretical advances can be found in our recent preprint [7].

# 4 General conclusions

In this short paper we have offered a brief overview of the main difficulties and some recent solution methods now available to solve Helmholtz equations in the mid and high frequency regimes, which occur in many applications and especially in geophysics. Although there is no established method as the go-to solver, we have proposed a number of different strategies based on two-level domain decomposition methods where the second level comes from the solution of local spectral problems. Indeed, spectral coarse spaces have shown excellent theoretically-proven results for symmetric positive definite problems and currently offer very promising directions to explore for the Helmholtz equation and other wave propagation problems.

The discretisation here is also intertwined with the solution method as solvability and accuracy are very important for wave propagation problems. However, problems in applications do not need to be over-resolved (for example, in full waveform inversion for a discretisation by a finite difference method minimising dispersion 4 ppwl are enough) as this can lead to increasingly large problems whose size is not fully justified by practical reasons. Further, multi-frontal direct solvers based on block low rank approximations have been developed in recent years and problems as large as 50 million unknowns can be tackled successfully by these methods. In this sense, domain decomposition solvers need to be designed with the idea to go beyond these limits while keeping the applicative context in mind.

Last but not least, while not of the same nature, positive Maxwell's equations present different challenges. Here, the auxiliary space preconditioner has successfully been applied to problems where the underlying operator has an infinite dimensional kernel. By exploiting the idea of subspace decomposition together with spectral methods of GenEO type, a new generation of preconditioners, capable of tackling heterogeneous problems, has been introduced. Future work includes an extensive numerical exploration of such an approach on realistic example problems.

Wave propagation problems have been a key source of difficult problems not just for domain decomposition but more widely in scientific computing. As large-scale computing infrastructure continues to evolve and practitioners become ever more ambitious, often driven by industrial challenges, robustness will remain a central theme when designing algorithms for the future. Our work here then contributes some of the most recent ideas towards achieving such desired robustness for domain decomposition methods applied to challenging applications in wave propagation.
**Acknowledgements** The first two authors gratefully acknowledge support from the EPSRC grant EP/S004017/1. The fifth author acknowledges support from the Wind project<fnref n="3" /> funded by Shell, Total and Chevron.

References

1. M. Ainsworth, H. A. Wajid: Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration, SIAM J. Numer. Anal., 48(1), pp. 346–371 (2010).

2. P. Amestoy, R. Brossier, A. Buttari, J.-Y. L'Excellent, T. Mary, L. Métivier, A. Miniussi, S. Operto: Fast 3D frequency-domain full-waveform inversion with a parallel block low-rank multifrontal direct solver: Application to OBC data from the North Sea, Geophysics, 81(6), pp. R363–R383 (2016).

3. I. M. Babuška, S. A. Sauter: Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?, SIAM J. Numer. Anal., 34(6), pp. 2392–2423 (1997).

4. J. Blanch, J. Jarvis, C. Hurren, Y. Liu, and L. Hu: Designing an exploration scale OBN: Acquisition design for subsalt imaging and velocity determination. In *SEG Technical Program Expanded Abstracts 2019*, pp. 192–196. Society of Exploration Geophysicists (2019).

5. M. Bonazzoli, V. Dolean, I. G. Graham, E. A. Spence, P.-H. Tournier: Domain decomposition preconditioning for the high-frequency time-harmonic Maxwell equations with absorption, Math. Comp. 86, pp. 2089–2127 (2017).

6. N. Bootland, V. Dolean, P. Jolivet, P.-H. Tournier: A comparison of coarse spaces for Helmholtz problems in the high frequency regime, arXiv preprint arXiv:2012.02678 (2020).

7. N. Bootland, V. Dolean, F. Nataf, P.-H. Tournier: Two-level DDM preconditioners for positive Maxwell equations, arXiv preprint arXiv:2012.02388 (2020).

8. L. Conen, V. Dolean, R. Krause, F. Nataf: A coarse space for heterogeneous Helmholtz problems based on the Dirichlet-to-Neumann operator, J. Comput. Appl. Math., 271, pp. 83–99 (2014).

9. B. Després: Domain decomposition method for the Helmholtz problem, C. R. Math. Acad. Sci. Paris. I Math., 311(6), pp. 313–316 (1990).

10. V. Dolean, P. Jolivet, P.-H. Tournier, S. Operto: Iterative frequency-domain seismic wave solvers based on multi-level domain-decomposition preconditioners, 82nd EAGE Annual Conference & Exhibition 2020 (1), pp. 1–5 (2020).

11. V. Dolean, P. Jolivet, P.-H. Tournier, S. Operto: Large-scale frequency-domain seismic wave modeling on h-adaptive tetrahedral meshes with iterative solver and multi-level domain-decomposition preconditioners, SEG Technical Program Expanded Abstracts 2020, pp. 2683–2688 (2020).

12. V. Dolean, P. Jolivet, P.-H. Tournier, L. Combe, S. Operto, S. Riffo:Large-scale finite-difference and finite-element frequency-domain seismic wave modelling with multi-level domain-decomposition preconditioner, arXiv:2103.14921, (2021).

13. V. Dolean, F. Nataf, R. Scheichl, N. Spillane: Analysis of a two-level Schwarz method with coarse spaces based on local Dirichlet-to-Neumann maps, Comput. Methods Appl. Math. 12(4), pp. 391–414 (2012).

14. Y. Du, H. Wu: Preasymptotic error analysis of higher order FEM and CIP-FEM for Helmholtz equation with high wave number, SIAM J. Numer. Anal., 53(2), pp. 782–804 (2015).

15. O. G. Ernst., M. J. Gander: Why it is difficult to solve Helmholtz problems with classical iterative methods. In: Graham I., Hou T., Lakkis O., Scheichl R. (eds) Numerical Analysis of Multiscale Problems. LNCSE, vol 83. Springer, Berlin, Heidelberg (2012).

<footnote n="3">³ https://www.geoazur.fr/WIND/</footnote>
16. M. J. Gander, H. Zhang: A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods, SIAM Rev., 61(1), pp. 3–76 (2019).

17. I. G. Graham, E. A. Spence, E. Vainikko: Domain decomposition preconditioning for high-frequency Helmholtz problems with absorption, Math. Comp. 86, pp. 2089–2127 (2017).

18. I. G. Graham, E. A. Spence, J. Zou: Domain Decomposition with local impedance conditions for the Helmholtz equation with absorption, SIAM J. Numer. Anal., 58(5), pp. 2515–2543 (2020).

19. R. Hiptmair, J. Xu: Nodal auxiliary space preconditioning in $H(\text{curl})$ and $H(\text{div})$ spaces, SIAM J. Numer. Anal., 45(6), pp. 2483–2509 (2007).

20. T. Mary: Block low-rank multifrontal solvers: complexity, performance and scalability. PhD Thesis, Université de Toulouse (2017).

21. J. M. Melenk, S. A. Sauter: Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtz equation, SIAM J. Numer. Anal., 49(3), pp. 1210–1243 (2011).

22. F. Nataf, H. Xiang, V. Dolean, N. Spillane: A coarse space construction based on local Dirichlet-to-Neumann maps, SIAM J. Sci. Comput., 33(4), pp. 1623–1642 (2011).

23. S. V. Nepomnyaschikh. Mesh theorems of traces, normalizations of function traces and their inversions. Sov. J. Numer. Anal. Math. Modeling, 6(3), pp. 223–242 (1991).

24. N. Spillane, V. Dolean, P. Hauret, F. Nataf, C. Pechstein, R. Scheichl: Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps, Numer. Math. 126(4), pp. 741–770 (2014).

25. R. Versteeg: The Marmousi experience: Velocity model determination on a synthetic complex data set, The Leading Edge 13(9), pp. 927–936, doi:10.1190/1.1437051 (1994).