api: fix indices inconsistency in custom FD

devito/finite_differences/coefficients.py

@@ -2,7 +2,7 @@
 from cached_property import cached_property
 
 from devito.finite_differences import Weights, generate_indices
-from devito.finite_differences.tools import numeric_weights, symbolic_weights
+from devito.finite_differences.tools import numeric_weights
 from devito.tools import filter_ordered, as_tuple
 
 __all__ = ['Coefficient', 'Substitutions', 'default_rules']
@@ -131,7 +131,7 @@ class Substitutions(object):
 
     Now define some partial d/dx FD coefficients of the Function u:
 
-    >>> u_x_coeffs = Coefficient(1, u, x, np.array([-0.6, 0.1, 0.6]))
+    >>> u_x_coeffs = Coefficient(2, u, x, np.array([-0.6, 0.1, 0.6]))
 
     And now create our Substitutions object to pass to equation:
 
@@ -141,7 +141,7 @@ class Substitutions(object):
     Now create a Devito equation and pass to it 'subs'
 
     >>> from devito import Eq
-    >>> eq = Eq(u.dt+u.dx, coefficients=subs)
+    >>> eq = Eq(u.dt+u.dx2, coefficients=subs)
 
     When evaluated, the derivatives will use the custom coefficients. We can
     check that by
@@ -231,7 +231,7 @@ def default_rules(obj, functions):
 
     from devito.symbolics.search import retrieve_dimensions
 
-    def generate_subs(deriv_order, function, index):
+    def generate_subs(deriv_order, function, index, indices):
         dim = retrieve_dimensions(index)[0]
 
         if dim.is_Time:
@@ -244,27 +244,15 @@ def generate_subs(deriv_order, function, index):
 
         subs = {}
 
-        mapper = {dim: index}
-        # Get full range of indices and weights
-        indices, x0 = generate_indices(function, dim,
-                                       fd_order, side=None, x0=mapper)
-        sweights = symbolic_weights(function, deriv_order, indices, x0)
-
         # Actual FD used indices and weights
         if deriv_order == 1 and fd_order == 2:
             fd_order = 1
 
-        indices, x0 = generate_indices(function, dim,
-                                       fd_order, side=None, x0=mapper)
-
-        coeffs = numeric_weights(function, deriv_order, indices, x0)
+        coeffs = numeric_weights(function, deriv_order, indices, index)
 
         for (c, i) in zip(coeffs, indices):
             subs.update({function._coeff_symbol(i, deriv_order, function, index): c})
 
-        # Set all unused weights to zero
-        subs.update({w: 0 for w in sweights if w not in subs})
-
         return subs
 
     # Determine which 'rules' are missing
@@ -274,7 +262,11 @@ def generate_subs(deriv_order, function, index):
         for w in i.weights:
             terms.update(w.find(sym))
 
-    args_present = filter_ordered(term.args[1:] for term in terms)
+    args_present = {}
+    for term in terms:
+        key = term.args[1:]
+        idx = term.args[0]
+        args_present.setdefault(key, []).append(idx)
 
     subs = obj.substitutions
     if subs:
@@ -288,7 +280,7 @@ def generate_subs(deriv_order, function, index):
     args_provided = list(set(args_provided))
 
     rules = {}
-    not_provided = []
+    not_provided = {}
     for i0 in args_present:
         if any(i0 == i1 for i1 in args_provided):
             # Perfect match, as expected by the legacy custom coeffs API
@@ -313,10 +305,10 @@ def generate_subs(deriv_order, function, index):
         if mapper:
             rules.update(mapper)
         else:
-            not_provided.append(i0)
+            not_provided.update({i0: args_present[i0]})
 
-    for i in not_provided:
-        rules = {**rules, **generate_subs(*i)}
+    for i, indices in not_provided.items():
+        rules = {**rules, **generate_subs(*i, indices)}
 
     return rules
 

devito/finite_differences/finite_difference.py

@@ -213,7 +213,7 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
     side = None
     # First order derivative with 2nd order FD is strongly discouraged so taking
     # first order fd that is a lot better
-    if deriv_order == 1 and fd_order == 2 and coefficients != 'symbolic':
+    if deriv_order == 1 and fd_order == 2:
         fd_order = 1
 
     # Zeroth order derivative is just the expression itself if not shifted

tests/test_staggered_utils.py

@@ -110,7 +110,6 @@ def test_stagg_fd_composite(expr, expected):
     y0 = y + y.spacing/2  # noqa
     a = Function(name="a", grid=grid, staggered=NODE)  # noqa
     b = Function(name="b", grid=grid, staggered=x)  # noqa
-
     assert eval(expr) == eval(expected)
 
 

tests/test_symbolic_coefficients.py

@@ -16,7 +16,7 @@ class TestSC(object):
     Class for testing symbolic coefficients functionality
     """
 
-    @pytest.mark.parametrize('order', [1, 2, 6])
+    @pytest.mark.parametrize('order', [2, 6])
     @pytest.mark.parametrize('stagger', [True, False])
     def test_default_rules(self, order, stagger):
         """
@@ -40,26 +40,27 @@ def test_default_rules(self, order, stagger):
                eq1.evaluate.evalf(_PRECISION).__repr__())
 
     @pytest.mark.parametrize('expr, sorder, dorder, dim, weights, expected', [
-        ('u.dx', 2, 1, 0, (-0.6, 0.1, 0.6),
-         '0.1*u(x, y) - 0.6*u(x - h_x, y) + 0.6*u(x + h_x, y)'),
+        ('u.dx2', 2, 2, 0, (-0.6, 0.1, 0.6),
+         '0.1*u - 0.6*u._subs(x, x - h_x) + 0.6*u._subs(x, x + h_x)'),
         ('u.dy2', 4, 2, 1, (0.121, -0.223, 1.648, -2.904, 0),
-         '1.648*u(x, y) + 0.121*u(x, y - 2*h_y) - 0.223*u(x, y - h_y) \
-- 2.904*u(x, y + h_y)')])
+         '1.648*u + 0.121*u._subs(y, y - 2*h_y) - 0.223*u._subs(y, y - h_y) \
+- 2.904*u._subs(y, y + h_y)')])
     def test_coefficients(self, expr, sorder, dorder, dim, weights, expected):
         """Test that custom coefficients return the expected result"""
         grid = Grid(shape=(10, 10))
         u = Function(name='u', grid=grid, space_order=sorder, coefficients='symbolic')
-        x = grid.dimensions
+        x, y = grid.dimensions
+        h_x, h_y = grid.spacing_symbols  # noqa
 
         order = dorder
-        dim = x[dim]
+        dim = grid.dimensions[dim]
         weights = np.array(weights)
 
         coeffs = Coefficient(order, u, dim, weights)
 
         eq = Eq(eval(expr), coefficients=Substitutions(coeffs))
         assert isinstance(eq.lhs, Differentiable)
-        assert expected == str(eq.evaluate.lhs)
+        assert sp.simplify(eval(expected) - eq.evaluate.lhs) == 0
 
     def test_function_coefficients(self):
         """Test that custom function coefficients return the expected result"""
@@ -108,20 +109,21 @@ def test_coefficients_w_xreplace(self):
         grid = Grid(shape=(4, 4))
         u = Function(name='u', grid=grid, space_order=2, coefficients='symbolic')
         x = grid.dimensions[0]
+        h_x = x.spacing
 
-        dorder = 1
+        dorder = 2
         weights = np.array([-0.6, 0.1, 0.6])
 
         coeffs = Coefficient(dorder, u, x, weights)
 
         c = sp.Symbol('c')
 
-        eq = Eq(u.dx+c, coefficients=Substitutions(coeffs))
+        eq = Eq(u.dx2+c, coefficients=Substitutions(coeffs))
         eq = eq.xreplace({c: 2})
 
-        expected = '0.1*u(x, y) - 0.6*u(x - h_x, y) + 0.6*u(x + h_x, y) + 2'
+        expected = 0.1*u - 0.6*u._subs(x, x - h_x) + 0.6*u.subs(x, x + h_x) + 2
 
-        assert expected == str(eq.evaluate.lhs)
+        assert sp.simplify(expected - eq.evaluate.lhs) == 0
 
     @pytest.mark.parametrize('order', [1, 2, 6, 8])
     @pytest.mark.parametrize('extent', [1., 10., 100.])
@@ -439,10 +441,10 @@ def test_nested_subs(self):
         coeffs0 = np.array([100, 100, 100])
         coeffs1 = np.array([200, 200, 200])
 
-        subs = Substitutions(Coefficient(1, p, x, coeffs0),
-                             Coefficient(1, p, y, coeffs1))
+        subs = Substitutions(Coefficient(2, p, x, coeffs0),
+                             Coefficient(2, p, y, coeffs1))
 
-        eq = Eq(p.forward, p.dx.dy, coefficients=subs)
+        eq = Eq(p.forward, p.dx2.dy2, coefficients=subs)
 
         mul = lambda e: sp.Mul(e, 200, evaluate=False)
         term0 = mul(p*100 +
@@ -456,7 +458,7 @@ def test_nested_subs(self):
                     p.subs({x: x+hx, y: y+hy})*100)
 
         # `str` simply because some objects are of type EvalDerivative
-        assert str(eq.evaluate.rhs) == str(term0 + term1 + term2)
+        assert sp.simplify(eq.evaluate.rhs - term0 - term1 - term2) == 0
 
     def test_compound_subs(self):
         grid = Grid(shape=(11,))
@@ -469,16 +471,16 @@ def test_compound_subs(self):
 
         coeffs0 = np.array([100, 100, 100])
 
-        subs = Substitutions(Coefficient(1, p, x, coeffs0))
+        subs = Substitutions(Coefficient(2, p, x, coeffs0))
 
-        eq = Eq(p.forward, (f*p).dx, coefficients=subs)
+        eq = Eq(p.forward, (f*p).dx2, coefficients=subs)
 
         term0 = f*p*100
         term1 = (f*p*100).subs(x, x-hx)
         term2 = (f*p*100).subs(x, x+hx)
 
         # `str` simply because some objects are of type EvalDerivative
-        assert str(eq.evaluate.rhs) == str(term0 + term1 + term2)
+        assert sp.simplify(eq.evaluate.rhs - (term0 + term1 + term2)) == 0
 
     def test_compound_nested_subs(self):
         grid = Grid(shape=(11, 11))
@@ -492,10 +494,10 @@ def test_compound_nested_subs(self):
         coeffs0 = np.array([100, 100, 100])
         coeffs1 = np.array([200, 200, 200])
 
-        subs = Substitutions(Coefficient(1, p, x, coeffs0),
-                             Coefficient(1, p, y, coeffs1))
+        subs = Substitutions(Coefficient(2, p, x, coeffs0),
+                             Coefficient(2, p, y, coeffs1))
 
-        eq = Eq(p.forward, (f*p.dx).dy, coefficients=subs)
+        eq = Eq(p.forward, (f*p.dx2).dy2, coefficients=subs)
 
         mul = lambda e, i: sp.Mul(f.subs(y, y+i*hy), e, 200, evaluate=False)
         term0 = mul(p*100 +
@@ -509,7 +511,7 @@ def test_compound_nested_subs(self):
                     p.subs({x: x+hx, y: y+hy})*100, 1)
 
         # `str` simply because some objects are of type EvalDerivative
-        assert str(eq.evaluate.rhs) == str(term0 + term1 + term2)
+        assert sp.simplify(eq.evaluate.rhs - (term0 + term1 + term2)) == 0
 
     def test_priority(self):
         grid = Grid(shape=(11,))