Subgrids

A subgrid is a coarser grid whose points are a subset of a parent Grid’s points. SubGrid is the API for building one.

This tutorial covers

  1. Creating a SubGrid – the constructor and what it derives
  2. Defining Functions on it
  3. Writing Operators – moving data between the two grids
  4. Finite differences
  5. What SubGrid does not cover
  6. The generated code – how the compiler lowers it

It covers integer subsampling, where subgrid points coincide exactly with parent points, so no interpolation is involved.

Creating a SubGrid

SubGrid(parent, factor) keeps every factor-th point of parent along each Dimension, so subgrid point i sits at parent index i*factor.

import numpy as np
from devito import (CondEq, ConditionalDimension, Eq, Function, Grid, Operator,
                    SubGrid, TimeFunction)
from devito.ir.iet import FindNodes, Iteration

grid = Grid(shape=(12, 12), extent=(11., 11.))
x, y = grid.dimensions

sub = SubGrid(grid, factor=3)
sub
SubGrid[shape=(4, 4), factor=(3, 3)]

Everything else is derived from the parent. shape is ceil(n/factor) per axis:

print('parent  shape', grid.shape, ' subgrid shape', sub.shape)
parent  shape (12, 12)  subgrid shape (4, 4)

spacing is scaled by the factor, while extent and origin are inherited – the subgrid spans the same physical domain with fewer points.

print('spacing  parent', tuple(float(v) for v in grid.spacing),
      ' subgrid', tuple(float(v) for v in sub.spacing))
print('extent   parent', tuple(float(v) for v in grid.extent),
      ' subgrid', tuple(float(v) for v in sub.extent))
print('origin   parent', tuple(float(v) for v in grid.origin),
      ' subgrid', tuple(float(v) for v in sub.origin))
spacing  parent (1.0, 1.0)  subgrid (3.0, 3.0)
extent   parent (11.0, 11.0)  subgrid (9.0, 9.0)
origin   parent (0.0, 0.0)  subgrid (0.0, 0.0)

parent and factor are the two things SubGrid adds over Grid:

sub.parent is grid, sub.factor
(True, (3, 3))

Its dimensions are ConditionalDimensions derived from the parent’s. These are what carry the subsampling factor into the generated code.

sub.dimensions, [d.parent for d in sub.dimensions]
((xs, ys), [x, y])

Pass name to control how they are named – useful when building several SubGrids off one parent.

SubGrid(grid, factor=2, name='c').dimensions
(xc, yc)

A factor per Dimension subsamples each axis independently, and 1 leaves an axis untouched:

mixed = SubGrid(grid, factor=(3, 1))
print('shape     ', mixed.shape)
print('dimensions', mixed.dimensions)        # `y` passed through unchanged
shape      (4, 12)
dimensions (xs, y)

Being a Grid subclass, a SubGrid also carries the usual Grid attributes:

sub.dim, sub.dtype, sub.time_dim, sub.stepping_dim
(2, numpy.float32, time, t)

What is rejected

Only factors that keep subgrid points on parent points are allowed. A non-integer factor would place them in between, which needs interpolation rather than subsampling.

for bad in (2.5, 0, (2, 2, 2)):
    try:
        SubGrid(grid, factor=bad)
    except ValueError as e:
        print(f'factor={bad!r:10} -> {e}')
factor=2.5        -> Non-integer subsampling factor `2.5`. A SubGrid's points must coincide with its parent's, so the factor must be a whole number of parent cells -- use `2` or `3`. Subsampling onto points that fall between parent points requires interpolation, which SubGrid does not perform
factor=0          -> Subsampling factor must be >= 1, got `0`
factor=(2, 2, 2)  -> Expected 1 or 2 factors for a 2D parent Grid, got 3

Nor can a SubGrid be subsampled again; combine the factors off the original parent instead.

try:
    SubGrid(sub, factor=2)
except ValueError as e:
    print(e)
Cannot subsample a SubGrid. `2` applied to a SubGrid of factor `(3, 3)` is the same as factor `(6, 6)` applied to its parent Grid directly, so use `SubGrid(parent.parent, factor=(6, 6))`

Defining Functions on a SubGrid

This is the main convenience: a Function on a SubGrid needs neither shape nor dimensions, both are taken from the grid – exactly as for a normal Grid.

u = Function(name='u', grid=grid)                      # on the parent
usub = Function(name='usub', grid=sub, space_order=0)  # on the subgrid

print('u    ', u.shape, u.dimensions)
print('usub ', usub.shape, usub.dimensions)
u     (12, 12) (x, y)
usub  (4, 4) (xs, ys)

TimeFunction works too, and steps over the parent’s time Dimension – a SubGrid coarsens space only.

usave = TimeFunction(name='usave', grid=sub, space_order=0)
usave.dimensions, usave.dimensions[0] is grid.stepping_dim
((t, xs, ys), True)

Equivalent, by hand

SubGrid is a convenience over ConditionalDimension, which can still be used directly. The manual form is equivalent but requires naming each Dimension and computing the shape yourself.

xs_manual = ConditionalDimension(name='xsm', parent=x, factor=3)
ys_manual = ConditionalDimension(name='ysm', parent=y, factor=3)
f_manual = Function(name='f_manual', grid=grid,
                    dimensions=(xs_manual, ys_manual),
                    shape=(4, 4), space_order=0)        # shape by hand

f_manual.shape == usub.shape
True

Writing Operators

An Eq between a parent Function and a SubGrid Function is a strided copy: Devito matches the subgrid points to their parent counterparts.

# NBVAL_IGNORE_OUTPUT
u.data[:] = np.arange(144).reshape(12, 12)
usub.data[:] = 0.

Operator(Eq(usub, u)).apply()

assert np.array_equal(usub.data, u.data[::3, ::3])
usub.data
Operator `Kernel` ran in 0.01 s
Data([[  0.,   3.,   6.,   9.],
      [ 36.,  39.,  42.,  45.],
      [ 72.,  75.,  78.,  81.],
      [108., 111., 114., 117.]], dtype=float32)

Each row and column of usub comes from every third row and column of u.

The reverse direction works the same way – writing the parent from the subgrid scatters onto the subsampled points, leaving the others untouched.

# NBVAL_IGNORE_OUTPUT
u.data[:] = -1.
usub.data[:] = 7.

Operator(Eq(u, usub)).apply()

assert np.all(u.data[::3, ::3] == 7.)
assert np.all(u.data[1::3, :] == -1.)      # untouched rows
u.data[:4, :4]
Operator `Kernel` ran in 0.01 s
Data([[ 7., -1., -1.,  7.],
      [-1., -1., -1., -1.],
      [-1., -1., -1., -1.],
      [ 7., -1., -1.,  7.]], dtype=float32)

Restricting the iteration range maps onto the subgrid. With x_m=2 and factor=3, the first subgrid point at or after x=2 is x=3, i.e. subgrid index 1, so index 0 is never written.

# NBVAL_IGNORE_OUTPUT
u.data[:] = 1.
usub.data[:] = -1.                     # sentinel: never written

Operator(Eq(usub, u)).apply(x_m=2)

usub.data[:, 0]
Operator `Kernel` ran in 0.01 s
Data([-1.,  1.,  1.,  1.], dtype=float32)

Several SubGrids can appear in one Operator, each with its own factor.

# NBVAL_IGNORE_OUTPUT
sub_a = SubGrid(grid, factor=2, name='a')
sub_b = SubGrid(grid, factor=4, name='b')
fa = Function(name='fa', grid=sub_a, space_order=0)
fb = Function(name='fb', grid=sub_b, space_order=0)

u.data[:] = np.arange(144).reshape(12, 12)
Operator([Eq(fa, u), Eq(fb, u)]).apply()

assert np.array_equal(fa.data, u.data[::2, ::2])
assert np.array_equal(fb.data, u.data[::4, ::4])
fa.shape, fb.shape
Operator `Kernel` ran in 0.01 s
((6, 6), (3, 3))

Finite differences on a SubGrid

Derivatives of a SubGrid Function are taken with respect to the coarse spacing, since the Dimension’s spacing is factor*h.

sub.dimensions[0].spacing, x.spacing
(3*h_x, h_x)

So the stencil steps factor*h_x and the weights scale accordingly:

fsub = Function(name='fsub', grid=sub, space_order=2)
fsub.dx.evaluate

\(\displaystyle - \frac{0.333333333 fsub(xs, ys)}{h_{x}} + \frac{0.333333333 fsub(xs + 3*h_x, ys)}{h_{x}}\)

Checking it numerically on a linear field, whose derivative is exact:

# NBVAL_IGNORE_OUTPUT
grid_l = Grid(shape=(21, 21), extent=(20., 20.))
sub_l = SubGrid(grid_l, factor=2)

fl = Function(name='fl', grid=sub_l, space_order=2)
dfl = Function(name='dfl', grid=sub_l, space_order=2)

# fl = 3*x, sampled on the coarse grid
for i in range(fl.shape[0]):
    fl.data[i, :] = 3.*i*float(sub_l.spacing[0])

Operator(Eq(dfl, fl.dx)).apply()

assert np.allclose(dfl.data[1:-1, 1:-1], 3., atol=1e-4)
float(dfl.data[5, 5])
Operator `Kernel` ran in 0.01 s
3.0

What SubGrid does not cover

Two related things are deliberately outside SubGrid and still want a ConditionalDimension built by hand.

Subsampling in time. A SubGrid coarsens space, so saving a wavefield every factor timesteps uses a ConditionalDimension on the time Dimension.

# NBVAL_IGNORE_OUTPUT
nt, factor = 9, 4
grid_t = Grid(shape=(4, 4))
ts = ConditionalDimension(name='ts', parent=grid_t.time_dim, factor=factor)

v = TimeFunction(name='v', grid=grid_t)
vsave = TimeFunction(name='vsave', grid=grid_t, time_dim=ts,
                     save=(nt + factor - 1)//factor)

Operator([Eq(v.forward, v + 1.), Eq(vsave, v)]).apply(time_M=nt-2)

[float(vsave.data[i, 0, 0]) for i in range(vsave.shape[0])]
Operator `Kernel` ran in 0.01 s
[0.0, 4.0, 0.0]

Extra conditions. SubGrid takes a factor, not an arbitrary predicate; for factor and a condition, construct the ConditionalDimension directly.

xs_cond = ConditionalDimension(name='xsc', parent=x, factor=3,
                               condition=CondEq(x % 2, 0))
xs_cond.factor, xs_cond.condition
(3, Eq(Mod(x, 2), 0))

Generic subsampling – subgrid points that fall between parent points, requiring interpolation – is a separate topic and is not covered here.

The generated code

The rest of this tutorial shows how Devito lowers the above, which is useful when reading or profiling the generated C.

Devito iterates the subgrid directly: one unit-stride loop per subsampled Dimension, with parent accesses rescaled – rather than walking every parent point under an if (x%factor == 0) guard.

def show(op):
    """Print the Operator's loop nest."""
    print(op.body.body[-1])


def steps(op):
    """The loop increment of each Iteration, by Dimension name."""
    return {i.dim.name: i.limits[2] for i in FindNodes(Iteration).visit(op)}

One loop per subsampled Dimension

# NBVAL_IGNORE_OUTPUT
op = Operator(Eq(usub, u))
show(op)
START(section0)
for (int xs = xs_m; xs <= xs_M; xs += 1)
{
  #pragma omp simd aligned(u,usub:64)
  for (int ys = ys_m; ys <= ys_M; ys += 1)
  {
    usub[xs][ys] = u[xs*xsf + 1][ys*ysf + 1];
  }
}
STOP(section0,timers)

Note what the loop iterates: xs, the subsampled Dimension itself, with a unit increment, over its own bounds xs_m .. xs_M. There is no if in the nest.

  • the subgrid access is simply usub[xs];
  • the parent access is rescaled to u[xs*xsf + 1].

The unit stride is deliberate. A x += factor loop would also remove the guard, but a non-unitary increment escapes the Lamport-style dependence analysis the compiler relies on, so it is avoided.

The bounds need no adjusting either: xs_m and xs_M are ConditionalDimension’s own runtime arguments, derived as ceil(x_m/factor) and floor(x_M/factor).

steps(op)
{'xs': 1, 'ys': 1}

The innermost Dimension

Subsampling the innermost (vectorized) Dimension is handled the same way. The alternative – a guard inside the SIMD loop – would mean a predicated store on every iteration.

# NBVAL_IGNORE_OUTPUT
sub_y = SubGrid(grid, factor=(1, 2))      # coarsen y only
usub_y = Function(name='usub_y', grid=sub_y, space_order=0)

op_inner = Operator(Eq(usub_y, u))
show(op_inner)
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
  #pragma omp simd aligned(u,usub_y:64)
  for (int ys = ys_m; ys <= ys_M; ys += 1)
  {
    usub_y[x][ys] = u[x + 1][ys*ysf + 1];
  }
}
STOP(section0,timers)
steps(op_inner)
{'x': 1, 'ys': 1}

The y loop is strided and the #pragma omp simd is retained.

# NBVAL_IGNORE_OUTPUT
u.data[:] = np.arange(144).reshape(12, 12)
Operator(Eq(usub_y, u)).apply()

assert np.array_equal(usub_y.data, u.data[:, ::2])
usub_y.data[:3]
Operator `Kernel` ran in 0.01 s
Data([[ 0.,  2.,  4.,  6.,  8., 10.],
      [12., 14., 16., 18., 20., 22.],
      [24., 26., 28., 30., 32., 34.]], dtype=float32)

Several Dimensions at once

Each subsampled Dimension gets its own loop, independently.

# NBVAL_IGNORE_OUTPUT
sub_xy = SubGrid(grid, factor=(3, 2))
usub_xy = Function(name='usub_xy', grid=sub_xy, space_order=0)

op_both = Operator(Eq(usub_xy, u))
show(op_both)
START(section0)
for (int xs = xs_m; xs <= xs_M; xs += 1)
{
  #pragma omp simd aligned(u,usub_xy:64)
  for (int ys = ys_m; ys <= ys_M; ys += 1)
  {
    usub_xy[xs][ys] = u[xs*xsf + 1][ys*ysf + 1];
  }
}
STOP(section0,timers)
# NBVAL_IGNORE_OUTPUT
u.data[:] = np.arange(144).reshape(12, 12)
Operator(Eq(usub_xy, u)).apply()

assert np.array_equal(usub_xy.data, u.data[::3, ::2])
steps(op_both)
Operator `Kernel` ran in 0.01 s
{'xs': 1, 'ys': 1}

When a guard is used instead

Striding is applied only where it provably preserves semantics. Where it does not, Devito visits every point under an if (x%factor == 0) guard. Results are the same either way; only the generated loop differs. Four cases, each with its code.

1. The loop nest is shared with an unguarded equation

Eq(f, 1) must run at every x, while Eq(g, f) runs every fourth. Both land in one loop nest, so striding it would also stride Eq(f, 1). The guard is kept.

# NBVAL_IGNORE_OUTPUT
grid1 = Grid(shape=(12,))
sub1 = SubGrid(grid1, factor=4)

f = Function(name='f', grid=grid1)
g = Function(name='g', grid=sub1, space_order=0)

op_shared = Operator([Eq(f, 1), Eq(g, f)])
show(op_shared)
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
  f[x + 1] = 1;
  if ((x)%(xsf) == 0)
  {
    g[x / xsf] = f[x + 1];
  }
}
STOP(section0,timers)
# NBVAL_IGNORE_OUTPUT
f.data[:] = 0.
g.data[:] = -1.
op_shared.apply()

# `f` written at every point, `g` only on the subgrid
assert np.all(f.data == 1.)
assert np.all(g.data == 1.)
f.data[:3], g.data[:]
Operator `Kernel` ran in 0.01 s
(Data([1., 1., 1.], dtype=float32), Data([1., 1., 1.], dtype=float32))

f[x + 1] = 1 sits outside the if, g[x / xs1f] inside it – which is exactly what is required.

2. The subsampled Dimension is sequential

Striding skips points, which is only sound when no dependence is carried along the Dimension. Here g reads its own previous subgrid value, so x is sequential.

# NBVAL_IGNORE_OUTPUT
grid2 = Grid(shape=(12, 12))
sub2 = SubGrid(grid2, factor=(3, 1))
xs3 = sub2.dimensions[0]

f2 = Function(name='f2', grid=grid2)
g2 = Function(name='g2', grid=sub2, space_order=0)

op_seq = Operator(Eq(g2, g2.subs(xs3, xs3-1) + f2))
show(op_seq)
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
  if ((x)%(xsf) == 0)
  {
    for (int y = y_m; y <= y_M; y += 1)
    {
      g2[x / xsf][y] = f2[x + 1][y + 1] + g2[(x / xsf) - 1][y];
    }
  }
}
STOP(section0,timers)
# NBVAL_IGNORE_OUTPUT
f2.data[:] = 1.
g2.data[:] = 0.
op_seq.apply()

# Each subgrid row accumulates the one before it; row 0 reads the halo (zero)
assert np.allclose(g2.data[:, 0], [0., 1., 2., 3.])
g2.data[:, 0]
Shrinking bounds (`x_m`, `x_M`); some `x` points will not be computed. Likely insufficient space_order for the derivatives.
/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:718: DevitoWarning: Shrinking bounds (`x_m`, `x_M`); some `x` points will not be computed. Likely insufficient space_order for the derivatives.
  warn(f"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); "
Shrinking bounds (`xs_m`, `xs_M`); some `xs` points will not be computed. Likely insufficient space_order for the derivatives.
/home/mloubout/DevitoCodes/devito/devito/operator/operator.py:718: DevitoWarning: Shrinking bounds (`xs_m`, `xs_M`); some `xs` points will not be computed. Likely insufficient space_order for the derivatives.
  warn(f"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); "
Operator `Kernel` ran in 0.01 s
Data([0., 1., 2., 3.], dtype=float32)

Time subsampling is the everyday instance of this: the time loop is sequential, so a saved-wavefield Eq keeps its guard. Note this one uses a ConditionalDimension directly – a SubGrid coarsens space, so it has no say over the time axis.

# NBVAL_IGNORE_OUTPUT
nt, factor = 9, 4
grid3 = Grid(shape=(4, 4))
ts = ConditionalDimension(name='ts', parent=grid3.time_dim, factor=factor)

v = TimeFunction(name='v', grid=grid3)
vsave = TimeFunction(name='vsave', grid=grid3, time_dim=ts,
                     save=(nt + factor - 1)//factor)

op_time = Operator([Eq(v.forward, v + 1.), Eq(vsave, v)])
show(op_time)
for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
{
  START(section0)
  for (int x = x_m; x <= x_M; x += 1)
  {
    #pragma omp simd aligned(v:64)
    for (int y = y_m; y <= y_M; y += 1)
    {
      v[t1][x + 1][y + 1] = v[t0][x + 1][y + 1] + 1.0F;
    }
  }
  STOP(section0,timers)

  if ((time)%(tsf) == 0)
  {
    START(section1)
    for (int x = x_m; x <= x_M; x += 1)
    {
      #pragma omp simd aligned(v,vsave:64)
      for (int y = y_m; y <= y_M; y += 1)
      {
        vsave[time / tsf][x + 1][y + 1] = v[t0][x + 1][y + 1];
      }
    }
    STOP(section1,timers)
  }
}

The saved slices are correct – time subsampling behaves exactly as it always has, it simply does not take the strided form.

# NBVAL_IGNORE_OUTPUT
op_time.apply(time_M=nt-2)

nwritten = (nt - 2)//factor + 1
for i in range(nwritten):
    assert np.allclose(vsave.data[i], i*factor)

[float(vsave.data[i, 0, 0]) for i in range(vsave.shape[0])]
Operator `Kernel` ran in 0.01 s
[0.0, 4.0, 0.0]

3. A condition is attached as well

A ConditionalDimension may carry an arbitrary condition alongside a factor. The combined test is no longer “every factor-th point”, so it is evaluated pointwise.

# NBVAL_IGNORE_OUTPUT
# `SubGrid` does not expose `condition`, so this one is built by hand
xs4 = ConditionalDimension(name='xs4', parent=x, factor=3,
                           condition=CondEq(x % 2, 0))
w = Function(name='w', grid=grid, dimensions=(xs4, y), shape=(4, 12),
             space_order=0)

u.data[:] = 1.
w.data[:] = -1.
op_cond = Operator(Eq(w, u))
op_cond.apply()
show(op_cond)
Operator `Kernel` ran in 0.01 s
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
  if ((x & 1) == 0 && (x)%(xs4f) == 0)
  {
    #pragma omp simd aligned(u,w:64)
    for (int y = y_m; y <= y_M; y += 1)
    {
      w[x / xs4f][y] = u[x + 1][y + 1];
    }
  }
}
STOP(section0,timers)

4. Accesses are not affine in the subsampled Dimension

An indirectly indexed read (u[idx[x, y], y]) is not affine in x, so the points a stride would skip cannot be shown to be unneeded.

# NBVAL_IGNORE_OUTPUT
idx = Function(name='idx', grid=grid, dtype=np.int32)
gi = Function(name='gi', grid=sub, space_order=0)

idx.data[:] = 0                      # read row 0 of `u` everywhere
u.data[:] = np.arange(144).reshape(12, 12)
op_indirect = Operator(Eq(gi, u[idx[x, y], y]))
op_indirect.apply()

assert np.array_equal(gi.data, u.data[0, ::3] * np.ones((4, 1)))
show(op_indirect)
Operator `Kernel` ran in 0.01 s
START(section0)
for (int x = x_m; x <= x_M; x += 1)
{
  if ((x)%(xsf) == 0)
  {
    #pragma omp simd aligned(gi,idx,u:64)
    for (int y = y_m; y <= y_M; y += 1)
    {
      if ((y)%(ysf) == 0)
      {
        int r0 = idx[x + 1][y + 1];
        gi[x / xsf][y / ysf] = u[r0 + 1][y + 1];
      }
    }
  }
}
STOP(section0,timers)

Summary

SubGrid(parent, factor) builds a coarser grid from a parent Grid, keeping every factor-th point. It is a Grid subclass, so Functions, TimeFunctions, Operators and finite differences all work as usual:

sub  = SubGrid(grid, factor=2)            # or factor=(2, 3); 1 = passthrough
usub = Function(name='usub', grid=sub)    # shape and dimensions derived
Operator(Eq(usub, u)).apply()             # strided copy from the parent
derived from the parent
shape ceil(n/factor) per axis
spacing factor * parent.spacing
extent, origin, dtype, time_dim inherited unchanged
dimensions ConditionalDimensions on the parent’s
parent, factor the two attributes SubGrid adds

Use a ConditionalDimension directly for subsampling in time, or to combine a factor with an extra condition – SubGrid covers neither.

Under the hood Devito iterates the subgrid with a unit stride and rescales parent accesses, falling back to a pointwise if (x%factor == 0) guard in the four cases above. Both compute the same thing; only the generated loop differs.

Generic subsampling – subgrid points that fall between parent grid points and so require interpolation – is a separate topic, not covered here.

Back to top