Symbolic coefficients on staggered stencils not working correctly

[EdCaunt]

Setting coefficents=symbolic on a function with staggered and not specifying stencil coefficients results in derivatives defaulting to incorrect values. It appears that specifying coefficients does not change them from this incorrect default.

A short example:

from devito import Grid, Function, NODE, Eq, Operator, Coefficient, Substitutions
import numpy as np
shape = (11,)
extent = (10.,)
grid = Grid(shape=shape, extent=extent)
x = grid.dimensions[0]
f = Function(name='f', grid=grid, space_order=2, staggered=NODE)
g = Function(name='g', grid=grid, space_order=2, staggered=NODE,
             coefficients='symbolic')
h = Function(name='h', grid=grid, space_order=2, staggered=NODE,
             coefficients='symbolic')
vf = Function(name='v_f', grid=grid, space_order=2, staggered=x)
vg = Function(name='v_g', grid=grid, space_order=2, staggered=x)
vh = Function(name='v_h', grid=grid, space_order=2, staggered=x)
h_coeffs = Coefficient(1, h, x, np.array([0., -1., 1.]))
coeffs = Substitutions(h_coeffs)
eqf = Eq(vf, f.dx)
eqg = Eq(vg, g.dx)
eqh = Eq(vh, h.dx, coefficients=coeffs)
# Should all result in the same kernel
opf = Operator(eqf)
print(opf.ccode)
# Different stencil to the normal one. Should default to same
opg = Operator(eqg)
print(opg.ccode)
# Different stencil to normal one. Manually set to same
oph = Operator(eqh)
print(oph.ccode)

One would expect all of these operators to have the same ccode (ignoring differences in variable names). However, this is not the case. The latter two do not have the same stencil operation. Eq(vf, f.dxc).evaluate returns Eq(v_f(x + h_x/2), -f(x)/h_x + f(x + h_x)/h_x) whilst Eq(vg, g.dx).evaluate and Eq(vh, h.dx, coefficients=coeffs) both return Eq(v_g(x + h_x/2), g(x + h_x)/(2*h_x))

[FabioLuporini]

some more discussion here: https://devitocodes.slack.com/archives/C3W8GF3H9/p1634720363011200?thread_ts=1634717124.009900&cid=C3W8GF3H9

[EdCaunt]

This was also found to affect substitutions at forward and backward time steps (although this may have been fixed since?). However, a new MFE is:

import numpy as np
from devito import *

grid = Grid(shape=(10, 10))
f = TimeFunction(name='f', grid=grid, space_order=2, coefficients='symbolic')

x, y = grid.dimensions
subs = Substitutions(
    Coefficient(1, f, x, np.array([2, 2, 2]) / x.spacing),
    Coefficient(1, f, y, np.array([2, 2, 2]) / y.spacing)
)

print(grad(f)._evaluate().subs(subs.rules))  
# --> Vector((2/h_x)*f(t, x, y) + (2/h_x)*f(t, x - h_x, y) + (2/h_x)*f(t, x + h_x, y), (2/h_y)*f(t, x, y) + (2/h_y)*f(t, x, y - h_y) + (2/h_y)*f(t, x, y + h_y))

print(grad(f, shift=0.5)._evaluate().subs(subs.rules))
# --> Vector(f(t, x, y)*W(x, 1, f(t, x, y), x + 0.5*h_x) + f(t, x + h_x, y)*W(x + h_x, 1, f(t, x, y), x + 0.5*h_x), f(t, x, y)*W(y, 1, f(t, x, y), y + 0.5*h_y) + f(t, x, y + h_y)*W(y + h_y, 1, f(t, x, y), y + 0.5*h_y))