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

The ten figures that follow all turn on a single pattern: a residual written symbolically, rearranged into an explicit update, then handed to the compiler. The three propagators benchmarked in Fig. 09 are three instances of it, and the stencil radius that drives every legality argument on this page is set here, as space_order.

Nothing on this page depends on the acoustic equation specifically. The tiling argument needs only that the update is explicit, local and affine — properties the compiler reads off the symbolic form rather than off any particular PDE. The same pipeline serves TTI and elastic in Fig. 09 without modification, and the precomputation of Fig. 06 observes the injection operator instead of modelling it, for the same reason.

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 pays off 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 strip-mine and interchange, sliding time inward:

Sliding time inside the block loop — the naive interchange shown on the right — is not legal: a block's right-hand edge reads points in the next block, which runs after it and has not reached that timestep yet. The families above are the repairs, each restoring legality a different way — by leaning the tiles, by recomputing the shared margin, or by cutting cone-shaped tiles whose boundaries follow the dependence cone, so every tile can start at once — concurrent start, in the diamond panel below. This paper uses skewed wavefront tiling, which Fig. 03 covers. The problem it addresses is orthogonal to the repair you choose: none of these families can be applied unmodified once the simulation contains sources that are not on the grid: injection is additive, so tiles that overlap or lean across a seam would need ownership masks to avoid double-counting a contribution — the machinery Fig. 06 builds.

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.

Both panels share one clock: 2r + 1 neighbours at t feed one point at t+1 on the left, and the same radius grows a 3D star on the right — 6r + 1 points along its six arms. Every blocking scheme below has to respect those arrows. Watch the two grow together as the order cycles: 5 and 13 points at space order 4, 9 and 25 at order 8, 13 and 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 full waveform inversion (FWI) and reverse time migration (RTM) — the two seismic-imaging workloads this solver is built for — 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. They are not symmetric as schedules, though: injection writes, so it pins the order, while interpolation only reads and can be deferred — as far as the time buffer allows.

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. The same geometry explains the paper's numbers: 1.6× at space order 4, about 1.13× at order 8, down to a token 5–10% or none at all by order 12, depending on the kernel (Fig. 09).

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 the polyhedral frameworks — 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 algebraically unchanged — the weights are pre-applied, so results may differ in the last bits — 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 interpolation hoisted into the precomputation, and a single iteration space that Fig. 03's skew can now be applied to. Legality rests on three conditions the paper states: coordinates fixed across the time domain, the mask and ID structures read-only once built, and the injection ordered after the stencil at the same timestep — which is what Inc preserves. This is the classic inspector/executor split from the sparse-compilation literature: Band A is the inspector, a cheap pass that resolves the data-dependent structure once; the fused nest is the executor, which runs many times against what the inspector found. 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. In the published release the operator is then built as op2 = Operator([eqxb, eqyb, eqxb2, eqyb2, stencil_2, eq0, eq1, eq2], opt=('advanced')), where each eqx… copies a runtime scalar out of a block_sizes Function — the tile and block shapes become runtime arguments, so Fig. 09's tuning sweeps need no C edits and no recompiles. That line is the entire compiler-facing surface of the technique.

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. The block sizes in the drawing are not constants: in the released code they arrive as scalars read from a block_sizes Function, so the autotuner re-shapes the nest at runtime on a single compiled binary. Note the order: t sits outside the block loops, so every block of a tile advances in lockstep at each step — the wave order is between tiles, which is both why the OpenMP parallelism on blocks is legal and where the pipeline fill and drain come from.

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: this technique adds no redundant computation (Fig. 08), so GFLOP/s moves with the same speedup shown here — but flop count isn't fixed across every tiling family (overlapped tiling trades redundant flops for concurrency), so points/s stays the metric of record, and the roofline is read as a diagnostic rather than a scoreboard. Temporal blocking is often reported at larger multiples, but against baselines tuned to a lesser degree and without the sparse-operator constraint that had to be lifted here; tuned against tuned is the conservative reading.

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 left unblocked on purpose, 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 the reading offered here 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 page does not claim, and where the work goes next.