Staggered grids and the interp-mode option

When fields in an equation live on different sub-grid offsets - a so-called staggered grid - Devito automatically interpolates the right-hand side onto the location of the left-hand side when it discretises the Eq. Most of the time this is transparent: you write the equation as usual, build an Operator, and Devito picks reasonable interpolations.

For some applications (especially elastic wave propagation and other vector PDEs) the way products of staggered fields are interpolated matters. Devito exposes this choice as a symbolic option on the Operator, passed via sym_opt:

Operator(eq, sym_opt={'interp-mode': 'direct'})     # default
Operator(eq, sym_opt={'interp-mode': 'symmetric'})  # alternative

sym_opt is the place where mathematical choices made during expression lowering live, separate from opt, which controls code generation and performance. interp-mode accepts two values: 'direct' (the default) and 'symmetric'. Both modes produce the same answer when at least one factor already matches the target staggering - the choice only matters when every factor sits elsewhere.

This tutorial walks through:

You will not need to call any interpolation routine directly - we work entirely with Function, Eq, and Operator. For a tour of how _eval_at actually rewrites derivatives and functions under the hood, see 09_fd_evaluation.ipynb.

Staggered locations

A staggered grid stores different fields at different sub-grid offsets of a single cell. In 1D with spacing \(h\) the two canonical locations are NODE (\(x_i\)) and the half-shifted location \(x_{i + 1/2} = x_i + h/2\). A first derivative naturally lives between them: the centred two-point stencil

\[ u'\!\left(x_{i + 1/2}\right) \;\approx\; \frac{u(x_{i+1}) - u(x_i)}{h} \]

is second-order accurate at \(x_{i+1/2}\) but only first-order at \(x_i\). In higher dimensions any subset of axes can be staggered, giving up to \(2^n\) locations in \(n\) dimensions.

In Devito, the staggered keyword controls the offset. The usual 2D cases:

from devito import Function, Grid, NODE

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

# cell centre, x-face, y-face, corner
fn = Function(name='fn', grid=grid, space_order=4, staggered=NODE)
fx = Function(name='fx', grid=grid, space_order=4, staggered=x)
fy = Function(name='fy', grid=grid, space_order=4, staggered=y)
fxy = Function(name='fxy', grid=grid, space_order=4, staggered=(x, y))

Each Function knows its natural sample location, which we can read off the indices_ref property:

for f in (fn, fx, fy, fxy):
    print(f"{f.name:<3s} indices_ref = {tuple(f.indices_ref)}")
fn  indices_ref = (x, y)
fx  indices_ref = (x + h_x/2, y)
fy  indices_ref = (x, y + h_y/2)
fxy indices_ref = (x + h_x/2, y + h_y/2)

Mixing staggerings in an equation

Suppose we want \(f_\mathrm{xy} = f_\mathrm{n} \cdot \partial_x f_\mathrm{x}\). The three pieces live at three different locations: fn at the cell centre, fx.dx on the \(x\)-face (a derivative inherits its operand’s grid), and fxy at the corner. We just write it down:

from devito import Eq

eq = Eq(fxy, fn * fx.dx)
eq

\(\displaystyle fxy(x + h_x/2, y + h_y/2) = fn(x, y) \frac{\partial}{\partial x} fx(x + h_x/2, y)\)

So far this is purely symbolic. The interpolation happens when we build an Operator. To see the discretised update statement, we can walk the Operator’s IR and pick out its Expression nodes - the symbolic form is much easier to read than the generated C code:

from devito import Operator
from devito.ir.iet import Expression, FindNodes


def updates(op):
    # Symbolic update equations carried by ``op``
    return [n.expr for n in FindNodes(Expression).visit(op)]


op = Operator(eq)
updates(op)
[Eq(r0, 1/h_x),
 Eq(fxy[x + 4, y + 4], (0.0416666667*(r0*(fx[x + 2, y + 4] + fx[x + 2, y + 5]) - r0*(fx[x + 6, y + 4] + fx[x + 6, y + 5])) + 0.333333333*(-r0*(fx[x + 3, y + 4] + fx[x + 3, y + 5]) + r0*(fx[x + 5, y + 4] + fx[x + 5, y + 5])))*(0.25*fn[x + 4, y + 4] + 0.25*fn[x + 5, y + 4] + 0.25*fn[x + 4, y + 5] + 0.25*fn[x + 5, y + 5]))]

Reading the right-hand side: fn enters through a four-corner average (the 2D linear interpolation from cell centres to corners), while fx.dx is computed by a centred fourth-order stencil already at the corner location, with a 2-point \(y\)-shift baked in. Each factor has been moved to the corner independently - this is what we mean by “direct” interpolation.

(r0 is a hoisted reciprocal of h_x, which is why it appears in the stencil.)

The interp-mode option

For a product whose factors all live at a staggering different from the target, Devito offers two ways of building the discretisation. Because this is a mathematical choice (it changes the discrete operator, not the generated C code) it sits in the sym_opt dict on the Operator, not in opt:

Operator(eq, sym_opt={'interp-mode': 'direct'})     # default
Operator(eq, sym_opt={'interp-mode': 'symmetric'})  # alternative

'direct' is the default because it produces the smallest stencil. For operators that need to remain self-adjoint under discretisation, however, 'symmetric' is the right choice.

To see why, let us recall that each shift between two staggered grids is a linear map and therefore has a transpose. For the 1D two-point average,

\[ \mathbf{I}\,u\,(x_{i+1/2}) = \tfrac{1}{2}\bigl[u(x_i) + u(x_{i+1})\bigr], \]

\(\mathbf{I}\) maps cell centres to half-steps, and \(\mathbf{I}^{\!\top}\) goes the other way. Many PDE operators decompose into block matrices of the form \(\mathbf{I}\,\mathbf{A}\,\mathbf{I}^{\!\top}\). Preserving that structure in the discretisation makes the discrete operator self-adjoint whenever \(\mathbf{A}\) is. The two modes correspond to two ways of building such a product:

  • 'direct' interpolates every factor to the target independently. The result is a product of averages, easy to read and cheap to evaluate.
  • 'symmetric' first gathers every factor at a common block location (NODE if any factor sits there), forms the product there, then projects the result to the target with a single \(\mathbf{I}\). This preserves the \(\mathbf{I}\,\mathbf{A}\,\mathbf{I}^{\!\top}\) matrix structure.

The two agree whenever at least one factor already matches the target - in that case there is no “block” to gather to and no choice to make.

Here is the same equation built in both modes side by side:

op_dir = Operator(eq, sym_opt={'interp-mode': 'direct'})
op_sym = Operator(eq, sym_opt={'interp-mode': 'symmetric'})

print('direct   :')
for u in updates(op_dir):
    print('  ', u)
print()
print('symmetric:')
for u in updates(op_sym):
    print('  ', u)
direct   :
   Eq(r0, 1/h_x)
   Eq(fxy[x + 4, y + 4], (0.0416666667*(r0*(fx[x + 2, y + 4] + fx[x + 2, y + 5]) - r0*(fx[x + 6, y + 4] + fx[x + 6, y + 5])) + 0.333333333*(-r0*(fx[x + 3, y + 4] + fx[x + 3, y + 5]) + r0*(fx[x + 5, y + 4] + fx[x + 5, y + 5])))*(0.25*fn[x + 4, y + 4] + 0.25*fn[x + 5, y + 4] + 0.25*fn[x + 4, y + 5] + 0.25*fn[x + 5, y + 5]))

symmetric:
   Eq(r0, 1/h_x)
   Eq(r1[y + 1], r0*(0.5*((0.0416666667*(fx[x + 1, y + 5] - fx[x + 6, y + 5]) + 0.291666667*(-fx[x + 2, y + 5] + fx[x + 5, y + 5]) + 0.333333333*(-fx[x + 3, y + 5] + fx[x + 4, y + 5]))*fn[x + 4, y + 5] + (0.0416666667*(fx[x + 2, y + 5] - fx[x + 7, y + 5]) + 0.291666667*(-fx[x + 3, y + 5] + fx[x + 6, y + 5]) + 0.333333333*(-fx[x + 4, y + 5] + fx[x + 5, y + 5]))*fn[x + 5, y + 5])))
   Eq(fxy[x + 4, y + 4], 0.5*(r1[y] + r1[y + 1]))

A worked example: the elastic stiffness

The constitutive law of linear elasticity, in Voigt notation, reads \(\boldsymbol{\sigma} = \mathbf{C}\,\boldsymbol{\varepsilon}\) with \(\boldsymbol{\sigma}\) and \(\boldsymbol{\varepsilon}\) six-component vectors and \(\mathbf{C}\) a symmetric \(6 \times 6\) stiffness. On the standard staggered grid for elastodynamics:

Voigt index Field Location staggered
1, 2, 3 normal cell centre NODE
4 shear \(yz\) \(yz\)-edge (y, z)
5 shear \(xz\) \(xz\)-edge (x, z)
6 shear \(xy\) \(xy\)-edge (x, y)

and the stiffness coefficients \(C_{ij}\) live at the cell centre.

'symmetric' reproduces the matrix factorisation

\[ \sigma_i \;=\; \mathbf{I}\Bigl(\sum_j C_{ij}\,\mathbf{I}^{\!\top}\!\varepsilon_j\Bigr), \]

with \(\mathbf{I}\) averaging from the cell centre up to \(\sigma_i\)’s location and \(\mathbf{I}^{\!\top}\) averaging \(\varepsilon_j\) back to the centre. Let’s build it in Devito.

from devito import Function, Grid, NODE, Eq, Operator

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

locs = {1: NODE, 2: NODE, 3: NODE,
        4: (y, z), 5: (x, z), 6: (x, y)}


def F(name, stag):
    return Function(name=name, grid=grid, space_order=4, staggered=stag)


sigma = {i: F(f's{i}', locs[i]) for i in range(1, 7)}
eps = {i: F(f'e{i}', locs[i]) for i in range(1, 7)}
C = {(i, j): F(f'C{i}{j}', NODE) for i in range(1, 7) for j in range(1, 7)}

A small helper that compiles a single-equation operator in either mode and returns its symbolic update statement, so we can compare modes side by side:

def show(eq, mode):
    op = Operator(eq, sym_opt={'interp-mode': mode})
    [u] = [n.expr for n in FindNodes(Expression).visit(op)
           if n.expr.lhs.function is eq.lhs.function]
    return u

Normal-normal block \(i, j \in \{1, 2, 3\}\). Everything is at the cell centre, so no interpolation is needed and the two modes agree:

eq = Eq(sigma[1], C[(1, 1)] * eps[1])
print('direct   :', show(eq, 'direct'))
print('symmetric:', show(eq, 'symmetric'))
direct   : Eq(s1[x + 4, y + 4, z + 4], C11[x + 4, y + 4, z + 4]*e1[x + 4, y + 4, z + 4])
symmetric: Eq(s1[x + 4, y + 4, z + 4], C11[x + 4, y + 4, z + 4]*e1[x + 4, y + 4, z + 4])

Normal \(\times\) shear \(i \in \{1, 2, 3\}\), \(j \in \{4, 5, 6\}\). The target \(\sigma_1\) already sits where \(C_{14}\) does (the cell centre), so the symmetric mode falls back to the direct one - only \(\varepsilon_4\) needs to be averaged down:

eq = Eq(sigma[1], C[(1, 4)] * eps[4])
print('direct   :', show(eq, 'direct'))
print('symmetric:', show(eq, 'symmetric'))
direct   : Eq(s1[x + 4, y + 4, z + 4], (0.25*e4[x + 4, y + 3, z + 3] + 0.25*e4[x + 4, y + 4, z + 3] + 0.25*e4[x + 4, y + 3, z + 4] + 0.25*e4[x + 4, y + 4, z + 4])*C14[x + 4, y + 4, z + 4])
symmetric: Eq(s1[x + 4, y + 4, z + 4], (0.25*e4[x + 4, y + 3, z + 3] + 0.25*e4[x + 4, y + 4, z + 3] + 0.25*e4[x + 4, y + 3, z + 4] + 0.25*e4[x + 4, y + 4, z + 4])*C14[x + 4, y + 4, z + 4])

Shear \(\times\) normal \(i \in \{4, 5, 6\}\), \(j \in \{1, 2, 3\}\). Now \(\sigma_4\) sits at the \(yz\)-edge while \(C_{41}\) and \(\varepsilon_1\) are both at the cell centre. The two modes differ:

  • 'direct' averages each factor to the edge independently - giving a product of averages.
  • 'symmetric' forms the product at the cell centre first and averages the result once - giving an average of products.
eq = Eq(sigma[4], C[(4, 1)] * eps[1])
print('direct   :', show(eq, 'direct'))
print()
print('symmetric:', show(eq, 'symmetric'))
direct   : Eq(s4[x + 4, y + 4, z + 4], (0.25*C41[x + 4, y + 4, z + 4] + 0.25*C41[x + 4, y + 5, z + 4] + 0.25*C41[x + 4, y + 4, z + 5] + 0.25*C41[x + 4, y + 5, z + 5])*(0.25*e1[x + 4, y + 4, z + 4] + 0.25*e1[x + 4, y + 5, z + 4] + 0.25*e1[x + 4, y + 4, z + 5] + 0.25*e1[x + 4, y + 5, z + 5]))

symmetric: Eq(s4[x + 4, y + 4, z + 4], 0.25*C41[x + 4, y + 4, z + 4]*e1[x + 4, y + 4, z + 4] + 0.25*C41[x + 4, y + 5, z + 4]*e1[x + 4, y + 5, z + 4] + 0.25*C41[x + 4, y + 4, z + 5]*e1[x + 4, y + 4, z + 5] + 0.25*C41[x + 4, y + 5, z + 5]*e1[x + 4, y + 5, z + 5])

Shear \(\times\) shear off-diagonal \(i \neq j\), both in \(\{4, 5, 6\}\). Every factor sits at a different location. This is where 'symmetric' shines: the discrete operator keeps the \(\mathbf{I}\,(C_{ij}\,\mathbf{I}^{\!\top}\varepsilon_j)\) structure that preserves the self-adjointness of \(\mathbf{C} = \mathbf{C}^{\!\top}\).

eq = Eq(sigma[4], C[(4, 5)] * eps[5])
print('direct   :', show(eq, 'direct'))
print()
print('symmetric:', show(eq, 'symmetric'))
direct   : Eq(s4[x + 4, y + 4, z + 4], (0.25*C45[x + 4, y + 4, z + 4] + 0.25*C45[x + 4, y + 5, z + 4] + 0.25*C45[x + 4, y + 4, z + 5] + 0.25*C45[x + 4, y + 5, z + 5])*(0.25*e5[x + 3, y + 4, z + 4] + 0.25*e5[x + 4, y + 4, z + 4] + 0.25*e5[x + 3, y + 5, z + 4] + 0.25*e5[x + 4, y + 5, z + 4]))

symmetric: Eq(s4[x + 4, y + 4, z + 4], 0.5*(r0[z] + r0[z + 1]))

Shear \(\times\) shear diagonal \(i = j\). The target and \(\varepsilon_i\) already share their staggering; only \(C_{ii}\) needs to be brought to the edge, and that happens implicitly through the standard substitution. The two modes agree:

eq = Eq(sigma[4], C[(4, 4)] * eps[4])
print('direct   :', show(eq, 'direct'))
print('symmetric:', show(eq, 'symmetric'))
direct   : Eq(s4[x + 4, y + 4, z + 4], (0.25*C44[x + 4, y + 4, z + 4] + 0.25*C44[x + 4, y + 5, z + 4] + 0.25*C44[x + 4, y + 4, z + 5] + 0.25*C44[x + 4, y + 5, z + 5])*e4[x + 4, y + 4, z + 4])
symmetric: Eq(s4[x + 4, y + 4, z + 4], (0.25*C44[x + 4, y + 4, z + 4] + 0.25*C44[x + 4, y + 5, z + 4] + 0.25*C44[x + 4, y + 4, z + 5] + 0.25*C44[x + 4, y + 5, z + 5])*e4[x + 4, y + 4, z + 4])

Self-adjointness of the discrete stiffness

The point of 'symmetric' is that the discrete operator inherits the self-adjointness of its continuous counterpart. We can verify this with the standard adjoint dot-product test: if the discrete \(\mathbf{C}\) is its own transpose, then for any two strain-shaped fields \(e_1\) and \(t_2\),

\[ \langle e_1,\,\mathbf{C}\,t_2 \rangle \;=\; \langle \mathbf{C}\,e_1,\,t_2 \rangle. \]

We assemble random symmetric \(C_{ij}\), random \(e_1\) and \(t_2\), compute both sides, and compare:

# NBVAL_IGNORE_OUTPUT
import numpy as np

np.random.seed(1234)

# Random symmetric stiffness
for i in range(1, 7):
    for j in range(i, 7):
        C[(i, j)].data[:] = np.random.rand(*C[(i, j)].shape)
        if j != i:
            C[(j, i)].data[:] = C[(i, j)].data

# Two right-hand-side strain fields
e1 = {i: F(f'e1_{i}', locs[i]) for i in range(1, 7)}
t2 = {i: F(f't2_{i}', locs[i]) for i in range(1, 7)}
# ... and the outputs
t1 = {i: F(f't1_{i}', locs[i]) for i in range(1, 7)}
e2 = {i: F(f'e2_{i}', locs[i]) for i in range(1, 7)}

for i in range(1, 7):
    e1[i].data[:] = 2 * np.random.rand(*e1[i].shape) - 1
    t2[i].data[:] = 2 * np.random.rand(*t2[i].shape) - 1

eqns = []
for i in range(1, 7):
    eqns.append(Eq(t1[i], sum(C[(i, j)] * e1[j] for j in range(1, 7))))
    eqns.append(Eq(e2[i], sum(C[(i, j)] * t2[j] for j in range(1, 7))))


def run(mode):
    Operator(eqns, sym_opt={'interp-mode': mode}).apply()
    lhs = sum(float(np.dot(e1[i].data.flatten(), e2[i].data.flatten()))
              for i in range(1, 7))
    rhs = sum(float(np.dot(t1[i].data.flatten(), t2[i].data.flatten()))
              for i in range(1, 7))
    return lhs, rhs
sym_lhs, sym_rhs = run('symmetric')
sym_rel = abs(sym_lhs - sym_rhs) / max(abs(sym_lhs), abs(sym_rhs))

dir_lhs, dir_rhs = run('direct')
dir_rel = abs(dir_lhs - dir_rhs) / max(abs(dir_lhs), abs(dir_rhs))
NUMA domain count autodetection failed, assuming 1
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
# NBVAL_IGNORE_OUTPUT
print(f"symmetric: <e1, C*t2> = {sym_lhs:.6f}, <C*e1, t2> = {sym_rhs:.6f}")
print(f"           relative difference = {sym_rel:.2e}")
print(f"direct   : <e1, C*t2> = {dir_lhs:.6f}, <C*e1, t2> = {dir_rhs:.6f}")
print(f"           relative difference = {dir_rel:.2e}")
symmetric: <e1, C*t2> = 29.060263, <C*e1, t2> = 29.060250
           relative difference = 4.35e-07
direct   : <e1, C*t2> = 27.701068, <C*e1, t2> = 16.359724
           relative difference = 4.09e-01

The 'symmetric' discretisation passes the dot-product test to float precision; the 'direct' one does not, because each factor was interpolated independently and the resulting discrete operator is no longer the transpose of itself. A quick assertion makes the contrast machine-checkable (and lets CI catch regressions even though the printed floats above can drift in the last digit between platforms):

assert sym_rel < 1e-5, sym_rel    # adjoint identity holds for 'symmetric'
assert dir_rel > 1e-2, dir_rel    # 'direct' is not self-adjoint

When to use which

Situation Mode
Acoustic / scalar-wave equations 'direct'
Elastic stress-strain or any \(\mathbf{I}\,\mathbf{A}\,\mathbf{I}^{\!\top}\) operator 'symmetric'
Adjoint-state inversion needing exact discrete adjoint 'symmetric'
Any equation where one factor already matches the target staggering either

'direct' is the default because it is the cheaper and smaller stencil; pick 'symmetric' deliberately when you need the adjoint structure preserved.

Back to top