Closing the performance gap between an iterative frequency-domain solver and an explicit time-domain scheme for 3D migration on parallel architectures

Abstract

Three-dimensional reverse-time migration with the constantdensity acoustic wave equation requires an efficient numerical scheme for the computation of wavefields. An explicit finitedifference scheme in the time domain is a common choice. However, it requires a significant amount of disk space for the imaging condition. The frequency-domain approach simplifies the correlation of the source and receiver wavefields, but requires the solution of a large sparse linear system of equations. For the latter, we use an iterative Krylov solver based on a shifted Laplace multigrid preconditioner with matrixdependent prolongation. The question is whether migration in the frequency domain can compete with a time-domain implementation when both are performed on a parallel architecture. Both methods are naturally parallel over shots, but the frequency-domain method is also parallel over frequencies. If we have a sufficiently large number of compute nodes, we can compute the result for each frequency in parallel and the required time is dominated by the number of iterations for the highest frequency. As a parallel architecture, we consider a commodity hardware cluster that consists of multicore central processing units (CPUs), each of them connected to two graphics processing units (GPUs). Here, GPUs are used as accelerators and not as an independent compute node. The parallel implementation of the 3D migration in frequency domain is compared to a time-domain implementation. We optimize the throughput of the latter with dynamic load balancing, asynchronous I/O, and compression of snapshots. Because the frequency-domain solver uses matrix-dependent prolongation, the coarse-grid operators require more storage than available on GPUs for problems of realistic size. Due to data transfer, there is no significant speedup using GPU-accelerators. Therefore, we consider an implementation on CPUs only. Nevertheless, with the parallelization over shots and frequencies, this approach could compete with the time-domain implementation on multiple GPUs.

Reduction of computing time for least-squares migration based on the Helmholtz equation by graphics processing units

“…The right-hand side is the sum over each source of its subsurface image, that is obtained by migration F H of the data at the receiver position corresponding to the given source. Note that migration in the frequency domain is described in detail in our previous work Knibbe et al [17]. The lefthand side consists of a sum over the forward modeling (6) for a given set of reflectivity coefficients for each source, consecutively followed by the migration.…”

“…Here, the data is transferred to and from the GPU for each new task. This approach has been investigated for the wave equation in the time domain in Knibbe et al [17]. While the simplicity of the time domain algorithm makes it easy to use GPUs of modest size to accelerate the computations, it is not trivial to use GPUs as accelerators for the Helmholtz solver.…”

“…For example, in migration, the wave equation has been traditionally solved by an explicit time discretization scheme in the time domain requiring large amounts of disk space. In Knibbe et al [17], we have shown that solving the wave equation in the frequency domain, i.e., the Helmholtz equation, can compete with a time domain solver given a sufficient number of parallel computational nodes with a limited usage of disk space. The Helmholtz equation is solved using iterative methods.…”

Higher-order Source-wavefield Reconstruction for Reverse-time Migration from Stored Values in a Boundary Strip Just One Point Wide

“…Three-dimensional RTM typically uses higher frequencies than FWI, and storage of source wavefields on disk or solid-state drives or even in main memory may lead to performance bottlenecks, in particular on manycore or GPU hardware. Migration in the frequency domain avoids the storage problem and outperforms migration in the time domain in two dimensions (Marfurt and Shin, 1989;Pratt, 1990;Østmo et al, 2002;Mulder and Plessix, 2004) but not yet in three dimensions (Riyanti et al, 2006;Plessix, 2009;Wang et al, 2010Wang et al, , 2011Knibbe et al, 2014;Amestoy et al, 2016). For 3D applications, there are several ways to reduce or circumvent the storage of the source wavefields (Dussaud et al, 2008;Nguyen and McMechan, 2015).…”

“…For 3D applications, there are several ways to reduce or circumvent the storage of the source wavefields (Dussaud et al, 2008;Nguyen and McMechan, 2015). Data compression with subsampling or with wavelet or Fourier transforms is one approach (Araya-Polo et al, 2011;Sun and Fu, 2013;Knibbe et al, 2014). Checkpointing, storing the wavefield at selected times to recompute the wavefields in small time intervals when needed for correlation with the reverse time computations, is another approach (Griewank, 1992;Symes, 2007).…”