Devito provides a robust API for sparse operations, encompassing two crucial functionalities: injection
and interpolation
. These operations are fundamental for handling sparse data, such as point sources and measurements, within numerical simulations.
The injection
operation involves injecting values at specified sparse positions, simulating scenarios like point sources in physical systems. Mathematically, if we denote a SparseTimeFunction
as $ S(t, x_s) $, where $ t $ represents time and $ x_s $ denotes sparse spatial coordinates, the injection operation can be expressed as:
Here, $ \delta(x - x_s) $ is the Dirac delta functions, and $ t $ and $ x_s $ represent the time and sparse spatial coordinates of the injection point, respectively.
On the other hand, the interpolation
operation reads the values of a field at sparse positions, mimicking point measurements. If $ F(t, x) $ is a field ('Function') defined on the full grid, the interpolation operation can be expressed as:
Here, $ I(t, r) $ represents the interpolated values at time $ t $ and sparse spatial coordinates $ x_r $.
In Devito, sparse operations are defined as methods for various sparse function types, including SparseFunction
, SparseTimeFunction
, PrecomputedSparseFunction
, and PrecomputedSparseTimeFunction
. For practicality, this tutorial focuses on time-dependent functions, but it's crucial to note that all operations discussed here apply to spatial-only functions as well.
from devito import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator, FixedLocator
grid = Grid((51, 51))
nt = 11
f = TimeFunction(name="f", grid=grid, space_order=4)
u = TimeFunction(name="u", grid=grid, space_order=4, save=nt)
We first define five sparse positions for the purpose of this tutorial. We consider four points in a 2D grid intentionally not aligned with the grid points to highlight the interpolation and injection operations.
#NBVAL_IGNORE_OUTPUT
npoint = 5
coords = np.random.rand(npoint, 2)/2 + .25
base = np.floor(coords / grid.spacing)*grid.spacing
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.scatter(coords[:, 0], coords[:, 1], s=10, label="Sparse positions")
ax.grid(which = "major")
ax.grid(which = "minor", alpha = 0.2)
ax.xaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.yaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.set_title("Off the grid sparse positions")
for i in range(npoint):
ax.annotate("(%.3f, %.3f)" % (coords[i, 0], coords[i, 1]), coords[i, :])
ax.legend()
plt.show()
A SparseFunction
is a devito object representing a Function
defined at sparse positions. It contains the coordinates of the sparse positions and the data at those positions. The coordinates are stored in a SubFunction
object as a Function
of shape (npoints, ndim)
where npoints
is the number of sparse positions and ndim
is the number of dimension of the grid.
A SparseFunction
comes with the two main methods inject(field, expr)
and interpolate(field)
that respectively inject expr
into field
at the sparse positons and interpolate field
at the sparse positions.
print(SparseTimeFunction.__doc__)
Tensor symbol representing a space- and time-varying sparse array in symbolic equations. Like SparseFunction, SparseTimeFunction carries multi-dimensional data that are not aligned with the computational grid. As such, each data value is associated some coordinates. A SparseTimeFunction provides symbolic interpolation routines to convert between TimeFunctions and sparse data points. These are based upon standard [bi,tri]linear interpolation. Parameters ---------- name : str Name of the symbol. npoint : int Number of sparse points. nt : int Number of timesteps along the time Dimension. grid : Grid The computational domain from which the sparse points are sampled. coordinates : np.ndarray, optional The coordinates of each sparse point. space_order : int, optional, default=0 Discretisation order for space derivatives. time_order : int, optional, default=1 Discretisation order for time derivatives. shape : tuple of ints, optional, default=(nt, npoint) Shape of the object. dimensions : tuple of Dimension, optional Dimensions associated with the object. Only necessary if the SparseFunction defines a multi-dimensional tensor. dtype : data-type, optional, default=np.float32 Any object that can be interpreted as a numpy data type. initializer : callable or any object exposing the buffer interface, default=None Data initializer. If a callable is provided, data is allocated lazily. allocator : MemoryAllocator, optional Controller for memory allocation. To be used, for example, when one wants to take advantage of the memory hierarchy in a NUMA architecture. Refer to `default_allocator.__doc__` for more information. Examples -------- Creation >>> from devito import Grid, SparseTimeFunction >>> grid = Grid(shape=(4, 4)) >>> sf = SparseTimeFunction(name='sf', grid=grid, npoint=2, nt=3) >>> sf sf(time, p_sf) Inspection >>> sf.data Data([[0., 0.], [0., 0.], [0., 0.]], dtype=float32) >>> sf.coordinates sf_coords(p_sf, d) >>> sf.coordinates_data array([[0., 0.], [0., 0.]], dtype=float32) Symbolic interpolation routines >>> from devito import TimeFunction >>> f = TimeFunction(name='f', grid=grid) >>> exprs0 = sf.interpolate(f) >>> exprs1 = sf.inject(f, sf) Notes ----- The parameters must always be given as keyword arguments, since SymPy uses ``*args`` to (re-)create the Dimension arguments of the symbolic object.
from itertools import product
s = SparseTimeFunction(name="s", grid=grid, npoint=npoint, nt=nt, coordinates=coords)
d1, d2 = grid.spacing
interp_points = np.concatenate([base+(s1*d1, s2*d2) for (s1, s2) in s._point_support])
#NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.scatter(coords[:, 0], coords[:, 1], s=10, label="Sparse positions")
ax.scatter(interp_points[:, 0], interp_points[:, 1], s=10, label="Interpolation support")
ax.grid(which = "major")
ax.grid(which = "minor", alpha = 0.2)
ax.xaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.yaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.set_title("Off the grid sparse positions")
for i in range(npoint):
ax.annotate("(%.3f, %.3f)" % (coords[i, 0], coords[i, 1]), coords[i, :])
ax.legend()
plt.show()
We first look at the interpolation section and how it is internally lowered into a linear interpolation method.
print("\n".join([str(s) for s in s.interpolate(f).evaluate]))
Eq(posx, (int)(floor((-o_x + s_coords(p_s, 0))/h_x))) Eq(posy, (int)(floor((-o_y + s_coords(p_s, 1))/h_y))) Eq(px, -floor((-o_x + s_coords(p_s, 0))/h_x) + (-o_x + s_coords(p_s, 0))/h_x) Eq(py, -floor((-o_y + s_coords(p_s, 1))/h_y) + (-o_y + s_coords(p_s, 1))/h_y) Eq(sum, 0.0) Inc(sum, (rsx*px + (1 - rsx)*(1 - px))*(rsy*py + (1 - rsy)*(1 - py))*f(t, rsx + posx, rsy + posy)) Eq(s(time, p_s), sum)
f.data.fill(0)
op = Operator([Eq(f.forward, f+1)] + s.interpolate(f))
#NBVAL_IGNORE_OUTPUT
op()
s.data
Operator `Kernel` ran in 0.01 s
Data([[0., 0., 0., 0., 0.], [1., 1., 1., 1., 1.], [2., 2., 2., 2., 2.], [3., 3., 3., 3., 3.], [4., 4., 4., 4., 4.], [5., 5., 5., 5., 5.], [6., 6., 6., 6., 6.], [7., 7., 7., 7., 7.], [8., 8., 8., 8., 8.], [9., 9., 9., 9., 9.], [0., 0., 0., 0., 0.]], dtype=float32)
#NBVAL_IGNORE_OUTPUT
op = Operator(s.inject(u, expr=s))
op()
Operator `Kernel` ran in 0.01 s
PerformanceSummary([(PerfKey(name='section0', rank=None), PerfEntry(time=4.6e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
#NBVAL_IGNORE_OUTPUT
plt.figure(figsize=(10, 10))
plt.imshow(u.data[1], vmin=0, vmax=1, cmap="jet", extent=[0,1,0,1])
plt.colorbar(fraction=0.046, pad=0.04)
plt.title("Linear weights")
plt.show()
Graham J. Hicks pioneered an approach to address one of the significant challenges in seismic simulation: the accurate and flexible positioning of sources and receivers in finite-difference computational grids. By integrating Kaiser windowed sinc functions into the interpolation process, Hicks' method facilitates arbitrary placement of these elements, overcoming the limitations imposed by the inherent discreteness of computational grids. This technique is especially beneficial in seismic wave propagation models, where precise source and receiver locations are crucial for accurate subsurface imaging.
The interpolation of a field value, $f$, at an arbitrary point using the Kaiser windowed sinc function can be mathematically expressed as follows:
$$ f(x) = \sum_{n=-r}^{r} f_n \cdot \text{sinc}\left( \frac{x - x_n}{\Delta x} \right) \cdot w\left( \frac{x - x_n}{\alpha} \right) $$where:
This interpolation formula is pivotal for improving the fidelity of finite-difference simulations in seismic applications. By allowing for the precise positioning of sources and receivers, it enables more accurate modeling of wave propagation through complex geological structures. The flexibility offered by this method significantly enhances the quality of seismic images, providing clearer insights into the Earth's subsurface.
s = SparseTimeFunction(name="s", grid=grid, npoint=npoint, nt=nt, coordinates=coords, interpolation='sinc')
interp_points = np.concatenate([base+(s1*d1, s2*d2) for (s1, s2) in s._point_support])
#NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.scatter(coords[:, 0], coords[:, 1], s=10, label="Sparse positions")
ax.scatter(interp_points[:, 0], interp_points[:, 1], s=10, label="Interpolation support")
ax.grid(which = "major")
ax.grid(which = "minor", alpha = 0.2)
ax.xaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.yaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.set_title("Off the grid sparse positions")
for i in range(npoint):
ax.annotate("(%.3f, %.3f)" % (coords[i, 0], coords[i, 1]), coords[i, :])
ax.legend()
plt.show()
print("\n".join([str(s) for s in s.interpolate(f).evaluate]))
Eq(posx, (int)(floor((-o_x + s_coords(p_s, 0))/h_x))) Eq(posy, (int)(floor((-o_y + s_coords(p_s, 1))/h_y))) Eq(sum, 0.0) Inc(sum, wsincrsx(p_s, rsx + 3)*wsincrsy(p_s, rsy + 3)*f(t, rsx + posx, rsy + posy)) Eq(s(time, p_s), sum)
f.data.fill(0)
op = Operator([Eq(f.forward, f+1)] + s.interpolate(f))
#NBVAL_IGNORE_OUTPUT
op()
s.data
Operator `Kernel` ran in 0.01 s
Data([[0. , 0. , 0. , 0. , 0. ], [0.9967038 , 0.9920868 , 0.99406236, 0.9936394 , 0.99242455], [1.9934076 , 1.9841737 , 1.9881247 , 1.9872788 , 1.9848491 ], [2.9901118 , 2.9762604 , 2.9821863 , 2.9809184 , 2.9772737 ], [3.9868152 , 3.9683473 , 3.9762495 , 3.9745576 , 3.9696982 ], [4.9835196 , 4.9604344 , 4.9703097 , 4.968197 , 4.962124 ], [5.9802237 , 5.952521 , 5.9643726 , 5.961837 , 5.9545474 ], [6.976928 , 6.9446087 , 6.958436 , 6.9554768 , 6.9469705 ], [7.9736304 , 7.9366946 , 7.952499 , 7.9491153 , 7.9393964 ], [8.970336 , 8.9287815 , 8.946562 , 8.942754 , 8.931824 ], [0. , 0. , 0. , 0. , 0. ]], dtype=float32)
#NBVAL_IGNORE_OUTPUT
op = Operator(s.inject(u, expr=s))
op()
Operator `Kernel` ran in 0.01 s
PerformanceSummary([(PerfKey(name='section0', rank=None), PerfEntry(time=3.9999999999999996e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
#NBVAL_IGNORE_OUTPUT
plt.figure(figsize=(10, 10))
plt.imshow(u.data[1], vmin=0, vmax=1, cmap="jet", extent=[0,1,0,1])
plt.colorbar(fraction=0.046, pad=0.04)
plt.title("Sinc weights")
plt.show()
In some cases, simple linear interpolation may not be sufficient for two main reasons:
PrecomputedSparseFunction
offer the interface to answer these two points by allowing the user to provide arbitrary precomputed interpolation weights over an arbitrary large support.
To illustrate this capability, we show in the following how to use PrecomputedSparseFunction
to compute a simple local average over a 4x4
window centered on a sparse point ( average over [x-1, x, x+1, x+2]
in 1D)
print(PrecomputedSparseFunction.__doc__)
Tensor symbol representing a sparse array in symbolic equations; unlike SparseFunction, PrecomputedSparseFunction uses externally-defined data for interpolation. Parameters ---------- name : str Name of the symbol. npoint : int Number of sparse points. grid : Grid The computational domain from which the sparse points are sampled. r : int Number of gridpoints in each Dimension to interpolate a single sparse point to. E.g. `r=2` for linear interpolation. coordinates : np.ndarray, optional The coordinates of each sparse point. gridpoints : np.ndarray, optional An array carrying the *reference* grid point corresponding to each sparse point. Of all the gridpoints that one sparse point would be interpolated to, this is the grid point closest to the origin, i.e. the one with the lowest value of each coordinate Dimension. Must be a two-dimensional array of shape `(npoint, grid.ndim)`. interpolation_coeffs : np.ndarray, optional An array containing the coefficient for each of the r^2 (2D) or r^3 (3D) gridpoints that each sparse point will be interpolated to. The coefficient is split across the n Dimensions such that the contribution of the point (i, j, k) will be multiplied by `interp_coeffs[..., i]*interp_coeffs[...,j]*interp_coeffs[...,k]`. So for `r=6`, we will store 18 coefficients per sparse point (instead of potentially 216). Must be a three-dimensional array of shape `(npoint, grid.ndim, r)`. space_order : int, optional, default=0 Discretisation order for space derivatives. shape : tuple of ints, optional, default=(npoint,) Shape of the object. dimensions : tuple of Dimension, optional Dimensions associated with the object. Only necessary if the SparseFunction defines a multi-dimensional tensor. dtype : data-type, optional, default=np.float32 Any object that can be interpreted as a numpy data type. initializer : callable or any object exposing the buffer interface, optional Data initializer. If a callable is provided, data is allocated lazily. allocator : MemoryAllocator, optional Controller for memory allocation. To be used, for example, when one wants to take advantage of the memory hierarchy in a NUMA architecture. Refer to `default_allocator.__doc__` for more information. Notes ----- The parameters must always be given as keyword arguments, since SymPy uses `*args` to (re-)create the Dimension arguments of the symbolic object.
coeffs = np.ones((5, 2, 5))
s = PrecomputedSparseTimeFunction(name="s", grid=grid, npoint=npoint, nt=nt,
interpolation_coeffs=coeffs,
coordinates=coords, r=2)
pos = tuple(product((-grid.spacing[1], 0, grid.spacing[1],2*grid.spacing[1]),
(-grid.spacing[1], 0, grid.spacing[1],2*grid.spacing[1])))
interp_points = np.concatenate([base+p for p in pos])
#NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.scatter(coords[:, 0], coords[:, 1], s=10, label="Sparse positions")
ax.scatter(interp_points[:, 0], interp_points[:, 1], s=10, label="Interpolation support")
ax.grid(which = "major")
ax.grid(which = "minor", alpha = 0.2)
ax.xaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.yaxis.set_minor_locator(FixedLocator(np.linspace(0, 1, 51)))
ax.set_title("Off the grid sparse positions")
for i in range(npoint):
ax.annotate("(%.3f, %.3f)" % (coords[i, 0], coords[i, 1]), coords[i, :])
ax.legend()
plt.show()
f.data.fill(0)
op = Operator([Eq(f.forward, f+1)] + s.interpolate(f/16))
#NBVAL_IGNORE_OUTPUT
op()
s.data
Operator `Kernel` ran in 0.01 s
Data([[0., 0., 0., 0., 0.], [1., 1., 1., 1., 1.], [2., 2., 2., 2., 2.], [3., 3., 3., 3., 3.], [4., 4., 4., 4., 4.], [5., 5., 5., 5., 5.], [6., 6., 6., 6., 6.], [7., 7., 7., 7., 7.], [8., 8., 8., 8., 8.], [9., 9., 9., 9., 9.], [0., 0., 0., 0., 0.]], dtype=float32)
#NBVAL_IGNORE_OUTPUT
u.data.fill(0)
op = Operator(s.inject(u, expr=s))
op()
Operator `Kernel` ran in 0.01 s
PerformanceSummary([(PerfKey(name='section0', rank=None), PerfEntry(time=4.1e-05, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
#NBVAL_IGNORE_OUTPUT
plt.figure(figsize=(10, 10))
plt.imshow(u.data[1], vmin=0, vmax=2, cmap="jet", extent=[0,1,0,1])
plt.colorbar(fraction=0.046, pad=0.04)
plt.title("Precomputed weights")
plt.show()