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})
    matches = [
        n.expr
        for n in FindNodes(Expression).visit(op)
        if n.expr.lhs.function is eq.lhs.function
    ]
    [u] = matches
    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

Why a gradient needs it

The stiffness above is one \(\mathbf{I}\,\mathbf{A}\,\mathbf{I}^{\!\top}\) operator among others. The case that forces the issue in practice is a gradient: an adjoint-state method accumulates \(\partial F/\partial m\) into a field that lives where the model parameter lives – the node – out of wavefields that live on faces and corners.

That accumulation is a transpose, and it is the transpose of something the forward already does. Where a node-centred parameter multiplies a staggered quantity, the forward reads it through an interpolation,

\[ v_i \mathrel{+}= \Delta t \; (\mathbf{I}\,b)_i \; (\nabla\!\cdot\tau)_i , \]

so the derivative with respect to the nodal \(b\) spreads the staggered residual back over the nodes that fed it, \(\mathbf{I}^{\!\top}\). Forming the product where its operands are and interpolating once is exactly that. Interpolating the operands one at a time instead computes \(\langle a\rangle\langle b\rangle\) where the transpose needs \(\langle ab\rangle\) – second-order accurate, and not a transpose. The symptom is not a wrong-looking gradient; it is a gradient whose Taylor test converges at first order instead of second.

Here is that pattern on its own: a node-centred m read at a staggered location on the way in, and the matching accumulation on the way back.

# NBVAL_IGNORE_OUTPUT
import numpy as np
from devito import Grid, Function, Eq, Operator, NODE

np.random.seed(0)
grid = Grid(shape=(24, 24), extent=(1., 1.))
x, y = grid.dimensions


def adjoint_test(loc, mode):
    m = Function(name='m', grid=grid, space_order=4, staggered=NODE)
    g = Function(name='g', grid=grid, space_order=4, staggered=NODE)
    S = Function(name='S', grid=grid, space_order=4, staggered=loc)
    t = Function(name='t', grid=grid, space_order=4, staggered=loc)
    q = Function(name='q', grid=grid, space_order=4, staggered=loc)

    # support strictly inside, so no edge term enters the identity
    inner = (slice(6, -6), slice(6, -6))
    for f in (m, S, q):
        f.data[:] = 0.
        f.data[inner] = np.random.rand(*f.data[inner].shape) - .5

    # forward: the propagation, always the compact `'direct'` stencil
    Operator([Eq(t, m * S)], sym_opt={'interp-mode': 'direct'}).apply()
    # adjoint: the accumulation back onto the node, in the mode under test
    Operator([Eq(g, q * S)], sym_opt={'interp-mode': mode}).apply()

    lhs = float(np.sum(np.array(t.data) * np.array(q.data)))   # <m*S, q>
    rhs = float(np.sum(np.array(m.data) * np.array(g.data)))   # <m, (q*S)^T>
    return abs(lhs - rhs) / max(abs(lhs), abs(rhs))


for name, loc in (('x-face', (x,)), ('corner', (x, y))):
    for mode in ('direct', 'symmetric'):
        print(
            f"{name:7s} {mode:9s} relative adjoint error "
            f"{adjoint_test(loc, mode):.2e}"
        )
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
x-face  direct    relative adjoint error 4.06e-01
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
x-face  symmetric relative adjoint error 2.07e-07
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
corner  direct    relative adjoint error 8.79e-01
Operator `Kernel` ran in 0.01 s
corner  symmetric relative adjoint error 7.36e-07

Only the accumulation changes mode – the forward keeps the compact 'direct' stencil throughout, as a propagation should. With that, 'symmetric' satisfies \(\langle \mathbf{I}m\,S, q\rangle = \langle m, \mathbf{I}^{\!\top}(qS)\rangle\) to float precision on faces and on corners alike, while 'direct' is out by tens of per cent. That discrepancy is what turns a second-order gradient into a first-order one.

for loc in ((x,), (x, y)):
    assert adjoint_test(loc, 'symmetric') < 1e-5
    assert adjoint_test(loc, 'direct') > 1e-2
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s
Operator `Kernel` ran in 0.01 s

Asking for it on one equation

An Operator usually wants 'direct' throughout: it is the compact stencil, and re-associating a product around a derivative that already sits on the target would replace that derivative with an interpolated one twice as wide. A gradient accumulated inside the time loop of the propagation it is the adjoint of therefore needs the two modes at once – compact stencils, exact transpose on the accumulation. Eq takes the override:

# NBVAL_IGNORE_OUTPUT
u = Function(name='u', grid=grid, space_order=4, staggered=(x, y))
w = Function(name='w', grid=grid, space_order=4, staggered=(x, y))
acc = Function(name='acc', grid=grid, space_order=4, staggered=NODE)

# the Operator is 'direct'; this one equation is not
op = Operator(
    [Eq(acc, acc + u * w, interp_mode='symmetric')],
    sym_opt={'interp-mode': 'direct'},
)
print(
    'symmetric placement on a direct Operator:',
    str(op) == str(
        Operator(
            [Eq(acc, acc + u * w)],
            sym_opt={'interp-mode': 'symmetric'},
        )
    ),
)
symmetric placement on a direct Operator: True

When to use which

Situation Mode
Acoustic / scalar-wave equations 'direct'
Elastic stress-strain or any \(\mathbf{IAI}^{\!\top}\) operator 'symmetric'
Adjoint-state inversion needing exact discrete adjoint 'symmetric'
A gradient accumulated from staggered fields onto a nodal parameter '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 – per Operator with sym_opt, or per equation with Eq(..., interp_mode=...) when only the accumulation needs it.

Back to top