Devito (v3.1.0): an embedded domain-specific language for finite differences and geophysical exploration

Abstract

We introduce Devito, a new domain-specific language for implementing high-performance finite difference partial differential equation solvers. The motivating application is exploration seismology where methods such as Full-Waveform Inversion and Reverse-Time Migration are used to invert terabytes of seismic data to create images of the earth's subsurface. Even using modern supercomputers, it can take weeks to process a single seismic survey and create a useful subsurface image. The computational cost is dominated by the numerical solution of wave equations and their corresponding adjoints. Therefore, a great deal of effort is invested in aggressively optimizing the performance of these wave-equation propagators for different computer architectures. Additionally, the actual set of partial differential equations being solved and their numerical discretization is under constant innovation as increasingly realistic representations of the physics are developed, further ratcheting up the cost of practical solvers. By embedding a domain-specific language within Python and making heavy use of SymPy, a symbolic mathematics library, we make it possible to develop finite difference simulators quickly using a syntax that strongly resembles the mathematics. The Devito compiler reads this code and applies a wide range of analysis to generate highly optimized and parallel code. This approach can reduce the development time of a verified and optimized solver from months to days. arXiv:1808.01995v1 [cs.DM] 6 Aug 2018 maintained or extended. In contrast, the use of high-level abstractions and symbolic reasoning provided by domain-specific languages (DSL) can significantly reduce the time it takes to implement and verify individual operators for use in inversion frameworks, as has already been shown for the finite element method (Logg et al., 2012;Rathgeber et al., 2015;Farrell et al., 2013).State-of-the-art seismic imaging is primarily based upon explicit finite difference schemes due to their relative simplicity and ease of implementation (Virieux, 1986;Symes, 2015a;Weiss and Shragge, 2013). When considering how to design a DSL for explicit finite difference schemes, it is useful to recognize the algorithm as being primarily a sub-class of stencil algorithms or polyhedral computation (Henretty et al., 2013; Andreolli et al., 2015;Yount, 2015). However, stencil compilers lack two significant features required to develop a DSL for finite differences: symbolic computational support required to express finite difference discretizations at a high level, enabling these expressions to be composed and manipulated algorithmically; support for algorithms that are not stencil-like, such as source and receiver terms that are both sparse and unaligned with the finite difference grid. Therefore, the design aims behind the Devito DSL can be summarized as:create a high-level mathematical abstraction for programming finite differences to enable composability and algorithmic optimization, insofar as possible use existing compiler technologies to o...

Time-domain sparsity-promoting least-squares migration with source estimation

“…Compared to earlier work, addressing the memory demands of (LS-)RTM via on-fly-Fourier transforms (Witte et al, 2019c), we tackle the problem of on-the-fly source estimation in the time domain. Because we use industry-strength time-domain finitedifference propagators provided by Devito (Luporini et al, 2018;Louboutin et al, 2019) and exposed in the Julia programming language by JUDI (Witte et al, 2019b), our approach scales in principle to large 3D industrial problems. While we address the importance of estimating source function, we believe that the sensitivity of LS-RTM to errors in the background velocity model needs to be studied as well albeit early work on time-harmonic LS-RTM showing some robustness with respect to these errors (Tu and Herrmann, 2015a).…”

Time‐domain sparsity promoting least‐squares reverse time migration with source estimation

“…To keep our time‐domain wave‐equation solvers with finite differences numerically stable (in our implementation, we used Devito, https://www.devitoproject.org) for our time‐domain finite‐difference simulations and gradient computations (Luporini et al ., 2018; Louboutin et al ., 2019) and JUDI (https://github.com/slimgroup/JUDI.jl) as an abstract linear algebra interface to our algorithms (Witte et al ., 2019b), we introduce an initial guess for the source‐time function with a bandwidth‐limited spectrum that is flat over the frequency range of interest. Under some assumptions on the source‐time function, we can write the true source‐time function as the convolution between the initial guess and the unknown filter – that is, we have where the symbol denotes the temporal convolution.…”

A dual formulation for time-domain wavefield reconstruction inversion

“…The implementation of the method described in this paper leverages on the Devito framework, which allows the automatic generation of highly-optimized finite-difference C code, starting from a symbolic representation of the wave equation Louboutin et al, 2018). In Julia, we interface to Devito through the JUDI package (Witte et al, 2019).…”