Fast finite-difference convolution for 3D problems in layered media

TL;DR

This paper introduces a fast direct solver for 3D Helmholtz and Maxwell equations in layered media, significantly improving computational efficiency over traditional Fourier transform-based methods, especially for non-uniform grids.

Contribution

The authors develop a novel cyclic reduction-based algorithm that reduces computational complexity for layered media problems, applicable to both uniform and non-uniform grids.

Findings

Achieves $O(N_xN_ylog(N_xN_y)N_z)$ complexity for uniform grids.

Outperforms discrete Fourier transform-based algorithms in efficiency.

Extensible to elasticity problems.

Abstract

We developed fast direct solver for 3D Helmholtz and Maxwell equations in layered medium. The algorithm is based on the ideas of cyclic reduction for separable matrices. For the grids with major uniform part (within the survey domain in the problems of geophysical prospecting, for example) and small non-uniform part (PML and coarsening to approximate problems in infinite domain) the computational cost of our approach is $O(N_xN_ylog(N_xN_y)N_z)$. For general non-uniform grids the cost is $O(N^{3/2}_xN^{3/2}_yN_z)$. The first asymptotics coincide with the cost of FFT-based methods, which can be applied for uniform gridding (in x and y) only. Our approach is significantly more efficient compared to the algorithms based on discrete Fourier transform which cost is $O(N^2_xN^2_yN_z)$. The algorithm can be easily extended for solving the elasticity problems as well.