1 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 EZ12 .
Recently, sweeping preconditioners have emerged as a promising approach to overcome this problem GZ19 .
Since the introduction of the moving perfectly matched layer (PML) preconditioner by Engquist and Ying EY11 , numerous impressive results and further developments
of this technique have been published.
We refer to GZ19 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 DZ16 or section 10 of GZ19 .
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 IW95 ; JTCI08 in which spherical coordinates are used. Assuming axisymmetry of all fields in the propagation of shear waves (SHwaves) between the core mantle boundary (CMB) at and the surface of the Earth at
in the frequency domain is described by the equation
(1) 
where
(2) 
and , subject to boundary conditions . The boundary operator is defined piecewise on with and by
(3) 
Here, is the mass density, is the shear modulus and the frequency. The background coefficients for and are provided by the spherically symmetric PREM model DA81 .
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 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 BB02 .
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 DirichlettoNeumann (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.
2 General framework of sweeping preconditioners
In GZ19 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 nonoverlapping 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.
2.1 Double sweep optimized Schwarz method
Let be a nonoverlapping decomposition of the domain into horizontal layers , see Figure 1. In the numerical examples below, we make the choice that one layer contains two finite elements in vertical direction () and elements in the horizontal direction () yielding a fixed aspect ratio for all finite elements. We denote by the interfaces between the layers and write for a function defined in layer . Based on these definitions we state the specialized sweeping algorithm.
Forward sweep: Given the last iterate in solve successively for
Backward sweep: Solve successively for  
The transmission operator is an (approximation) of the DtN map
(4)  
where solves  
(5)  
with . If is well defined for and the original problem and subdomain problems are uniquely solvable then DOSM converges in one double sweep to the exact solution for GZ19 . In practice, using the exact DtN map as a transmission operator is computationally too expensive. Therefore, is chosen as an approximation of the DtN. The algorithm is then usually used as a preconditioner for GMRES with initial guess .
2.2 Moving PML approximation of the DtN
The problem for calculating the exact is posed on the whole exterior domain . Usually, it is assumed that on an absorbing boundary condition implemented by a PML is present. As the name “moving PML” suggests this PML is shifted closer to . This replaces the original problem posed on by a modified problem on the (usually smaller) domain . In practice, the PML is usually started right at the coupling interface and . This leads to the operator
(6)  
where solves  
(7)  
Here, the original differential operators , have been replaced by modified versions , due to the complex scaling applied in the PML region.
3 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 at on which the shall be computed.
Let be the finite element space in .
A tensor product discretization is assumed.
Hence, , where and are onedimensional finite element spaces.
In particular, is a linear combination of terms of the form with and .
Given Dirichlet data the DtN in the finite element setting is computed as follows.
First find with such that
(8) 
Then .
This problem will be solved in two steps.
In section 3.1 we treat the case in which
is a discrete eigenfunction of the weighted Laplacian on
. 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.3.1 The DtN applied to a discrete eigenfunction
Let for denote the discrete eigenfunctions of the discretized Laplacian on
with eigenvalue
, i.e.(9) 
for all .
The next proposition shows that the are also eigenfunctions of the discretized DtN map.
Proposition 1
Let with be the solution to
(10) 
for all .
Then .
Proof
Define . Note that . So, if we can show that for all then follows. Using (9) for yields:
(11) 
since solves (10).
Remark 1
The finite element space for the ODE (10) has to coincide with the first factor of the tensor product . In other words, the ODE discretization is chosen as the restriction of the two dimensional discretization to the radial direction.
3.2 General Dirichlet data
To apply the DtN to general Dirichlet data we express it in the eigenbasis and apply the DtN as
(12) 
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.
rjasanow1994effective . 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.4 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 th order conforming finite elements and have been implemented in the finite element library Netgen/NGSolve, see JS97 ; JS14 . Scripts for reproducing the numerical results are provided at DOI: 10.5281/zenodo.3686802.
4.1 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.
4.1.1 Academic example
Consider the Helmholtz equation on the unit disk with discontinuous coefficients. The bilinear form in polar coordinates is given by
Periodic boundary conditions in are used.
At an absorbing boundary condition implemented by a PML is set, while at we
simply use natural boundary conditions.
The geometrical setup is as shown in Figure 1.
Similar to the experiments from GZ19 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 , on the next then again continuing in this fashion.
The density is simply .
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}^{1}1The original problem to which we apply the preconditioner is still the one without damping., i.e. the operator in equation (7) describes a Helmholtz problem with complex frequency , where . Adding additional damping to the preconditioner leads to nearly robust iteration numbers for . However, for both versions the iteration numbers increase drastically as the contrast is increased. As a result, the preconditioner becomes completely inefficient.
8  3  8 / 9  9 / 10  9 / 11  8 / 10  8 / 10 

16  6  13 / 11  17 / 17  30 / 32  39 / 40  46 / 45 
32  12  15 / 12  46 / 37  88 / 93  154 / 162   /  
64  24  19 / 13  189 / 185   /    /    /  
4.1.2 SHwaves 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 . 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 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 and .
Despite these difficulties, the preconditioner performs well.
The iteration numbers grow only very mildly and the low tolerance of is achieved for 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.
(a) with PML. 
(b) free surface. 
Dirac source  random source  

# iter (PML)  # iter (free)  # iter (PML)  # iter (free)  
256  3  7  13  7  18 
512  6  10  37  9  53 
1024  12  14  73  13  119 
2048  24  15  203  15  307 
4.2 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 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.
error (Dirac)  error (random)  

256  3  
512  6  
1024  12  
2048  24 
4.2.1 Sound speed perturbations
So far only the radially symmetric shear velocity of PREM has been considered. Let us now introduce a velocity perturbation
(13) 
Setting results in a relative perturbation
of strength .
In the following, the tensor product DtN based on is employed to precondition the system for .
In Table 4 the iteration numbers obtained with the tensor product DtN are compared to the moving PML^{2}^{2}2The results for the moving PML with are slightly different compared to Table 2 because the random source is not set to zero in the uppermost layer. 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 . 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 .
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 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.
5 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 modebymode 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.
5.1 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. with constant , then the perturbed problem is separable as well. Denote by the DtN numbers as computed in Proposition 3.1 for perturbation . We fix the interface closest to the CMB. In Figure 4 the relative DtN error
(14) 
with is shown for . The same mesh, in particular the same discrete modes, have been used for all .
The relative error on the guided modes 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.
5.2 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
(15T)  
and reflecting  
(15R) 
boundary conditions for and .
Let us consider the functions
which may be interpreted as “local DtN numbers”. It is wellknown and easy to check that these functions satisfy the Riccati equation
(16) 
Moreover, we have the initial conditions
(17) 
For simplicity, let us switch to the perturbation . The Fréchet derivatives
(formally) satisfy the linear initial value problems
(18) 
which follows by differentiating (16) with respect to ^{3}^{3}3
As we only aim for a heuristic explanation rather then a theorem in this section, we do not go into the technicalities of justifying this statement rigorously in the presence of singularities of
.. Our aim is to compare the sizes of and .The solutions to the initial value problems (18) can be expressed in terms of the solutions
to the homogeneous equations with initial conditions as
(19) 
The main difference between transparent and reflecting boundary conditions is that is typically close to purely imaginary (for it is identically ) whereas for realvalued and is realvalued and has singularities at zeros of . For we have . Therefore, is close to whereas the kernel has singularities, and its modulus is often much larger than . In view of (19), this explains, why is typically much greater than . This conclusion is consistent with the relative DtN error for the guided modes as shown in Figure 4. In contrast, leads to exponentially decaying modes and removes the singularities from the corresponding kernel. In this case, and are of comparable size.
5.3 Comparison of DtN numbers for constant perturbation
Let us illustrate the qualitative statements of the previous section by considering only constant perturbations which allows us to compute the solutions to (15T) and (15R). These problems become constant coefficient problems in the perturbed wavenumber . The solutions for for transparent
and reflecting boundary conditions yield the perturbed DtN numbers
Comparing this with the background profile gives the relative errors
and  
In Figure 5 these functions are plotted for constant and . 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 .
6 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 approximation – based on a tensor product structure or not – is doomed to fail in practice.
Acknowledgements.
J. Preuß is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen;Author contributions T.H. initiated the research on sweeping preconditioners for problems with reflection. All authors designed and performed the research. Selection of the model and derivation of the discrete formulation of the SHwaves example is due to J.P.. All implementations and numerical computations were performed by J.P. using the finite element library NGSolve. J.P. drafted the paper and all authors contributed to the final manuscript.
Conflict of interest
The authors declare that they have no conflict of interest.
Appendix A Derivation of the variational formulation
In IW95 the equation for SHwaves in the time domain assuming axial symmetry is given as
(A1) 
Multiplying with and transforming to the frequency domain we obtain
(A2) 
To obtain a variational formulation in spherical coordinates this equation has to be integrated taking the volume element into account. Due to axisymmetry the integration over only yields a constant factor, which can be omitted. We obtain the linear form:
(A3) 
To obtain the bilinear form we integrate by parts in and . The boundary terms in vanish because of the free surface boundary condition. Since homogeneous Dirichlet boundary conditions are imposed at and the boundary terms in vanish as well. Let be the subspace of with homogeneous boundary values at and . The variational problem is: Find so that for all with
(A4) 
References
 (1) Becker, T.W., Boschi, L.: A comparison of tomographic and geodynamic mantle models. Geochemistry, Geophysics, Geosystems 3(1) (2002). DOI 10.1029/2001GC000168. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2001GC000168
 (2) Dziewonski, A.M., Anderson, D.L.: Preliminary reference earth model. Physics of the Earth and Planetary Interiors 25(4), 297 – 356 (1981). DOI https://doi.org/10.1016/00319201(81)900467. URL http://www.sciencedirect.com/science/article/pii/0031920181900467
 (3) Engquist, B., Ying, L.: Sweeping Preconditioner for the Helmholtz Equation: Moving Perfectly Matched Layers. Multiscale Modeling & Simulation 9(2), 686–710 (2011). DOI 10.1137/100804644. URL https://doi.org/10.1137/10080464
 (4) Ernst, O.G., Gander, M.J.: Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods, pp. 325–363. Springer Berlin Heidelberg, Berlin, Heidelberg (2012). DOI 10.1007/9783642220616_10. URL https://doi.org/10.1007/9783642220616_10
 (5) Gander, M.J., Zhang, H.: A Class of Iterative Solvers for the Helmholtz Equation: Factorizations, Sweeping Preconditioners, Source Transfer, Single Layer Potentials, Polarized Traces, and Optimized Schwarz Methods. SIAM Review 61(1), 3–76 (2019). DOI 10.1137/16M109781X. URL https://doi.org/10.1137/16M109781X
 (6) Igel, H., Weber, M.: SHwave propagation in the whole mantle using highorder finite differences. Geophysical Research Letters 22(6), 731–734 (1995). DOI 10.1029/95GL00312. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/95GL00312
 (7) Jahnke, G., Thorne, M.S., Cochard, A., Igel, H.: Global SHwave propagation using a parallel axisymmetric spherical finitedifference scheme: Application to whole mantle scattering. Geophysical Journal International 173(3), 815–826 (2008). DOI 10.1111/j.1365246X.2008.03744.x. URL https://doi.org/10.1111/j.1365246X.2008.03744.x
 (8) Rjasanow, S.: Effective algorithms with circulantblock matrices. Linear Algebra and its applications 202, 55–69 (1994)
 (9) Schöberl, J.: NETGEN An advancing front 2D/3Dmesh generator based on abstract rules. Computing and Visualization in Science 1(1), 41–52 (1997)
 (10) Schöberl, J.: C++11 implementation of finite elements in NGSolve. Tech. rep., ASC201430, Institute for Analysis and Scientific Computing; Karlsplatz 13, 1040 Vienna, Austria (2014)
 (11) ZepedaNúñez, L., Demanet, L.: The method of polarized traces for the 2D Helmholtz equation. Journal of Computational Physics 308, 347 – 388 (2016). DOI https://doi.org/10.1016/j.jcp.2015.11.040. URL http://www.sciencedirect.com/science/article/pii/S0021999115007809
Comments
There are no comments yet.