Efficient Helmholtz equation solver for frequency domain waveform inversion based on the decomposition into one-way wave equations

When using a ring array to perform ultrasound tomography, the most computationally intensive component of frequency-domain full waveform inversion (FWI) is the Helmholtz equation solver. The Helmholtz equation is an elliptic partial differential equation (PDE) whose discretization leads to a large system of equations; in many cases, the solution of this large system is itself the inverse problem and requires an iterative method. Our current solution relies on discretizing the 2D Helmholtz equation based on a 9-point stencil and using the resulting block tridiagonal structure to efficiently compute a block LU factorization. Conceptually, the L and U systems are equivalent to a forward and backward wave propagation along one of the spatial dimensions of the grid, resulting in a direct non-iterative solution to the Helmholtz equation based on a single forward and backward sweep. Based on this observation, the PDE representation of the Helmholtz equation is split into two one-way wave equations prior to discretization. The numerical implementations of these one-way wave equations are highly parallelizable and lend themselves favorably to accelerated GPU implementations. We consider the Fourier split-step and phase-shift-plus-interpolation (PSPI) methods from seismic imaging as numerical solutions to the one-way wave equations. We examine how each scheme affects the numerical accuracy of the final Helmholtz equation solution and present its impact on FWI with breast imaging examples.