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.



There are no comments yet.


page 10


Boundary element methods for Helmholtz problems with weakly imposed boundary conditions

We consider boundary element methods where the Calderón projector is use...

Quantitative bounds on Impedance-to-Impedance operators with applications to fast direct solvers for PDEs

We prove quantitative norm bounds for a family of operators involving im...

A Stable Boundary Integral Formulation of an Acoustic Wave Transmission Problem with Mixed Boundary Conditions

In this paper, we consider an acoustic wave transmission problem with mi...

Approximating the ground state eigenvalue via the effective potential

In this paper, we study 1-d random Schrödinger operators on a finite int...

Modeling electrochemical systems with weakly imposed Dirichlet boundary conditions

Finite element modeling of charged species transport has enabled analysi...

Physics-Informed Machine Learning Method for Large-Scale Data Assimilation Problems

We develop a physics-informed machine learning approach for large-scale ...

On the spectral asymptotics of waves in periodic media with Dirichlet or Neumann exclusions

We consider homogenization of the scalar wave equation in periodic media...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 (SH-waves) between the core mantle boundary (CMB) at and the surface of the Earth at

in the frequency domain is described by the equation




and , subject to boundary conditions . The boundary operator is defined piecewise on with and by


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 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.

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 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.

2.1 Double sweep optimized Schwarz method

Figure 1: Decomposition of the domain into layers.

Let be a non-overlapping 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

where solves

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

where solves

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 one-dimensional 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


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.


for all .

The next proposition shows that the are also eigenfunctions of the discretized DtN map.

Proposition 1

Let with be the solution to


for all .

Then .


Define . Note that . So, if we can show that for all then follows. Using (9) for yields:


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


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 preconditioner111The 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 - / - - / - - / -
Table 1: GMRES iteration numbers x / y for moving PML approximation of the DtN for the academic example. For case y the preconditioner is constructed based on a damped problem with wavenumber for . A random source which is set to zero in the PML layer is used as the right hand side. The dash means that the desired tolerance was not achieved within 200 iterations.

4.1.2 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 . 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.
Figure 2: The real part of the solution for with a Dirac source. A PML was applied at the Earth’s surface for (a) while (b) uses a free surface condition.
Figure 3: GMRES residuals for moving PML approximation of the DtN with a Dirac source. On the left a free surface boundary condition is used at the Earth’s surface while on the left a PML is implemented.
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
Table 2: Influence of the boundary condition on GMRES iteration numbers for moving PML approximation of the DtN. The case denoted by ’PML’ uses a radiation condition at the Earth’s surface implemented by a PML while ’free’ sets the free surface boundary condition. Two different sources are considered: a Dirac source (left columns) and a random source which is set to zero in the PML layer (right columns).

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
Table 3: Relative -error after one application of DOSM with tensor product DtN transmission conditions for a tensor product discretization. Two different sources are considered: a Dirac source (left column) and a random source (right column).

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


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 PML222The 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.

Free surface condition at Earth’s surface
256 3 1/  18 4/  18 4/  18 6/  18 7/  18 8/  18 11/  18 512 6 1/  53 9/  53 12/  53 17/  53 26/  53 39/  53 62/  53 1024 12 2/ 119 16/ 120 23/ 120 35/ 122 53/ 122 84/ 118 158/ 123 2048 24 2/ 292 38/ 292 102/ 293 215/ 274 422/ 291 735/ 238 -/ 274 PML at Earth’s surface
256 3 1/ 7 3/ 7 3/ 7 3/ 7 3/ 7 4/ 7 4/ 7 5/ 7 512 6 1/ 9 3/ 9 4/ 9 4/ 9 4/ 9 5/ 9 6/ 9 7/ 9 1024 12 1/ 13 4/ 14 5/ 14 6/ 14 7/ 14 8/ 14 10/ 14 14/ 13 2048 24 2/ 15 5/ 15 6/ 15 7/ 15 8/ 16 11/ 16 15/ 16 22/ 16

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 was not reached after 1000 iterations.

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 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.

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


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.

Figure 4: The relative DtN error on the discrete eigenmodes as computed in (14) for a linear velocity perturbation of fixed strength . The blue line is for a free surface boundary condition at the Earth’s surface while the red line is obtained with an absorbing boundary condition implemented by a PML. The modes have been ordered by increasing magnitude of .

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

and reflecting

boundary conditions for and .

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 . The Fréchet derivatives

(formally) satisfy the linear initial value problems


which follows by differentiating (16) with respect to 333

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


The main difference between transparent and reflecting boundary conditions is that is typically close to purely imaginary (for it is identically ) whereas for real-valued and is real-valued 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 (15-T) and (15-R). 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


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 .

Figure 5: Comparison of analytic DtN errors for a constant perturbation and . Here, and denote reflecting and transparent boundary conditions respectively.

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.

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 SH-waves 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 SH-waves in the time domain assuming axial symmetry is given as


Multiplying with and transforming to the frequency domain we obtain


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:


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