Sweeping preconditioners for stratified media in the presence of reflections

In this paper we consider sweeping preconditioners for stratified media, especially in the presence of reflections. In the most famous class of sweeping preconditioners Dirichlet-to-Neumann operators for half-space problems are approximated through absorbing boundary conditions. In the presence of reflections absorbing boundary conditions are not accurate resulting in an unsatisfactory performance of these sweeping preconditioners. We explore the potential of using more accurate Dirichlet-to-Neumann operators within the sweep. To this end, we make use of the separability of the equation for the background model. While this improves the accuracy of the Dirichlet-to-Neumann operator, we find both from numerical tests and analytical arguments that it is very sensitive to perturbations in the presence of reflections. This implies that even if accurate approximations to Dirichlet-to-Neumann operators can be devised for a stratified medium, sweeping preconditioners are limited to very small perturbations.


Introduction
Time harmonic wave equations arise in various applications. Their numerical discretization leads to large linear systems which are difficult to solve with classical iterative methods [4]. Recently, sweeping preconditioners have emerged as a promising approach to overcome this problem [5]. Since the introduction of the moving perfectly matched layer (PML) preconditioner by Engquist and Ying [3], numerous impressive results and further developments of this technique have been published. We refer to [5] for a comprehensive review.
Unfortunately, the range of problems in which sweeping preconditioners can be used is limited. We are not aware of any publication in which sweeping preconditioners have been successfully applied to media that contain strong resonant cavities. In fact, numerical experiments indicate that the established sweeping methods are not suitable for treating such problems, see e.g. section 7.4 of [11] or section 10 of [5]. Additionally, sweeping preconditioners require an absorbing boundary condition on at least one boundary of the domain at which the process of sweeping can be started. Since these assumptions are violated in many practically relevant problems, e.g. from global seismology, it is important to explore if these limitations of the sweeping technique can be overcome.
This paper investigates this question for the case of stratified media. As a concrete example we consider a problem from [6,7] in which spherical coordinates {r, ϕ, θ} are used. Assuming axisymmetry of all fields in ϕ the propagation of shear waves (SH-waves) u between the core mantle boundary (CMB) at r = R CMB and the surface of the Earth at r = R ⊕ in the frequency domain is described by the equation Lu = f, for (r, θ) ∈ Ω := (R CMB , R ⊕ ) × (0, π), where Lu := −ρω 2 ur 2 sin 2 (θ) − sin 2 (θ) r 2 ∂ ∂r r 4 µ ∂u ∂r − 1 sin(θ) ∂ ∂θ sin 3 (θ)µ ∂u ∂θ (2) and f = f s r 2 sin 2 (θ), subject to boundary conditions Bu = 0. The boundary operator is defined piecewise on ∂Ω = Γ D ∪ Γ N with Γ D = {θ = 0} ∪ {θ = π} and Γ N = {r = R CMB } ∪ {r = R ⊕ } by Here, ρ is the mass density, µ = ρv 2 SH is the shear modulus and ω the frequency. The background coefficients for v SH and ρ are provided by the spherically symmetric PREM model [2]. In the appendix, cf. Section A, we give a derivation of the variational formulation that we base the finite element discretization on. Conventional sweeping preconditioners cannot be applied to this problem since absorbing boundary conditions are missing. In this paper an extension of the sweeping preconditioner is presented which overcomes this limitation for the spherically symmetric background model. It will further be investigated to which extent this preconditioner can then also be used for other models in which the coefficients ρ and v SH are small perturbations from the spherically symmetric case. This would be realistic, since 3D tomographic models of the Earth deviate only a few percent from the background model [1].
The remainder of this paper is structured as follows. In Section 2 we recall the general framework of sweeping preconditioners and the commonly employed moving PML approximation of the Dirichlet-to-Neumann (DtN) operator. As a prototype of an improved approximation of the DtN operator we construct a DtN operator based on the separability of the background problem on the discrete level in Section 3. The potential of this method is explored in Section 4 with numerical experiments. Here it is observed that the DtN for the free surface boundary condition is very sensitive to perturbations. This observation is discussed in detail in Section 5. We draw some conclusions regarding the applicability of sweeping preconditioners in the presence of reflections in Section 6.

General framework of sweeping preconditioners
In [5] several sweeping methods have been described in the framework of the double sweep optimized Schwarz method (DOSM). Here, we adapt DOSM to our specific setting. This includes the restriction to special cases, e.g. only non-overlapping domain decompositions and DtN transmission conditions are considered. The simplified version is more appropriate for the purpose of this paper since it allows to focus attention on the approximation of the DtN map.

Double sweep optimized Schwarz method
Let Ω = J j=1 Ω j be a non-overlapping decomposition of the domain into horizontal layers Ω j , see Figure 1. In the numerical examples below, we make the choice that one layer Ω j contains two finite elements in vertical direction (r) and 2 · J elements in the horizontal direction (θ) yielding a fixed aspect ratio for all finite elements. We denote by Γ j,j±1 := ∂Ω j ∩ ∂Ω j±1 the interfaces between the layers and write u j for a function defined in layer Ω j . Based on these definitions we state the specialized sweeping algorithm.
Forward sweep: Given the last iterate u (n−1) j in Ω j , j = 1, . . . , J solve successively for j = 1, . . . , J − 1 Backward sweep: Solve successively for j = J, . . . , 1 The transmission operator P j is an (approximation) of the DtN map where v solves with Ω ext If DtN j is well defined for j = 2, . . . , J and the original problem and subdomain problems are uniquely solvable then DOSM converges in one double sweep to the exact solution for P j = DtN j [5]. In practice, using the exact DtN map as a transmission operator is computationally too expensive. Therefore, P j is chosen as an approximation of the DtN. The algorithm is then usually used as a preconditioner for GMRES with initial guess u (0) j = 0, j = 1, . . . , J.

Moving PML approximation of the DtN
The problem for calculating the exact DtN j is posed on the whole exterior domain Ω ext j . Usually, it is assumed that on ∂Ω ∩ ∂Ω 1 an absorbing boundary condition implemented by a PML is present. As the name "moving PML" suggests this PML is shifted closer to Ω j . This replaces the original problem posed on j−1 i=1 Ω i by a modified problem on the (usually smaller) domain Ω PML j . In practice, the PML is usually started right at the coupling interface Γ j,j−1 and Ω PML j = Ω j−1 . This leads to the operator where Here, the original differential operators L, B have been replaced by modified versionsL,B due to the complex scaling applied in the PML region.

A tensor product approximation of the DtN
In this section a new approach to approximate the DtN for tensor product discretizations will be presented. Let us note that the purpose of this approach is primarily to explore the potential of sweeping preconditioners under the assumption of accurate DtN approximations in the sequel. Let us fix an interface Γ j,j−1 at r = R j on which the DtN j shall be com- . This problem will be solved in two steps. In section 3.1 we treat the case in which g h is a discrete eigenfunction of the weighted Laplacian on Γ j,j−1 . The action of the DtN on such data is given by multiplication with a number which can be computed by solving an ODE. In other words, the DtN is diagonal in the basis of the discrete eigenfunctions. This allows to treat the case of general Dirichlet data in section 3.2 by a simple change of basis.

The DtN applied to a discrete eigenfunction
for all ϑ h ∈ V θ h . The next proposition shows that the ψ are also eigenfunctions of the discretized DtN map.

Remark 1 The finite element space V
Rj h for the ODE (10) has to coincide with the first factor of the tensor product W In other words, the ODE discretization is chosen as the restriction of the two dimensional discretization to the radial direction.

General Dirichlet data
To apply the DtN to general Dirichlet data g h ∈ V θ h we express it in the eigenbasis g h = L =1 g ψ and apply the DtN as Then we transform back to the finite element basis. The transformation to the discrete eigenbasis involves the solution of a dense linear system which is composed of the eigenvectors in the finite element basis. In case of a uniform 2D mesh with periodic boundary conditions, the boundary mass and stiffness matrices are discrete block circulant matrices, and the transformation is given by the tensor (Kronecker) product of a small matrix and an FFT matrix, cf. [8]. Possibly the inversion of dense matrices can also be avoided for more general grids in some cases by designing analogues of infinite elements adapted to this setting, but the results of this paper may discourage from starting this effort.

Numerical experiments
In this section numerical experiments for the model problem will be presented. First, the shortcomings of the moving PML approximation of the DtN are demonstrated. Then we apply our new approximation to the spherically symmetric model and investigate its performance in case of small perturbations. All experiments are carried out using 4th order H 1 -conforming finite elements and have been implemented in the finite element library Netgen/NGSolve, see [9,10]. Scripts for reproducing the numerical results are provided at DOI: 10.5281/zenodo.3686802.

Moving PML
The moving PML preconditioner relies on the assumption that the computational domain is truncated by an absorbing boundary condition in at least one direction and that the medium is free of large resonant cavities. The following two examples demonstrate that these assumptions are crucial. Periodic boundary conditions in θ are used. At r = 1.0 an absorbing boundary condition implemented by a PML is set, while at r = 0 we simply use natural boundary conditions. The geometrical setup is as shown in Figure 1. Similar to the experiments from [5] we let the wavespeed vary discontinuously between the layers where the strength of the discontinuity is given by a factor α. To this end, we set on the first layer c = 1/(1 + α), on the next c = 1/(1 − α) then again c = 1/(1+α) continuing in this fashion. The density is simply ρ = 1.
We apply the DOSM with moving PML approximation of the DtN to this problem. The results are shown in Table 1. We also consider the case in which additional damping is added to the preconditioner 1 , i.e. the operatorL in equation (7) describes a Helmholtz problem with complex frequency ω + iγ, where γ = 1. Adding additional damping to the preconditioner leads to nearly robust iteration numbers for α = 0. However, for both versions the iteration numbers increase drastically as the contrast α is increased. As a result, the preconditioner becomes completely inefficient.

SH-waves in frequency domain
In the academic example reflections generated by the strongly discontinuous coefficients lead to the method's breakdown. Severe reflections can also be generated by boundary conditions as is the case for the model problem (1) from seismology. This can be demonstrated by comparing the GMRES iteration numbers for two different boundary conditions: -In the first case, a PML is implemented at the Earth's surface.
-In the second case, the PML is removed, which realizes the physically desired free surface condition.
The tolerance is set to 10 −7 . Two sources are considered: a Dirac and a random source.
The iteration numbers for the case of a PML at the Earth's surface are shown in Table 2 in the columns 'PML' . The wave field for ω = 2048 with the Dirac source is shown in Figure 2 (a). The waves generated from the point source bounce off from the CMB. Additional reflections are caused by the discontinuous coefficients and the Dirichlet boundary conditions at θ = 0 and θ = π. Despite these difficulties, the preconditioner performs well. The iteration numbers grow only very mildly and the low tolerance of 10 −7 is achieved for ω = 2048 in 15 iterations for both sources.
The situation changes drastically when the PML at the Earth's surface is removed. The solver now has to capture a very complex wavefield created by additional reflections from the Earth's surface as Figure 2 (b) shows. The results for this case are given in Table 2 in the columns 'free'. The iteration numbers now grow drastically with increasing wavenumber. In Figure 3 it can be seen that the residual stagnates for many iterations. GMRES first has to filter out the reflections until a convergence to the solution occurs. This renders the moving PML preconditioner unsuitable for our purpose. Let us note that in additional computational studies we also tried adding damping to the preconditioning problem which however did not improve the situation.

Tensor product DtN
We tackle problem (1) with the free surface boundary condition using a tensor product discretization of the DtN as described in Section 3. The L 2 -error after one application of the preconditioner with respect to a solution computed with a direct solver is shown in Table 3. Apparently, the preconditioner can be used as a direct solver. Analogous results are obtained for the academic example from Section 4.1.1. This confirms the theoretical result derived in Proposition 3.1. Setting v(r, θ) = v PREM (r) (1 + εv pert (r, θ)) results in a relative perturbation of strength ε. In the following, the tensor product DtN based on v PREM is employed to precondition the system for v(r, θ).
In Table 4 the iteration numbers obtained with the tensor product DtN are compared to the moving PML 2 based approximation. The moving PML yields iteration numbers which are essentially independent of the strength of the perturbation since it approximates the DtN based on the perturbed shear velocity. The tensor product DtN instead is obtained from the separable background velocity. Thus, iteration numbers grow with the strength of the perturbation.
The performance of both approaches is highly dependent on the boundary conditions imposed at the Earth's surface. The tensor product DtN is very robust against perturbations for the PML boundary condition. As a result, it outperforms the moving PML method for perturbations up to 2%. When the PML is replaced by a free surface boundary condition then the DtN apparently becomes very sensitive to perturbations. Hence, the tensor product DtN only yields acceptable iteration numbers in the high frequency regime for perturbations smaller than 0.1%.

Remark 2
Even in cases where the model is perfectly captured in the solution of the exterior problem in terms of the data, i.e. assuming no perturbations in the model coefficients, one often needs to deal with numerical perturbations. For instance, we considered problems without data perturbation where, however, adaptive meshes that violate the tensor product structure or for the Table 4 Comparison of GMRES iteration numbers x / y: The numbers x are obtained when the tensor product DtN based on the background model is used to set up a preconditioner for the system with a perturbed shear velocity. For the case y the moving PML approximation of the DtN is employed. The strength of the perturbation ε is given in %. A random source is used for the right hand side (set to zero where PML is present). The dash ' -' means that the desired tolerance of 10 −7 was not reached after 1000 iterations. solution of (10) different ODE solvers have been used. These numerical perturbations resulted in the same effect, a dramatic amplification thereof, rendering the approach practically useless in the presence of reflections.
The tensor product DtN approximation considered above realizes an exact solution of the exterior problem in the case of no perturbation of the background model. The missing robustness of this approach is not due to the tensor product construction, but rather due to the high sensitivity of the exact DtN operator itself. The investigation of this effect will be the subject of the subsequent section.

On the sensitivity of DtN operators
The experiments from the previous section suggest that the DtN for the free surface boundary condition is much more sensitive to perturbations than for an absorbing boundary condition. Performing a mode-by-mode analysis, see section 5.1, confirms this and reveals that the issue is already observed in one dimension. This leads to an analytical sensitivity analysis of scattering problems on the half line via Riccati equations in Section 5.2. Additional illustrations for this case are provided in Section 5.3.

Modal analysis
If in the study above in Section 4.2.1 we replace the perturbation model (13) with only linear velocity perturbations, i.e. v(r) = v PREM (r) (1 + ε) with constant ε ≥ 0, then the perturbed problem is separable as well. Denote by DtN j (ψ , ε) the DtN numbers as computed in Proposition 3.1 for perturbation ε. We fix the interface j = J − 1 closest to the CMB. In Figure 4 the relative DtN error with ε ≈ 3.9 · 10 −5 is shown for ω ∈ {512, 1024, 2048}. The same mesh, in particular the same discrete modes, have been used for all ω.
The relative error on the guided modes |λ | ω 2 for the free surface boundary condition is highly oscillatory and even in the best case almost two orders of magnitude larger than for the PML boundary condition. It also grows linearly in ω and ε. The first statement can be seen in the Figure while the second was observed in other experiments not shown here. These findings show that the DtN for the free surface boundary condition is indeed highly sensitive to perturbations which explains the poor performance of the tensor product DtN based on the background model as observed in Section 4.2.1.

Analysis via Riccati equations
The high sensitivity already occurs in one dimension as demonstrated by the modal analysis. Hence, to get to the core of the problem, we will analyze the one dimensional scattering problems with transparent boundary conditions for a > 0 and ω ≥ 0. Let us consider the functions which may be interpreted as "local DtN numbers". It is well-known and easy to check that these functions satisfy the Riccati equation Moreover, we have the initial conditions For simplicity, let us switch to the perturbation E(x) := ω 2 ε(x). The Fréchet derivatives (formally) satisfy the linear initial value problems which follows by differentiating (16) with respect to E 3 . Our aim is to compare the sizes of k T (0) and k R (0). The solutions to the initial value problems (18) can be expressed in terms of the solutions The main difference between transparent and reflecting boundary conditions is that v T is typically close to purely imaginary (for E = 0 it is identically iω) whereas v R for real-valued ω and E is real-valued and has singularities at zeros of u R . For E = 0 we have v R (x) = −ω cot(ω(a − x)). Therefore, |k y T | is close to 1 whereas the kernel |k y R | has singularities, and its modulus is often much larger than 1. In view of (19), this explains, why |k R (0)| is typically much greater than |k T (0)|. This conclusion is consistent with the relative DtN error for the guided modes as shown in Figure 4. In contrast, |λ | ω 2 leads to exponentially decaying modes and removes the singularities from the corresponding kernel. In this case, |k R (0)| and |k T (0)| are of comparable size.

Comparison of DtN numbers for constant perturbation
Let us illustrate the qualitative statements of the previous section by considering only constant perturbations ε(x) = ε which allows us to compute the solutions to (15-T) and (15-R). These problems become constant coefficient problems in the perturbed wavenumber ω ε = ω √ 1 + ε. The solutions for ω > 0 for transparent and reflecting boundary conditions u ε R (x) = cos(ω ε x) − cot(ω ε a) sin(ω ε x) yield the perturbed DtN numbers In Figure 5 these functions are plotted for constant ε = 10 −3 and a = 1. In this case, the relative error for the transparent boundary conditions can be bounded independently of ω while the relative error for reflecting boundary conditions is much larger, highly oscillatory and grows like ω.

Conclusion
This paper investigates the potential of sweeping preconditioners for stratified media in absence of an absorbing boundary condition. For such a problem the DtN cannot be reasonably approximated by a moving PML. To resolve this issue, a tensor product discretization of the DtN, which is based on separability of the equation for the background model, has been introduced yielding a direct solver for the unperturbed background model. Despite its perfect approximation of the DtN the applicability of the resulting sweeping preconditioner is limited due to a very high sensitivity of the DtN to perturbations. As a conclusion we can state that in the presence of reflections any sweeping preconditioner relying on an accurate, but not perfect, DtN approximationbased on a tensor product structure or not -is doomed to fail in practice.