Understanding devito conventions and defaults

[cdelavegamartin]

I'm just getting started with devito and I have some questions:

Switch of fd_order for first derivative

I would have expected the four following expressions to produce the expression: $\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$

from devito import Grid, Function
grid = Grid(shape=(10, ))
u = Function(name='u', grid=grid, space_order=2)
display(u.dx.evaluate)
display(u.dx(fd_order=2).evaluate)
display(u.dx(fd_order=2, side=centered).evaluate)
display(u.dxc(fd_order=2).evaluate)

$\displaystyle - \frac{u(x)}{h_{x}} + \frac{u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{u(x)}{h_{x}} + \frac{u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$

I understand this is in order to avoid potential instability with the centered scheme but why is not done when higher space_order or fd_order is used? (only happens when not specifying centered, even though the default is centered) eg

u = Function(name='u', grid=grid, space_order=3)
display(u.dx.evaluate)
display(u.dx(fd_order=2).evaluate)
display(u.dx(fd_order=3).evaluate)
display(u.dxc(fd_order=2).evaluate)
display(u.dxc(fd_order=3).evaluate)

$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{u(x)}{h_{x}} + \frac{u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$

right vs left one sided fd convention

For the left side first order accurate first derivative, the stencils are what I would expect

u = Function(name='u', grid=grid, space_order=1)
display(u.dx(side=left,).evaluate)
display(u.dxl.evaluate)
display(u.dxl(fd_order=1).evaluate)

$\displaystyle \frac{u(x)}{h_{x}} - \frac{u(x - h_x)}{h_{x}}$
$\displaystyle \frac{u(x)}{h_{x}} - \frac{u(x - h_x)}{h_{x}}$
$\displaystyle \frac{u(x)}{h_{x}} - \frac{u(x - h_x)}{h_{x}}$

But the right sided ones are displaced one grid point

u = Function(name='u', grid=grid, space_order=1)
display(u.dx(side=right,).evaluate)
display(u.dxr.evaluate)
display(u.dxr(fd_order=1).evaluate)

$\displaystyle - \frac{u(x + h_x)}{h_{x}} + \frac{u(x + 2h_x)}{h_{x}}$
$\displaystyle - \frac{u(x + h_x)}{h_{x}} + \frac{u(x + 2h_x)}{h_{x}}$
$\displaystyle - \frac{u(x + h_x)}{h_{x}} + \frac{u(x + 2h_x)}{h_{x}}$

I would have expected
$\displaystyle - \frac{u(x)}{h_{x}} + \frac{u(x + h_x)}{h_{x}}$

Max differentiation order of u not throwing errors when space_order=1

u = Function(name='u', grid=grid, space_order=1)
display(u.dxc(fd_order=1).evaluate)
display(u.dx(fd_order=2).evaluate)
display(u.dxc(fd_order=2).evaluate)
display(u.dxc(fd_order=3).evaluate)
display(u.dxc(fd_order=4).evaluate)

$\displaystyle - \frac{u(x)}{h_{x}} + \frac{u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{u(x)}{h_{x}} + \frac{u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$
$\displaystyle - \frac{0.5 u(x - h_x)}{h_{x}} + \frac{0.5 u(x + h_x)}{h_{x}}$
$\displaystyle \frac{0.0833333333 u(x - 2h_x)}{h_{x}} - \frac{0.666666667 u(x - h_x)}{h_{x}} + \frac{0.666666667 u(x + h_x)}{h_{x}} - \frac{0.0833333333 u(x + 2h_x)}{h_{x}}$

But

u = Function(name='u', grid=grid, space_order=2)
display(u.dxc(fd_order=3).evaluate)

Produces ValueError: Function does not support x-derivative with `fd_order` 3

I would be very grateful to understand the reason for these conventions/defaults, as I keep questioning if what is happening under the hood for more complex operators is what I'm actually expecting.