Finite-difference evaluation: from Eq to stencil

When you write an equation like Eq(target, expr) in Devito and pass it to Operator, a fair bit happens before any C is generated. The right-hand side is evaluated at the location of the left-hand side: derivatives get their x0 set to the LHS’s sample point, and any factor that sits at a different staggering is interpolated.

This tutorial unpacks that step. The goal is to make the rules predictable:

We work with Function, Eq, and Operator only - no internal calls. For the deeper 'direct' vs 'symmetric' interpolation modes used on products of staggered fields, see 08_staggered_interp.ipynb.

The lowering pipeline

The compile-time hook is Eq._evaluate. In essence:

def _evaluate(self, **kw):
    lhs = self.lhs._evaluate(**kw)
    rhs = self.rhs._eval_at(self.lhs, **kw)._evaluate(**kw)
    return Eq(lhs, rhs)

Two passes happen on the RHS:

  1. _eval_at(lhs) moves the symbolic expression to the LHS’s sample location. For a Derivative, this sets the x0 keyword. For a Function, this returns either the same function (if it already lives at the right place) or a 0-order FD interpolation stencil.
  2. _evaluate() then expands every derivative into its concrete FD stencil at whatever x0 it now carries.

So the rule of thumb is: when you write Eq(target, expr), every sub-expression of expr is asked “please give me your value at target.indices_ref”. Whether that triggers averaging depends on what the sub-expression is.

from devito import Function, Grid, NODE, Eq, 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)]


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

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))

for f in (fn, fx, fy, fxy):
    print(f.name, '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)

Where a derivative naturally lives

A finite-difference derivative inherits the grid of its operand. Its natural sample point depends on the stencil it carries. With Devito’s default centred stencil, the derivative of a Function f is built around f’s indices: differentiating an f at NODE produces a derivative at NODE; differentiating an f at x + h_x/2 produces a derivative at x + h_x/2. The stencil offsets are what shifts, not the sample point.

Before any equation context, f.dx carries no x0 and Devito will use its natural location:

# NBVAL_IGNORE_OUTPUT
print('fn.dx symbolic :', fn.dx)
print()
print('fn.dx.evaluate :', fn.dx.evaluate)
fn.dx symbolic : Derivative(fn(x, y), x)

fn.dx.evaluate : 0.0833333333*fn(x - 2*h_x, y)/h_x - 0.666666667*fn(x - h_x, y)/h_x + 0.666666667*fn(x + h_x, y)/h_x - 0.0833333333*fn(x + 2*h_x, y)/h_x

_eval_at on a Derivative: shifting x0

When Devito sees Eq(target, f.dx), it calls f.dx._eval_at(target). This sets the derivative’s x0 to target.indices_ref so that the FD stencil is built around the target’s sample point.

Crucially, the operand isn’t interpolated - only the stencil offsets shift. That keeps the order of accuracy and the stencil shape clean. So fn.dx at x + h_x/2 is not an average of fn.dx values; it is a stencil over fn samples reshuffled to land on x + h_x/2.

# NBVAL_IGNORE_OUTPUT
# Same operand, two different targets, two different x0s
print('--- fn.dx evaluated naturally (no Eq context)')
print(fn.dx.evaluate)
print()
print('--- fn.dx as it lowers inside Eq(fx, fn.dx) (target at x + h_x/2)')
print(fn.dx._eval_at(fx).evaluate)
--- fn.dx evaluated naturally (no Eq context)
0.0833333333*fn(x - 2*h_x, y)/h_x - 0.666666667*fn(x - h_x, y)/h_x + 0.666666667*fn(x + h_x, y)/h_x - 0.0833333333*fn(x + 2*h_x, y)/h_x

--- fn.dx as it lowers inside Eq(fx, fn.dx) (target at x + h_x/2)
-1.125*fn(x, y)/h_x + 0.0416666667*fn(x - h_x, y)/h_x + 1.125*fn(x + h_x, y)/h_x - 0.0416666667*fn(x + 2*h_x, y)/h_x

The same shift happens implicitly inside an Operator. Comparing two equations whose only difference is the LHS:

# NBVAL_IGNORE_OUTPUT
op1 = Operator(Eq(fn, fn.dx))  # target at NODE
op2 = Operator(Eq(fx, fn.dx))  # target at x + h_x/2

print('--- Eq(fn, fn.dx)')
for u in updates(op1):
    print('  ', u)
print()
print('--- Eq(fx, fn.dx) (note the different fn-indices)')
for u in updates(op2):
    print('  ', u)
--- Eq(fn, fn.dx)
   Eq(r0, 1/h_x)
   Eq(fn[x + 4, y + 4], r0*(0.0833333333*(fn[x + 2, y + 4] - fn[x + 6, y + 4]) + 0.666666667*(-fn[x + 3, y + 4] + fn[x + 5, y + 4])))

--- Eq(fx, fn.dx) (note the different fn-indices)
   Eq(r0, 1/h_x)
   Eq(fx[x + 4, y + 4], r0*(0.0416666667*(fn[x + 3, y + 4] - fn[x + 6, y + 4]) + 1.125*(-fn[x + 4, y + 4] + fn[x + 5, y + 4])))

_eval_at on a Function: interpolation

A bare Function on the RHS is a different story. If its staggering does not match the target’s, Devito emits a 0-order FD interpolation operator that averages it onto the target’s sample points. In 1D this is the two-point average; in higher dimensions it tensors out.

# fx lives at x + h_x/2, fn at NODE - the assignment averages
op = Operator(Eq(fn, fx))
for u in updates(op):
    print('  ', u)
   Eq(fn[x + 4, y + 4], 0.5*(fx[x + 3, y + 4] + fx[x + 4, y + 4]))

When auto-interpolation doesn’t happen

Three things short-circuit the automatic interpolation. Knowing them is the difference between “magical” stencils and predictable ones.

1. The two sides already share a staggering

If lhs.indices_ref == rhs.indices_ref there is nothing to do, so the RHS comes through unchanged. This includes the common staggered=None case: passing None is not a “don’t interpolate” flag - it simply means “no staggering specified”, which Devito treats as NODE. So staggered=None and staggered=NODE produce identical lowerings, and both still trigger interpolation when the RHS lives elsewhere.

# Same staggering on both sides - no interpolation
fx2 = Function(name='fx2', grid=grid, space_order=4, staggered=x)
op = Operator(Eq(fx, fx2))
print('Eq(fx, fx2)   :', updates(op)[0])

# staggered=None really is NODE - interpolation still kicks in
fnone = Function(name='fnone', grid=grid, space_order=4, staggered=None)
op = Operator(Eq(fnone, fx))
print('Eq(fnone, fx) :', updates(op)[0])
print('-> identical to Eq(fn, fx) earlier')
Eq(fx, fx2)   : Eq(fx[x + 4, y + 4], fx2[x + 4, y + 4])
Eq(fnone, fx) : Eq(fnone[x + 4, y + 4], 0.5*(fx[x + 3, y + 4] + fx[x + 4, y + 4]))
-> identical to Eq(fn, fx) earlier

2. interp_order=0 on the Function

Passing interp_order=0 when you build a Function opts it out of averaging. The RHS is then sampled at the LHS’s grid points without any average. This is useful for piecewise-constant material parameters (a velocity map, a mask) that should not be smoothed.

fx0 = Function(name='fx0', grid=grid, space_order=4,
               staggered=x, interp_order=0)
op = Operator(Eq(fn, fx0))
print('Eq(fn, fx0) :', updates(op)[0])
print('-> fx0 sampled directly, no average')
Eq(fn, fx0) : Eq(fn[x + 4, y + 4], fx0[x + 4, y + 4])
-> fx0 sampled directly, no average

3. An explicit x0 on a Derivative

If you pass x0 to f.dx(x0=...) yourself, Devito will not overwrite it. Use this when you want to control where the derivative is sampled independently of the LHS - for example to keep a centred stencil even when the target is staggered.

# NBVAL_IGNORE_OUTPUT
# Default vs. explicit x0
default = Operator(Eq(fx, fn.dx))               # x0 = x + h_x/2
explicit = Operator(Eq(fx, fn.dx(x0={x: x})))   # x0 = x  (centred)
print('default  :', updates(default)[-1])
print('explicit :', updates(explicit)[-1])
default  : Eq(fx[x + 4, y + 4], r0*(0.0416666667*(fn[x + 3, y + 4] - fn[x + 6, y + 4]) + 1.125*(-fn[x + 4, y + 4] + fn[x + 5, y + 4])))
explicit : Eq(fx[x + 4, y + 4], r0*(0.0833333333*(fn[x + 2, y + 4] - fn[x + 6, y + 4]) + 0.666666667*(-fn[x + 3, y + 4] + fn[x + 5, y + 4])))

Mixed equations

Real equations combine all of the above. The lowering walks the RHS recursively: each factor in a product is evaluated at the LHS’s location; each summand in an Add is evaluated independently. Derivatives keep their x0-shift behaviour, Functions their averaging behaviour.

# NBVAL_IGNORE_OUTPUT
# A product of a Function and a derivative, target at the corner
eq = Eq(fxy, fn * fx.dx)
op = Operator(eq)
for u in updates(op):
    print('  ', u)
   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 is averaged to the corner (four-point average); fx.dx has its x0 set to the corner so the FD stencil is already there. The two interpolations happen independently - that is the 'direct' mode discussed in 08_staggered_interp.ipynb. The 'symmetric' mode rewires this product into an I (a I^T b) form when self-adjointness matters.

A sum behaves the same way, summand by summand:

# NBVAL_IGNORE_OUTPUT
# fn lives at NODE; fy.dy is naturally at NODE too (a y-staggered
# field's y-derivative lands at NODE). Target fx is at x + h_x/2.
eq = Eq(fx, fn + fy.dy)
op = Operator(eq)
for u in updates(op):
    print('  ', u)
   Eq(r0, 1/h_y)
   Eq(fx[x + 4, y + 4], 0.0208333333*(r0*(fy[x + 4, y + 2] + fy[x + 5, y + 2]) - r0*(fy[x + 4, y + 5] + fy[x + 5, y + 5])) + 0.5625*(-r0*(fy[x + 4, y + 3] + fy[x + 5, y + 3]) + r0*(fy[x + 4, y + 4] + fy[x + 5, y + 4])) + 0.5*(fn[x + 4, y + 4] + fn[x + 5, y + 4]))

Both summands had to be brought to x + h_x/2: fn by averaging in x, and the fy.dy stencil by an x0 shift on its derivative dim plus an interpolation in y (since fy lives at y + h_y/2).

Cheat sheet

Step What Devito does
Eq._evaluate Calls rhs._eval_at(lhs) then ._evaluate().
Derivative._eval_at(target) Sets x0 = target.indices_ref on the derivative.
Function._eval_at(target) Emits a 0-order FD average to target.indices_ref (when staggerings differ).
LHS and RHS share staggering No-op.
staggered=None Treated as NODE. Does not opt out of interpolation.
interp_order=0 Function opts out of averaging - sampled directly at the target.
Explicit x0 on a derivative Devito leaves it alone.

For the sym_opt={'interp-mode': ...} choice that governs how products of staggered fields are projected onto a target, see 08_staggered_interp.ipynb.

CI regression guards

The cells above print symbolic stencils for human reading; their textual form can drift between SymPy versions and is not load-bearing. The following cell pins the behavioural invariants of this tutorial so that a regression in the lowering will fail CI even when output comparison is disabled.

# 1. `staggered=None` lowers identically to `staggered=NODE`.
rhs_none = str(updates(Operator(Eq(fnone, fx)))[0].rhs)
rhs_node = str(updates(Operator(Eq(fn, fx)))[0].rhs)
assert rhs_none == rhs_node, (rhs_none, rhs_node)

# 2. `interp_order=0` opts out of averaging - the RHS is sampled directly.
rhs_io0 = str(updates(Operator(Eq(fn, fx0)))[0].rhs)
assert rhs_io0 == 'fx0[x + 4, y + 4]', rhs_io0

# 3. Same staggering on both sides is a pure pass-through.
rhs_same = str(updates(Operator(Eq(fx, fx2)))[0].rhs)
assert rhs_same == 'fx2[x + 4, y + 4]', rhs_same

# 4. An explicit `x0` on a derivative is not overwritten by `_eval_at` -
#    the two stencils below must differ.
default_rhs = str(updates(Operator(Eq(fx, fn.dx)))[-1].rhs)
explicit_rhs = str(updates(Operator(Eq(fx, fn.dx(x0={x: x}))))[-1].rhs)
assert default_rhs != explicit_rhs

print('all invariants hold')
all invariants hold
Back to top