api: support RSFD

devito/finite_differences/finite_difference.py

@@ -87,6 +87,10 @@ def first_derivative(expr, dim, fd_order=None, side=centered, matvec=direct, x0=
     fd_order = fd_order or expr.space_order
     deriv_order = 1
 
+    # Enforce sable time coefficients
+    if dim.is_Time and coefficients != 'symbolic':
+        coefficients = 'taylor'
+
     return make_derivative(expr, dim, fd_order, deriv_order, side,
                            matvec, x0, coefficients, expand)
 
@@ -202,6 +206,10 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
     if deriv_order == 1 and fd_order == 2 and coefficients != 'symbolic':
         fd_order = 1
 
+    # Enforce sable time coefficients
+    if dim.is_Time and coefficients != 'symbolic':
+        coefficients = 'taylor'
+
     return make_derivative(expr, dim, fd_order, deriv_order, side,
                            matvec, x0, coefficients, expand)
 

devito/finite_differences/rsfd.py

@@ -1,8 +1,9 @@
-from sympy import finite_diff_weights
+from functools import wraps
 
+from devito.types import NODE
 from devito.types.dimension import StencilDimension
 from .differentiable import Weights, DiffDerivative
-from .tools import generate_indices_staggered
+from .tools import generate_indices_staggered, fd_weights_registry
 from .elementary import sqrt
 
 __all__ = ['drot', 'dxrot', 'dyrot', 'dzrot']
@@ -41,22 +42,28 @@ def drot(expr, dim, dir=1, x0=None):
     # Spacing along diagonal
     r = sqrt(sum(d.spacing**2 for d in expr.grid.dimensions))
 
-    # Center point
+    # Center point and indices
     start, indices = generate_indices_staggered(expr, dim, expr.space_order, x0=x0)
-    adim_start = start.subs({dim: 0, dim.spacing: 1})
-    adim_indices = [i.subs({dim: 0, dim.spacing: 1}) for i in indices]
-    coeffs = finite_diff_weights(1, adim_indices, adim_start)[-1][-1]
 
-    i = StencilDimension('i', adim_indices[0]-adim_start, adim_indices[-1]-adim_start)
+    # a-dimensional center for FD coefficients.
+    adim_start = x0.get(dim, expr.indices_ref[dim]).subs({dim: 0, dim.spacing: 1})
+    mid = expr.indices_ref[dim].subs({dim: 0, dim.spacing: 1})
+
+    # a-dimensional indices
+    adim_indices = [i.subs({dim: 0, dim.spacing: 1}) for i in indices]
 
-    # FD weights
+    # FD coeffs
+    # Dispersion reduction weights currently not working as the lsqr
+    # system needs to be setup for the whole stencil
+    coeffs = fd_weights_registry['taylor'](expr, 1, adim_indices, adim_start)
+    i = StencilDimension('i', adim_indices[0]-mid, adim_indices[-1]-mid)
     w0 = Weights(name='w0', dimensions=i, initvalue=coeffs)
 
     # Skip y if 2D
     signs = smapper[dir][::(1 if ndim == 3 else 2)]
 
     # Direction substitutions
-    dim_mapper = {d: d + signs[di]*i*d.spacing - shift(signs[di], adim_start)*d.spacing
+    dim_mapper = {d: d + signs[di]*i*d.spacing - shift(signs[di], mid)*d.spacing
                   for (di, d) in enumerate(expr.grid.dimensions)}
 
     # Create IndexDerivative
@@ -67,6 +74,46 @@ def drot(expr, dim, dir=1, x0=None):
     return deriv/r
 
 
+grid_node = lambda grid: {d: d for d in grid.dimensions}
+all_staggered = lambda grid: {d: d + d.spacing/2 for d in grid.dimensions}
+
+
+def check_staggering(func):
+    """
+    Because of the very specific structure of the 45 degree stencil
+    only two cases can happen:
+
+    1. No staggering. In this case the stencil is center on the node where
+       the Function/Expr is defined and the diagonal is well defined.
+    2. Full staggering. The center node must be either NODE or grid.dimension so that
+       the diagonal align with the staggering of the expression. I.e a NODE center point
+       for a `grid.dimension` staggered expression
+       or a `grid.dimension` center point for a NODE staggered expression.
+
+    Therefore acceptable combinations are:
+    - NODE center point and NODE staggering
+    - grid.dimension center point and grid.dimension staggering
+    - NODE center point and grid.dimension staggering
+    - grid.dimension center point and NODE staggering
+    """
+    @wraps(func)
+    def wrapper(expr, x0=None, expand=True):
+        grid = expr.grid
+        x0 = {k: v for k, v in x0.items() if k.is_Space}
+        if expr.staggered is NODE or expr.staggered is None:
+            cond = x0 == {} or x0 == all_staggered(grid) or x0 == grid_node(grid)
+        elif expr.staggered == grid.dimensions:
+            cond = x0 == {} or x0 == all_staggered(grid) or x0 == grid_node(grid)
+        else:
+            cond = False
+        if cond:
+            return func(expr, x0=x0, expand=expand)
+        else:
+            raise ValueError('Invalid staggering or x0 for rotated finite differences')
+    return wrapper
+
+
+@check_staggering
 def dxrot(expr, x0=None, expand=True):
     x = expr.grid.dimensions[0]
     r = sqrt(sum(d.spacing**2 for d in expr.grid.dimensions))
@@ -80,6 +127,7 @@ def dxrot(expr, x0=None, expand=True):
         return dx45
 
 
+@check_staggering
 def dyrot(expr, x0=None, expand=True):
     y = expr.grid.dimensions[1]
     r = sqrt(sum(d.spacing**2 for d in expr.grid.dimensions))
@@ -93,6 +141,7 @@ def dyrot(expr, x0=None, expand=True):
         return dy45
 
 
+@check_staggering
 def dzrot(expr, x0=None, expand=True):
     z = expr.grid.dimensions[-1]
     r = sqrt(sum(d.spacing**2 for d in expr.grid.dimensions))

tests/test_derivatives.py

@@ -261,6 +261,56 @@ def test_fd_space(self, derivative, space_order):
                     Dpolyvalues[space_border:-space_border])
         assert np.isclose(np.mean(error), 0., atol=1e-3)
 
+    @pytest.mark.parametrize('staggered', [(True, True), (False, False),
+                                           (True, False), (False, True)])
+    @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
+    def test_fd_space_45(self, staggered, space_order):
+        """
+        Rotated derivatives require at least 2D to get access to diagonal points.
+        We create a simple 1D gradient over a 2D grid to check against a polynomial
+        """
+        # dummy axis dimension
+        nx = 100
+        ny = nx
+        xx = np.linspace(-1, 1, nx)
+        dx = xx[1] - xx[0]
+        if staggered[0] and not staggered[1]:
+            xx_s = xx + dx/2
+        elif not staggered[0] and staggered[1]:
+            xx_s = xx - dx/2
+        else:
+            xx_s = xx
+        # Symbolic data
+        grid = Grid(shape=(nx, ny), dtype=np.float32)
+        x = grid.dimensions[0]
+        u = Function(name="u", grid=grid, space_order=space_order,
+                     staggered=None if staggered[0] else grid.dimensions)
+        du = Function(name="du", grid=grid, space_order=space_order,
+                      staggered=None if staggered[1] else grid.dimensions)
+        # Define polynomial with exact fd
+        coeffs = np.ones((space_order,), dtype=np.float32)
+        polynome = sum([coeffs[i]*x**i for i in range(0, space_order)])
+        polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32)
+        # Fill original data with the polynomial values
+        for i in range(ny):
+            u.data[:, i] = polyvalues
+        # True derivative of the polynome
+        Dpolynome = diff(polynome)
+        Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx_s], np.float32)
+        # FD derivative, symbolic
+        u_deriv = getattr(u, 'dx45')
+        # Compute numerical FD
+        stencil = Eq(du, u_deriv)
+        op = Operator(stencil, subs={d.spacing: dx for d in grid.dimensions})
+        op.apply()
+
+        # Check exactness of the numerical derivative except inside space_brd
+        space_border = space_order
+        error = abs(du.data[space_border:-space_border, ny//2] -
+                    Dpolyvalues[space_border:-space_border])
+
+        assert np.isclose(np.mean(error), 0., atol=1e-3)
+
     @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
     @pytest.mark.parametrize('stagger', [centered, left, right])
     # Only test x and t as y and z are the same as x