A parallel sweeping preconditioner for heterogeneous 3D Helmholtz equations

Introduction

While definite elliptic partial differential equations can be efficiently solved by a wide variety of techniques (e.g., multigrid, ILU, or structured matrix factorizations), indefinite elliptic equations tend to be more challenging. This paper is concerned with three-dimensional heterogeneous Helmholtz equations of the form,

$$\mathcal{A} u \equiv \left[-\Delta - \frac{\omega^2}{c^2(x)}\right] u(x) =f(x),$$ (1)

where $c(x)$ is the spatially varying wave speed, and $u(x) e^{-i \omega t}$ is the time-harmonic response to an acoustic wave equation with forcing function $f(x) e^{-i \omega t}$ . It is important to recognize that $-\Delta$ is positive-definite and that its combination with the negative-definite $-\frac{\omega^2}{c^2}$ term results in an indefinite system.

Before discussing the overall asymptotic complexity of solution techniques, it is helpful to first motivate why high frequency problems require large numbers of degrees of freedom: Given the wave speed bounds $c_{\text{min}} \le c(x) \le c_{\text{max}}$ , we can define the minimum wavelength as $\lambda_{\text{min}} = 2 \pi c_{\text{min}}/\omega$ . In order to resolve oscillations in the solution using piecewise polynomial basis functions, e.g., with finite-difference and finite-element methods, it is necessary to increase the number of degrees of freedom in each direction at least linearly with the number of wavelengths spanned by the domain. In order to combat pollution effects (4), which are closely related to phase errors in the discrete solution, one must use asymptotically more than a constant number of grid points per wavelength with standard discretization schemes. Nevertheless, it is common practice to resolve the domain to as few as ten points per wavelength. In any case, piecewise polynomial discretizations require $\Omega(\omega^d)$ degrees of freedom in $d$ dimensions.

Until recently, doubling the frequency of Eq. (1) not only increased the size of the linear system by at least a factor of $2^d$ , it also doubled the number of iterations required for convergence with standard preconditioned Krylov methods (18,6,17). Thus, denoting the number of degrees of freedom in a three-dimensional finite-element or finite-difference discretization as $N = \Omega(\omega^3)$ , every linear solve required $\Omega(\omega^4)$ work with traditional iterative techniques. Engquist and Ying recently introduced two classes of sweeping preconditioners for Helmholtz equations without large cavities (14,15): Both approaches approximate a block $LDL^T$ factorization of the Helmholtz operator in block tridiagonal form in a manner which exploits a radiation boundary condition. The first approach performs a block tridiagonal factorization algorithm in $\mathcal{H}$ -matrix arithmetic (25,22), while the second approach approximates the Schur complements of the factorization using auxiliary problems with artificial radiation boundary conditions. Though the $\mathcal{H}$ -matrix sweeping preconditioner has theoretical support for two-dimensional problems (14,28), there is not yet justification for three-dimensional problems.

This paper therefore focuses on the second approach, which relies on multifrontal factorizations (21,27,12,34) of the approximate auxiliary problems in order to achieve an $O(\gamma ^2 N^{4/3})$ setup cost and an $O(\gamma N \log N)$ application cost, where $\gamma(\omega)$ denotes the number of grid points used for each Perfectly Matched Layer (PML) (29). While the sweeping preconditioner is competitive with existing techniques even for a single right-hand side, its main advantage is for problems with large numbers of right-hand sides, as the preconditioner appears to converge in $O(1)$ iterations for problems without large cavities. Thus, after setting up the preconditioner, typically only $O(\gamma N \log N)$ work is required for each solution.

A parallel sweeping preconditioner for heterogeneous 3D Helmholtz equations