Source-linked AI summary
Devito (v3.1.0): an embedded domain-specific language for finite differences and geophysical exploration
Mathias Louboutin, Michael Lange, Fabio Luporini, Navjot Kukreja, Philipp A. Witte, Felix J. Herrmann, Paulius Velesko, Gerard J. Gorman
TL;DR
Exploration-seismology inversion requires costly wave-equation computation and repeated development of optimized propagators as formulations and architectures change. Devito embeds a finite-difference DSL in Python with SymPy and compiles symbolic equations into optimized parallel code. The paper reports that the approach decouples solver development from performance optimization while maintaining generality and good computational performance, including applications beyond the motivating seismic problem.
Problem
Seismic inversion requires computationally demanding wave-equation and adjoint solves, while developing and optimizing new propagators for changing formulations and architectures is costly.
Method
Devito embeds a finite-difference DSL in Python, uses SymPy for symbolic mathematical expressions, and compiles them into optimized parallel code.
Results
Devito successfully decouples high-level solver development from low-level performance optimization while delivering good computational performance and maintaining generality.
Takeaways & Limitations
Devito supports concise, parameterized wave-propagator development for seismic inversion and can also be used for problems based on other equations.
Takeaways & Limitations
More advanced algorithms and constraints are necessary for more complex inversion problems involving less accurate initial models or noisy data.
Abstract
from arXiv · showhide
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.
1 Introduction
Devito addresses the computational and development costs of seismic inversion by embedding finite-difference programming in Python and using symbolic compilation to generate optimized parallel code. Its design combines high-level mathematical expression, compiler optimization of affine loops, and extensions for non-stencil operations.
- Motivation: Exploration-seismology inversion is computationally demanding, while productionizing new wave propagators can take many person-months or person-years.Traditional workflows also make mathematical software difficult to port across architectures.
- Problem: Existing stencil compilers lack symbolic support for composing finite-difference discretizations and flexibility for sparse, unaligned source and receiver terms.These limitations hinder their use in complex seismic-imaging workflows.
- Design aims: Devito aims to provide a high-level mathematical abstraction, optimize affine loop nests with compiler technologies, and extend support to non-affine computation.Source and receiver terms are examples of non-affine components.
- Implementation: Embedding the DSL in Python and leveraging SymPy generates an abstract syntax tree from which optimized and parallel C or stencil-DSL code can be produced.Just-in-time compilation permits static and dynamic analysis during code generation.
- Capabilities: Symbolic manipulation, code generation, and JIT compilation let users define wave propagators in a few Python lines and change spatial discretization order through a single parameter.The paper applies these capabilities to workflows including Full-Waveform Inversion and Reverse-Time Migration.
2 Background
The background reviews the evolution from general-purpose compilers and mathematical languages to PDE and stencil DSLs. It identifies persistent limitations in abstraction, extensibility, workflow integration, and portability that motivate Devito.
- Compiler evolution: Compiler development has repeatedly raised abstraction levels to manage increasing program and hardware complexity while automating performance optimization.Domain-specific languages emerged by restricting problem domains to enable more specialized optimizations.
- Motivation: Compiler techniques are presented as a way to replace much of the manual labor required to port optimized codes across rapidly evolving computer architectures.The motivation is especially relevant where improving runtime performance on a particular platform is non-trivial.
- Existing DSLs: PDE and stencil DSLs provide higher-level formulations, but prior finite-difference systems have had limited use in seismic imaging and broader scientific engineering.The paper attributes this limitation to factors beyond technology inertia.
- Abstraction limits: Many stencil DSLs abstract only polyhedral-like affine loops, leaving developers to write complex numerical kernels manually in low-level languages.Hand-tuned wave-propagator codes can reach thousands of lines and become time-consuming and error-prone.
- Extensibility: Existing DSLs often lack flexibility for sparse source terms and interpolation, limiting extension beyond their original scope and integration into complex workflows.Standalone DSLs may also omit composability with external optimization packages such as SciPy.
- Seismic frameworks: Seismic frameworks offer varying abstraction and equation coverage, but C++ implementations can limit extensibility while subroutine-based systems restrict high-level flexibility and interfacing.Some frameworks also limit architecture portability through C or Fortran implementations.
3 Symbolic definition of finite difference stencils with Devito
Devito represents finite-difference wave-equation operators symbolically, combining mathematical notation with data-carrying objects and automated code generation. Its compiler lowers these expressions into optimized C/OpenMP implementations while supporting multidimensional grids, temporal discretization, and additional PDE abstractions.
- Motivation: Devito targets the computational workload of solving discrete wave equations and adjoints in seismic inversion.The framework is motivated by the computational demands of seismic imaging and by the need to improve development and innovation time.
- Design objectives: The DSL provides concise symbolic definitions of explicit, matrix-free finite-difference operators with arbitrary spatial discretization order.The acoustic wave equation serves as a running example, with squared slowness, wavefield, and source represented in the formulation.
- Symbolic objects: Devito symbols combine SymPy-compatible symbolic manipulation with user data and compiler metadata such as grid, dimensions, and numerical type.Function objects represent state variables and associated data, while Grid defines the computational domain and discrete function shape.
- Compilation and optimization: The generated C99/OpenMP code executes from Python and interoperates with NumPy-based scientific software.Devito also manages aligned memory and NUMA-related optimizations through its own allocators.
- Compilation and optimization: Devito’s compiler progressively lowers symbolic expressions to semantically equivalent C code through multiple optimization passes.The symbolic engine applies transformations such as common-subexpression elimination, factorization, and loop-invariant code motion; the loop engine applies vectorization, parallelization, and blocking.
- Symbolic objects: Grid-derived dimensions allow the same symbolic operator definitions to support two- and three-dimensional problems.Function objects inherit the spatial dimensionality from the Grid, while TimeFunction adds a time dimension for buffered or stored histories.
- Stencil construction: Devito provides shorthand derivative and Laplacian expressions and a secondary high-level abstraction for PDE features not captured by basic finite-difference abstractions.Examples include u.dx, u.dx2, and dimension-agnostic u.laplace notation.
4 Seismic modeling and inversion
Seismic inversion reconstructs subsurface properties by minimizing mismatch between measured and modeled wavefields. Devito expresses forward, adjoint, and gradient operators symbolically, including checkpointing strategies for memory-intensive FWI.
- Seismic inversion reconstructs subsurface properties by minimizing the misfit between recorded measurements and numerically modeled seismic data.
- FWI formulates wave-speed recovery as a PDE-constrained optimization problem parameterized by velocity or square slowness.
- Discrete adjoint wave-equation and FWI gradient: Computing gradients requires both forward and adjoint wavefields, creating substantial memory demands that checkpointing reduces by recomputing forward states.
- Acoustic forward modelling operator: Devito generates forward-modeling stencils that update the wavefield, inject sources, and interpolate values at receiver locations.
- Discrete adjoint wave-equation and FWI gradient: The adjoint propagator reverses time, injects residual data at receiver locations, and samples the wavefield at the original source location.
- Discrete adjoint wave-equation and FWI gradient: Optimal checkpointing reduces stored history but requires log(nt)+2 extra PDE solves, while boundary reconstruction requires one extra forward solve.
5 Verification
Devito-generated operators are verified against analytical solutions and adjoint-based tests. Temporal and spatial convergence follow expected rates, supporting the correctness of the symbolic discretization before inversion tests.
- Verification tests assess wave-equation convergence, the discrete adjoint, and the computed gradient.
- The forward operator is compared with an analytical constant-medium solution using accuracy-versus-grid-size and accuracy-versus-time measures.
- The tests assume a sufficiently large domain and short propagation time so the wavefield does not reach the boundaries.
- Temporal convergence: 1.94 slope: temporal error decreases near quadratically with time-step size, matching the theoretical second-order rate.
- Spatial convergence: Spatial convergence follows theoretical rates for orders k = 2,4,6,8, with an exception at order 10 for small grid sizes due to numerical-accuracy limits.
- These time and space results verify the accuracy and correctness of Devito’s symbolic discretization and support subsequent operator verification.
5.2 Propagators verification for inversion
Devito verifies its forward, adjoint, and gradient operators for seismic inversion through dot and gradient tests. The tests show numerical-precision adjoints and theoretically matching Taylor-expansion slopes.
- Adjoint test: The dot test verifies the discrete adjoint operator in both 2D and 3D across discretization orders.The adjoint is accurate up to numerical precision, with error of order 1e−16.
- Adjoint test: The adjoint propagator is the exact numerical adjoint of the forward propagator and implements the adjoint wave equation.
- Gradient test: The gradient test evaluates the FWI objective gradient against its Taylor expansion by varying perturbation scale h.The test uses a small model perturbation and varies h between 10−6 and 100.
- Gradient test: 1.06 and 2.01 are the observed numerical convergence orders, confirming correct implementation of the FWI operators.
- Gradient test: The matching error slopes with theoretical h and h2 slopes verify the accuracy of the inversion operators.The gradient test was run in double precision with an 8th-order spatial discretization.
5.3 Validation: Full-Waveform Inversion
Devito is validated on a Marmousi-ii full-waveform inversion example using forward and adjoint propagators within a seismic inversion framework. After 15 iterations, the result provides a proof of concept while leaving more complex data and initialization settings for advanced algorithms.
- Marmousi-ii validation: The Marmousi-ii example reconstructs the true velocity model after 15 inversion iterations.The model spans 4km × 16km with a 10m grid, using a 10Hz Ricker wavelet and 4s recording.
- Scope: The Marmousi-ii result is presented as a good proof of concept for Devito.
- Scope: The wave propagator remains unaffected by the need for more advanced algorithms and constraints in harder inversion settings.Such settings include less accurate initial models, noisy data, and field-recorded data.
- Framework integration: Devito implements the propagators and FWI gradient in a few concise lines with a mathematical syntax.
- Framework integration: The example uses JUDI as a high-level abstraction layer over Devito for seismic inversion and linear algebra.
5.4 Computational Fluid Dynamics
Three computational-fluid-dynamics examples demonstrate Devito’s flexibility beyond wave propagation: convection, Burgers’ equations, and Poisson problems. The examples cover time-dependent, nonlinear coupled, and steady-state PDEs with solutions consistent with reference results where reported.
- Examples: The CFD examples include convection, Burgers’ equations, and Poisson equations implemented with Devito.
- Convection: The convection example discretizes time with forward differences and space with backward differences.Its solution evolution is consistent with the expected reference result.
- Burgers’ equation: The Burgers’ example demonstrates support for coupled and nonlinear equations with solution components u and v.The diffusion coefficient is denoted by ν, and the resulting solution corresponds to the reference solution.
- Poisson equation: The Poisson example solves a steady-state problem iteratively and provides two implementations producing the same result.The second implementation uses BufferedDimension to alternate automatically between time buffers.
- Conclusion: Together, the examples cover nonlinear equations, coupled PDE systems, and steady-state problems.
6 Performance
Devito’s performance evaluation examines discretization accuracy, runtime, and hardware utilization. Higher-order discretizations can reduce cost through coarser grids, while optimized runs achieve good portable performance but show inefficiencies at 16th order.
- Accuracy-cost trade-off: Higher-order discretizations produce more accurate solutions on coarser grids with smaller runtime.A grid size two times bigger implies a reduction of memory usage by a factor of 24 for 3D modeling.
- Accuracy-cost trade-off: Devito lets users choose discretization order and grid size according to desired accuracy and available computational resources.
- Roofline analysis: 17.1sec for a 512×512×512 domain and 162.6sec for a 1024×1024×1024 domain illustrate near-linear runtime growth with domain size.The larger domain is 8 times bigger and about 9 times slower.
- Roofline analysis: The computation is memory-bound because the working sets do not fit in the L3 cache.
- Roofline analysis: 16th-order discretization shows a performance drop relative to other orders and the attainable peak, especially for larger domains.Profiling supports inefficient memory usage, particularly misaligned data accesses, as the hypothesis for this behavior.
- Portability: Devito reaches 50-60% of the attainable peak, a level described as usually very good, while retaining portable application code across architectures.Testing new architectures requires compiler extensions, if any, while application code remains unchanged.
7 Future Work
Devito separates mathematical formulation from performance engineering, enabling quicker development, maintenance, portability, and performance across problems. Future work targets broader PDE support, large-scale computation, specialized backends, and applications still under development.
- Future Work: Compiler technology can be reused across problems, allowing HPC developer effort to focus on improvements that keep pace with technological evolution.The paper describes software reuse as fundamental to adapting to changing technology.
- Future Work: YASK integration provides a specialized backend in which Devito generates and compiles YASK code for Intel architectures.The backend is presented as a route toward long-term performance portability.
- Future Work: Future extensions require support for large-scale computations, more PDE types, staggered grids with half-node discretization, and domain decomposition.These requirements are intended to extend the DSL to a broader community.
- Future Work: Additional applications are still in development and are planned for future releases.The passage does not specify which applications will be added.
8 Conclusions
Devito introduces a Python-embedded DSL for finite-difference time-domain simulation and applies it to seismic inversion. The approach produces highly optimized parallel solvers in few lines while remaining general enough for other PDE-based problems.
- 8 Conclusions: Devito introduces a DSL for time-domain finite-difference simulation and applies it to a seismic inverse problem.The conclusion frames seismic inversion as the primary application.
- 8 Conclusions: A highly optimized and parallel finite-difference solver can be implemented in just a few lines of Python code.This directly links the DSL interface to solver implementation efficiency.
- 8 Conclusions: Devito already supports problems based on equations beyond those required for seismic imaging applications.The conclusion points to computational fluid dynamics examples as evidence of this broader use.
- 8 Conclusions: The traditional code for these problems is highly complex because mathematical and performance-engineering complexity are interleaved.The paper identifies this interleaving as the primary reason for code complexity.
- 8 Conclusions: Devito decouples these concerns while delivering good computational performance and maintaining generality.Both properties are stated as outcomes that will continue to be improved in future versions.
9 Code Availability
The Devito source code, examples, tests, installation instructions, publications, and documentation are available online. The stated installation procedure clones version v3.1.0, creates the Conda environment, and installs the package in editable mode.
- 9 Code Availability: Source code, examples, and test scripts are available on GitHub.The repository includes a README for installation.
- 9 Code Availability: Project documentation and publication information are available through the Devito website and Sphinx-generated documentation.The passage distinguishes the project overview site from the generated documentation.
- 9 Code Availability: The installation procedure clones the v3.1.0 branch, creates the Conda environment, and installs Devito in editable mode.The commands are `git clone -b v3.1.0`, `conda env create -f environment.yml`, and `pip install -e .`.