Contraction operator transformation for the complex heterogeneous Helmholtz equation

An efficient solution of the three-dimensional Helmholtz equation is known to be crucial in many applications, especially geophysics. In this paper, we present and test two preconditioning approaches for the discrete problem resulting from the second order finite-difference discretization of this equation. The first approach combines shifted-Laplacian preconditioner with inversion of a separable matrix, corresponding to the horizontally-layered velocity model, using fast Fourier based transforms. The second approach is novel and involves a special transformation resulting in a preconditioner with a contraction operator (CO preconditioner). The two approaches have near the same arithmetical complexity; however, the second approach, developed in this paper, provides a faster convergence of an iterative solver as illustrated by numerical experiments and analysis of the spectral properties of the preconditioned matrices. Our numerical experiments involve parallel modeling of highly heterogeneous lossy and lossless media at different frequencies. We show that the CO-based solver can tackle problems with hundreds of millions of unknowns on a conventional cluster node. The CO preconditioned solver demonstrates a very moderate increase of iteration count with the frequency. We have conducted a comparison of the performance of the developed method versus open-source parallel sweeping preconditioner. The results indicate that, the CO solver is several times faster with respect to the wall-clock time and consumes substantially less memory than the code based on the sweeping preconditioner at least in the example we tested.