Devito

Devito is a Python package for optimized stencil computation (e.g., finite differences, image processing, machine learning) using symbolic specification.

Key capabilities include:

• Automated Code Generation: Generates optimized kernels for CPUs (SIMD, OpenMP), GPUs (OpenACC), and multi-node clusters (MPI).
• Symbolic Math: Built on SymPy for high-level problem definitions.
• Advanced Operators: Support for sparse operators (interpolation), linear operators (convolutions), tensor contractions, and boundary conditions.
• Optimization: Includes an autotuning framework, blocking, and aggressive symbolic transformations for FLOP reduction.
• Integration: Works with NumPy, SymPy, Dask, SciPy, TensorFlow, and PyTorch.

Devito is a Python-based domain-specific language (DSL) and code generation framework for solving systems of partial differential equations (PDEs) through symbolic computation. It allows you to define equations in pure Python using SymPy-powered syntax and automatically generates optimized C/C++ code for various architectures (CPUs with OpenMP/MPI, and accelerators like CUDA, HIP, SYCL, and OpenACC).

Key Distinctions:

• It is NOT a seismic modeling framework: While often used in seismic imaging, Devito is physics-agnostic. It is a general-purpose PDE solver engine.
• It is NOT a set of pre-baked solvers: You provide the physics and the numerical schemes; Devito provides the high-performance implementation.

When writing code intended for MPI execution, you must distinguish between the global domain shape and the local domain shape (the data shape after decomposition).

• Always use grid.shape (or f.grid.shape for a function f) to refer to the actual global dimensions of the problem.
• Never use f.shape or f.data.shape to determine global dimensions, as these return the local shape specific to an MPI rank.

Devito uses two primary ways to express finite difference operators:

1. f() notation (Functional): The standard, safe way to express derivatives and stencil expressions. It is designed to ensure only legal stencil expressions are built.
2. f[] notation (Indexed): Used for accessing specific grid points or escaping the constraints of the functional language. Warning: This notation is not safe and can easily lead to out-of-bounds memory accesses. Use it judiciously.

You can mix both notations when constructing equations (Eq).

>>> from devito import Grid, Function
>>> grid = Grid(shape=(5, 6))
>>> f = Function(name='f', grid=grid, space_order=2)
>>> f.dx
Derivative(f(x, y), x)
>>> f.dx.evaluate
-f(x, y)/h_x + f(x + h_x, y)/h_x

>>> x, y = grid.dimensions
>>> f[x + 1000, y]
f[x + 1000, y]

>>> f.dx + f[x + 1000, y]
Derivative(f(x, y), x) + f[x + 1000, y]

When subclassing complex types like Function or SparseTimeFunction, you must handle Devito's 'reconstruction' mechanism. Because Devito objects are immutable, symbolic transformations create new objects.

To allow Devito to correctly rebuild your custom objects during compilation, your class must inherit from Reconstructable and define __rargs__ (for positional arguments) and __rkwargs__ (for keyword arguments) that match your __init__ method.

from devito import Reconstructable

class Foo(Reconstructable):
    __rargs__ = ('a', 'b')
    __rkwargs__ = ('c',)

    def __init__(self, a, b, c=4):
           self.a = a
           self.b = b
           self.c = c

    def __repr__(self):
        return "x(%d, %d)" % (self.a, self.b)

# Usage of reconstruction:
a = Foo(3, 5)
a._rebuild(c=5) # Returns a new object: x(3, 5) with c=5

Devito applies several symbolic optimizations to reduce the 'operation count' (FLOPs) in complex expressions. Key optimizations include:

• Factorization: Attempting to factor out common terms to reduce redundant floating-point operations.
• Common Sub-expression Elimination (CSE): Identifying and reusing identical computations.
• Loop-invariant hoisting: Moving computations out of loops if they don't change across iterations.
• Cross-iteration redundancy detection: Optimizing high-order derivatives.

While Devito aims for optimal factorization, some missed opportunities may exist. However, these typically only result in slightly higher runtimes if the code is not heavily compute-bound, and numerical outputs remain indistinguishable to machine precision.