Custom finite difference weights in DevitoPRO

Custom finite difference weights can be used in Devito to create more efficient operators. This is particularly useful for applications that require high accuracy. Devito provides a simple mathematical interface to define custom weights and apply them to a given derivative. However, in some cases, the user may prefer to define a global finite difference stencil to be used for all derivatives of a given derivative order and finite difference order. We provide an API in DevitoPRO to register your own global finite difference weights.

from devito import Grid, Function
from devitopro import *

custom_w = [-2, -1, 3.14, 1, 2]
grid = Grid((11, 11))

u = TimeFunction(name='u', grid=grid, space_order=4, time_order=1)

Let’s look at the default behavior for reference

u.dx

\(\displaystyle \frac{\partial}{\partial x} u(t, x, y)\)

u.dx.evaluate

\(\displaystyle \frac{0.0833333333 u(t, x - 2*h_x, y)}{h_{x}} - \frac{0.666666667 u(t, x - h_x, y)}{h_{x}} + \frac{0.666666667 u(t, x + h_x, y)}{h_{x}} - \frac{0.0833333333 u(t, x + 2*h_x, y)}{h_{x}}\)

Now we can use the standard devito API for a custom finite difference stencil

u.dx(w=custom_w).evaluate

\(\displaystyle \frac{3.14 u(t, x, y)}{h_{x}} - \frac{2.0 u(t, x - 2*h_x, y)}{h_{x}} - \frac{u(t, x - h_x, y)}{h_{x}} + \frac{u(t, x + h_x, y)}{h_{x}} + \frac{2.0 u(t, x + 2*h_x, y)}{h_{x}}\)

Now let’s instead define a global rule for the coefficents to be used by every function computing the same derivative. The function register_weights allows to register a global finite difference stencil for a given derivative order and finite difference order. The weights can be defined as a list of coefficients, and the function will automatically normalize them by the grid spacing to the power of the derivative order if needed (normalized=False).

#NBVAL_IGNORE_OUTPUT

print(register_weights.__doc__)

Register custom weights for a given derivative order and fd order.

Parameters
----------
deriv_order : int
    The order of the derivative.
fd_order : int
    The order of the finite difference.
weights : Function or Symbol
    The weights to register.
normalized : bool, optional
    whether the weights are normalized or not. If False, the weights
    are divided by spacing**order of the dimension at evaluation.

deriv_order = 1
fd_order = 4
register_weights(deriv_order, fd_order, custom_w, normalized=False)

Now let’s see if u properly picks up the right weights.

u.dx.evaluate

\(\displaystyle \frac{3.14 u(t, x, y)}{h_{x}} - \frac{2.0 u(t, x - 2*h_x, y)}{h_{x}} - \frac{u(t, x - h_x, y)}{h_{x}} + \frac{u(t, x + h_x, y)}{h_{x}} + \frac{2.0 u(t, x + 2*h_x, y)}{h_{x}}\)

We can now define any function and it will use these weights for a first order derivative of any dimension for a 4th order stencil

from IPython.display import display

for i in range(3):
    f = Function(name=f'f{i}', grid=grid, space_order=4, time_order=1)
    for dim in grid.dimensions:
        display(f.dx.evaluate)

\(\displaystyle \frac{3.14 f0(x, y)}{h_{x}} - \frac{2.0 f0(x - 2*h_x, y)}{h_{x}} - \frac{f0(x - h_x, y)}{h_{x}} + \frac{f0(x + h_x, y)}{h_{x}} + \frac{2.0 f0(x + 2*h_x, y)}{h_{x}}\)

\(\displaystyle \frac{3.14 f0(x, y)}{h_{x}} - \frac{2.0 f0(x - 2*h_x, y)}{h_{x}} - \frac{f0(x - h_x, y)}{h_{x}} + \frac{f0(x + h_x, y)}{h_{x}} + \frac{2.0 f0(x + 2*h_x, y)}{h_{x}}\)

\(\displaystyle \frac{3.14 f1(x, y)}{h_{x}} - \frac{2.0 f1(x - 2*h_x, y)}{h_{x}} - \frac{f1(x - h_x, y)}{h_{x}} + \frac{f1(x + h_x, y)}{h_{x}} + \frac{2.0 f1(x + 2*h_x, y)}{h_{x}}\)

\(\displaystyle \frac{3.14 f1(x, y)}{h_{x}} - \frac{2.0 f1(x - 2*h_x, y)}{h_{x}} - \frac{f1(x - h_x, y)}{h_{x}} + \frac{f1(x + h_x, y)}{h_{x}} + \frac{2.0 f1(x + 2*h_x, y)}{h_{x}}\)

\(\displaystyle \frac{3.14 f2(x, y)}{h_{x}} - \frac{2.0 f2(x - 2*h_x, y)}{h_{x}} - \frac{f2(x - h_x, y)}{h_{x}} + \frac{f2(x + h_x, y)}{h_{x}} + \frac{2.0 f2(x + 2*h_x, y)}{h_{x}}\)

\(\displaystyle \frac{3.14 f2(x, y)}{h_{x}} - \frac{2.0 f2(x - 2*h_x, y)}{h_{x}} - \frac{f2(x - h_x, y)}{h_{x}} + \frac{f2(x + h_x, y)}{h_{x}} + \frac{2.0 f2(x + 2*h_x, y)}{h_{x}}\)

You can also clear the registered weights by calling clear_registered_weights. You can either call it with the derivative order and finite difference order as arguments to remove a single entry or without arguments to clear all registered weights.

# Add a second rule
deriv_order = 2
fd_order = 4
register_weights(deriv_order, fd_order, custom_w, normalized=False)
u.dx2.evaluate

\(\displaystyle \frac{3.14 u(t, x, y)}{h_{x}^{2}} - \frac{2.0 u(t, x - 2*h_x, y)}{h_{x}^{2}} - \frac{u(t, x - h_x, y)}{h_{x}^{2}} + \frac{u(t, x + h_x, y)}{h_{x}^{2}} + \frac{2.0 u(t, x + 2*h_x, y)}{h_{x}^{2}}\)

# Clear the second order weights
clear_registered_weights(deriv_order=2, fd_order=4)
u.dx2.evaluate

\(\displaystyle - \frac{2.5 u(t, x, y)}{h_{x}^{2}} - \frac{0.0833333333 u(t, x - 2*h_x, y)}{h_{x}^{2}} + \frac{1.33333333 u(t, x - h_x, y)}{h_{x}^{2}} + \frac{1.33333333 u(t, x + h_x, y)}{h_{x}^{2}} - \frac{0.0833333333 u(t, x + 2*h_x, y)}{h_{x}^{2}}\)

# The first order weights are still registered
u.dx.evaluate

\(\displaystyle \frac{3.14 u(t, x, y)}{h_{x}} - \frac{2.0 u(t, x - 2*h_x, y)}{h_{x}} - \frac{u(t, x - h_x, y)}{h_{x}} + \frac{u(t, x + h_x, y)}{h_{x}} + \frac{2.0 u(t, x + 2*h_x, y)}{h_{x}}\)

# Clear all registered weights
clear_registered_weights()
# Now it should be back to standard taylor weights
u.dx.evaluate

\(\displaystyle \frac{0.0833333333 u(t, x - 2*h_x, y)}{h_{x}} - \frac{0.666666667 u(t, x - h_x, y)}{h_{x}} + \frac{0.666666667 u(t, x + h_x, y)}{h_{x}} - \frac{0.0833333333 u(t, x + 2*h_x, y)}{h_{x}}\)

Back to top