Temporal blocking of finite-difference stencil operators with sparse “off-the-grid” sources

Temporal blocking reorders the loop nest so that a cached block is advanced through several timesteps before the next block is touched — one trip to memory serving many updates instead of one. It is worth doing because the field dwarfs the cache: 512³ in single precision is 512 MiB against roughly 48 MiB of last-level cache, with hundreds of timesteps to get through. The transformation is a loop interchange, sliding time inward:

That interchange is not legal as written: a block's edge needs neighbours from the block beside it, which has not reached the timestep yet. Each family above restores legality differently — by leaning, by recomputing, or by choosing a tile shape that never asks. This paper uses skewed wavefront tiling, which Fig. 03 covers. The problem this paper addresses is orthogonal to that choice: none of these families can be applied at all once the simulation contains sources that are not on the grid.

Discretise a PDE with explicit finite differences and every grid point at timestep t+1 becomes a weighted sum of its neighbours at t. The reach of that sum is the space order. It is dense, regular and predictable, and therefore tileable. The dependence pattern below is what every legality argument on this page refers back to.

Left, the canonical case: three neighbours at t feed one point at t+1, and the update marches along the row. Every blocking scheme below has to respect those three arrows. Right, the same idea in 3D — a star reaching r = so/2 points along each axis, so the point count is 6r + 1. Watch the arms as the order cycles: 13 points at space order 4, 25 at order 8, 37 at order 12. The paper benchmarks all three, and Fig. 10 shows what that growth costs.

A seismic or medical-imaging simulation is not a clean stencil benchmark. Energy is injected by sources and sampled by receivers whose physical coordinates are not aligned with the discretisation. Both sit on the same grid at the same time — a shot record has one source and thousands of receivers — and each spreads its contribution across the surrounding points through weights that vary continuously with position. Drag either marker and watch the weights move.

Two iteration spaces share one field. The stencil walks the lattice in order; these two arrive through indirect accesses at coordinates that were never part of it. The arrows run opposite ways on purpose — injection scatters outward, interpolation gathers inward — because the receiver operator is the adjoint of the source operator, built from the same weights.

That symmetry matters later: in FWI and RTM both appear in the same operator, so anything that fixes one has to fix the other. In time the two are aligned, acting on the same timestep index; in space neither is.

Spatial blocking cuts the domain into cache-sized tiles, then throws the data away at the end of every timestep. Temporal blocking keeps going — advance one tile through several timesteps before touching the next. That is only legal if everything a point reads has already been computed. Tiles run left to right, so the dangerous dependence is the one reaching right, into a tile that has not started. Skewing leans the tiles to pull it back.

Every dot is a grid point at one timestep; every parallelogram is one space–time tile, numbered in execution order. Press play and the sweep runs the way the generated code does: tile ① is carried through all four timesteps before tile ② begins, so its data stays in cache for the whole time block. The white marker is the point being updated, and the arrows are the 2r+1 values it reads from the row below. At skew 0, watch the right-hand arrows go red every time the sweep reaches a tile's trailing edge: it is asking for values from a tile that has not run. Raise the skew to r and the counter stays at zero for the whole run.

Raise it beyond r and the run is still clean, but the tile now sweeps across more columns than it updates for no extra reuse. The minimum legal skew is r, so a wider stencil forces a steeper lean whether or not you want one. This is why the paper's 1.6× at space order 4 falls to about 1.13× at order 8 and vanishes by order 12.

Skewing settles the stencil's own dependences. It does nothing for the source, because the source is a separate loop over a separate iteration space — one pass per timestep, over all sources, across the whole domain. That was harmless without time-tiling, when the whole domain advanced in lockstep. Once tiles run ahead of each other, a single source whose footprint straddles a tile seam has to be applied at two different moments in the schedule. One loop cannot be in two places. It is also why a general-purpose tiler cannot rescue this: the injection loop’s accesses are non-affine, so PLUTO, Polly, Loo.py and CLooG tile the stencil beside it and stop at this one.

Left, one source scatters into four grid points, and the tile seam runs straight through them. Right, the same situation in execution order: tile A is advanced through its whole time block before tile B starts, so the single per-timestep injection loop would have to fire once inside A's pass and once inside B's pass. This is the applicability barrier the paper identifies. Nothing about wave propagation forbids temporal blocking; the sparse operator simply cannot be scheduled until it is folded into the grid's own iteration space.

Fig. 03 showed skewing repairing the stencil's own dependences. Fig. 04 showed the source failing for a different reason. Put them on the same axes and the asymmetry becomes visible: the grid and the source are fixed in physical space, but the tiles lean, so the seams sweep leftward as time advances. A source sitting at one position gets crossed. Turning up the skew does not remove the crossing — it just moves which timestep it happens at.

Grid points sit in fixed columns; the source's two affected points are marked at every timestep. Tile seams lean left as t grows, since that lean is the skew. As the sweep climbs, a seam eventually slides across the source. On the timesteps where it does, the two points that must receive the same contribution at the same simulated time belong to tiles running at different moments.

Watch the two lamps as you drag: the left one turns green at s ≥ 2 and stays green; the right one stays red no matter what you do to the skew, because the injection loop's placement is not a geometry problem. Only switching the injection to precomputed turns it green — which is Fig. 06.

The fix removes the second iteration space entirely. Inject into an empty field for a single timestep (a few more only if the wavelet starts at zero) to find out exactly which grid points a source can ever touch; record those points as on-grid arrays; then rewrite the injection as a masked term evaluated inside the stencil loop. The arithmetic is unchanged, but it is now aligned, vectorisable and legal to time-tile, and the compiler generates it from a symbolic description rather than the user writing it by hand.

Band A pays a one-off cost — a single timestep of injection into an empty field — to learn the footprint; the paper reports this as negligible next to the gain, and it holds for non-linear injection too, because it observes the operator rather than modelling it. Band B stores that footprint three ways: a binary mask, an identifier into the pre-weighted source values, and a per-column count of non-zeros so the sweep can skip empty z slices instead of testing every point.

Band C is the result: one loop nest, no indirect scatter, the same operation count, and a single iteration space that Fig. 03's skew can now be applied to. Fig. 07 shows the Devito that produces it, line by line.

This is what the transformation costs the person writing the solver. On the left, the textbook Devito acoustic propagator — five lines, and the source injection is a single method call. On the right, the same model rewritten so that the injection is expressible inside the grid's own loop nest. Click any highlighted line to see what it does and why it is there.

The physics is untouched — pde and pde_2 are the same equation. Everything added on the right exists to move one thing: the source's contribution, out of its own loop and into the grid's. Two mechanisms do the work. implicit_dims forces an expression into a loop nest it does not syntactically belong to; symbolic_max lets a loop bound be a value read from an array at runtime. Together they let a sparse, data-dependent operation be written as a dense one.

The loop nest after both transformations. Reading outwards: two innermost loops share one x, y pair, which is the fusion from Fig. 06; the skew appears as an index shift, which is the lean from Fig. 03; then blocks, tiles, and the time-tile loop. Note where the parallelism sits, and what is not blocked.

Two loops share the innermost x, y position: the dense stencil sweep over z, and a short injection loop whose length is read from nnz_mask. That adjacency turns the source's data dependence into a local, affine one that the outer tiling can reason about. It also explains why the sparsity reduction in Fig. 06 matters so much: without nnz_mask, the second loop would run the full nz and multiply by zero almost every time.

Three propagators of industrial significance — isotropic acoustic, anisotropic acoustic (TTI) and isotropic elastic — on 512³ grids over 512 ms, benchmarked on Azure Broadwell and Skylake VMs against Devito's own aggressively auto-tuned, vectorised, spatially blocked kernels. Both sides were swept over the full tile and block parameter space, so this is a tuned-versus-tuned comparison.

Up to 1.6× at space order 4 and around 1.13× or better at order 8, varying by propagator, on kernels that were already flop-optimised and tuned and that previously could not be time-tiled at all. The metric is GPoints/s, time to solution, rather than GFLOP/s: a temporally blocked kernel can post a lower flop rate and still finish sooner, which is why the roofline is read as a diagnostic rather than a scoreboard.

Fig. 09 shows the speedup collapsing between space order 4 and 12. The arithmetic behind it also accounts for the tile sizes in Table I. Two constraints pull the tile width in opposite directions. Skewing wants the tile wide, because a leaning tile touches W + r·T columns while only updating W of them. Cache wants the tile narrow, because z is never blocked, so a tile holds W² · nz elements for every field and every time buffer at once. The stencil radius sets how sharp the conflict is.

The cyan curve is skew efficiency: it climbs with tile width because the r·T columns of lean amortise over a wider tile. The violet curve is the working set: it grows as W² and crosses the L3 ceiling fast, because the z dimension is deliberately left unblocked for vectorisation. Between them sits a window of workable tile widths.

At space order 4 the window is comfortable. At order 12 the radius triples, efficiency needs a tile far wider than cache will hold, and the window shuts — which is our reading of why Table I's order-12 rows jump to 256×256: the autotuner has effectively backed away from temporal blocking rather than found a better tile. Note also that elastic, with nine state parameters, hits the ceiling at a much narrower tile than acoustic does, which is consistent with its tiles staying at 32×32 while acoustic moves to 64×64.

What this artifact deliberately does not claim, and where the work goes next.