Research

We present a ray-based finite element method for the high-frequency Helmholtz equation in smooth media, whose basis is learned adaptively from the medium and source. The method requires a fixed number of grid points per wavelength to represent the wave field; moreover, it achieves an asymptotic convergence rate of $\mathcal{O}(\omega^{-1/2})$, where $\omega$ is the frequency parameter in the Helmholtz equation. The local basis is motivated by the geometric optics ansatz and is composed of polynomials modulated by plane waves propagating in a few dominant ray directions. The ray directions are learned by processing a low-frequency wave field that probes the medium with the same source. Once the local ray directions are extracted, they are incorporated into the local basis to solve the high-frequency Helmholtz equation. This process can be continued to further improve the approximations for both local ray directions and high-frequency wave fields iteratively. Finally, a fast solver is developed for solving the resulting linear system with an empirical complexity $\mathcal{O}(\omega ^d)$ up to a poly-logarithmic factor. Numerical examples in 2D are presented to corroborate the claims.

Joint work with Hongkai Zhao, Jun Fang, and Jianliang Qiang.

We present a fast iterative solver for Lippmann-Schwinger equation for high frequency waves scattered by a smooth and compactly supported perturbation. The solver is based on the sparsifying preconditioner and the method of polarized traces. The algorithm has two levels: we first precondition the Lippmann-Schwinger equation using the sparsifying preconditioner that relies on solving a sparse linear system, which is then solved fast using a matrix-free variant of the method of polarized traces in a domain decomposition setting using four sweeps in different directions. The assymptotic cost of applying the preconditioner is $\mathcal{O}(N\log{N})$, where $N$ is the number of degrees of freedom. Numerical experiments in 2D indicate that the number of iterations in both levels depends weakly on the frequency.

Joint work with Hongkai Zhao.

We present a fast solver for the 3D high-frequency Helmholtz equation with heterogeneous, constant density, acoustic media. The solver is based on the nested polarized-traces method, coupled with distributed linear algebra libraries and pipelining to obtain a solver with online runtime $\mathcal{O}(\max(1,R/n)N)$ where $N = n^3$ is the total number of degrees of freedom and $R$ is the number of right-hand sides.

The solver has two levels, the outer level, in which we seek the solution in the whole domain; and the inner level, in which we seek the solution of the Helmholtz equation in a subdomain. We present two implementations of the method depending on the choice for the inner solver. We either use high-quality off-the-shelf distributed linear algebra libraries to solve the local problem at each subdomain, or we use the nested variant of the method of polarized traces, in which each layer is decomposed further in smaller subproblems that we call beams, and the local problem is solved iteratively within each layer. In order to speed up the nested approach we use distributed linear algebra to solve the problems within each beam.

Joint work with Laurent Demanet, Russell Hewett and Adrien Scheuer.

We present a solver for the 2D high-frequency Helmholtz equation in heterogeneous acoustic media, with online parallel complexity that scales optimally as $\mathcal{O}(N/L)$, where $N$ is the number of volume unknowns, and $L$ is the number of processors, as long as L grows at most like a small fractional power of N. The solver decomposes the domain into layers, and uses transmission conditions in boundary integral form to explicitly define "polarized traces", i.e., up- and down-going waves sampled at interfaces. Local direct solvers are used in each layer to precompute traces of local Green's functions in an embarrassingly parallel way (the offline part), and incomplete Green's formulas are used to propagate interface data in a sweeping fashion, as a preconditioner inside a GMRES loop (the online part). Adaptive low-rank partitioning of the integral kernels is used to speed up their application to interface data. The method uses second-order finite differences. The complexity scalings are empirical but motivated by an analysis of ranks of off-diagonal blocks of oscillatory integrals. They continue to hold in the context of standard geophysical community models such as BP and Marmousi 2, where convergence occurs in 5 to 10 GMRES iterations.

Another related work is Preconditioning the 2D Helmholtz equation with polarized traces; some extensions are presented in A short note on the nested-sweep polarized traces method for the 2D Helmholtz equation and in Nested domain decomposition with polariced traces for the 2D Helmholtz Equation. An extension to HDG discretizations can be found in A short note on a fast and high-order Hybridazable Discontinuous Galerkin solver for the 2D high-frequency Helmholtz equation.

Joint work with Laurent Demanet.

In this project we study the application of a timestepping method unconstrained by the CFL condition, for computational acoustic wave propagation in the context of seismic waveform inversion. The numerical scheme is a locally one-dimensional (LOD) variant of alternating dimension implicit (ADI) method. Its main feature is that it allows "upscaled timestepping", i.e., the time step is only restricted by the Nyquist sampling rate. The advantage over traditional explicit timestepping methods is threefold: (1) the complexity is automatically lowered in the low-frequency case, without having to coarsen the spatial mesh, (2) the size of the time step is not adversely affected by high velocity contrasts, and (3) perfectly matched layers (PML) can be much steeper, hence much thinner, without a penalty on the time step. The novelty of the project, from a numerical analysis viewpoint, is that a PML is added to the LOD scheme.

Joint work with Laurent Demanet and Russell Hewett.

In this project, we treat the extrapolation of wave-fields using limited aperture Dirichlet data on a broken line. We prove that, in uniform velocity and under geometrical assumptions, it is possible to extrapolate a harmonic wave-field several wavelengths within machine precision. This extrapolation is accurate inside a cone, whose slope we determine using an asymptotic development. We also provide bounds on the extrapolation error. The analysis relies on the properties of the prolate spheroidal wave functions (PSWF). The theoretical results are illustrated by numerical experiments, and we propose a solver for the Helmholtz equation based on this approach.