Merge pull request #2489 from devitocodes/custom-coeff-spacing

api: more robust processing of custom weight with spacing

devito/builtins/initializers.py

@@ -179,7 +179,9 @@ def fset(f, g):
     # Create the padded grid:
     objective_domain = ObjectiveDomain(lw)
     shape_padded = tuple([np.array(s) + 2*l for s, l in zip(shape, lw)])
-    grid = dv.Grid(shape=shape_padded, subdomains=objective_domain)
+    extent_padded = tuple([s-1 for s in shape_padded])
+    grid = dv.Grid(shape=shape_padded, subdomains=objective_domain,
+                   extent=extent_padded)
 
     f_c = dv.Function(name='f_c', grid=grid, space_order=2*max(lw), dtype=dtype)
     f_o = dv.Function(name='f_o', grid=grid, dtype=dtype)

devito/finite_differences/finite_difference.py

@@ -151,7 +151,7 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici
     expand = True if dim.is_Time else expand
 
     # The stencil indices
-    nweights, wdim = process_weights(weights, expr)
+    nweights, wdim, scale = process_weights(weights, expr, dim)
     indices, x0 = generate_indices(expr, dim, fd_order, side=side, matvec=matvec,
                                    x0=x0, nweights=nweights)
     # Finite difference weights corresponding to the indices. Computed via the
@@ -208,4 +208,7 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici
 
         deriv = EvalDerivative(*terms, base=expr)
 
+    if scale:
+        deriv = dim.spacing**(-deriv_order) * deriv
+
     return deriv

devito/finite_differences/tools.py

@@ -321,16 +321,18 @@ def make_shift_x0(shift, ndim):
                      "None, float or tuple with shape equal to %s" % (ndim,))
 
 
-def process_weights(weights, expr):
+def process_weights(weights, expr, dim):
     if weights is None:
-        return 0, None
+        return 0, None, False
     elif isinstance(weights, Function):
         if len(weights.dimensions) == 1:
-            return weights.shape[0], weights.dimensions[0]
+            return weights.shape[0], weights.dimensions[0], False
         wdim = {d for d in weights.dimensions if d not in expr.dimensions}
         assert len(wdim) == 1
         wdim = wdim.pop()
         shape = weights.shape
-        return shape[weights.dimensions.index(wdim)], wdim
+        return shape[weights.dimensions.index(wdim)], wdim, False
     else:
-        return len(list(weights)), None
+        # Adimensional weight from custom coeffs need to be multiplied by h^order
+        scale = not all(sympify(w).has(dim.spacing) for w in weights if w != 0)
+        return len(list(weights)), None, scale

examples/seismic/tutorials/07_DRP_schemes.ipynb

@@ -83,7 +83,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Eq(-u(t, x, y)/dt + u(t + dt, x, y)/dt + 0.1*u(t, x, y) - 0.6*u(t, x - h_x, y) + 0.6*u(t, x + h_x, y), 0)\n"
+      "Eq(-u(t, x, y)/dt + u(t + dt, x, y)/dt + (0.1*u(t, x, y) - 0.6*u(t, x - h_x, y) + 0.6*u(t, x + h_x, y))/h_x, 0)\n"
      ]
     }
    ],

tests/test_symbolic_coefficients.py

@@ -40,7 +40,8 @@ def test_coefficients(self, expr, sorder, dorder, dim, weights, expected):
         assert np.all(deriv.weights == weights)
 
         assert isinstance(eq.lhs, Differentiable)
-        assert sp.simplify(eval(expected).evalf(_PRECISION) - eq.evaluate.lhs) == 0
+        s = dim.spacing**(-dorder)
+        assert sp.simplify(eval(expected).evalf(_PRECISION) * s - eq.evaluate.lhs) == 0
 
     def test_function_coefficients(self):
         """Test that custom function coefficients return the expected result"""
@@ -118,7 +119,8 @@ def test_coefficients_w_xreplace(self):
         eq = Eq(u.dx2(weights=weights) + c)
         eq = eq.xreplace({c: 2})
 
-        expected = 0.1*u - 0.6*u._subs(x, x - h_x) + 0.6*u.subs(x, x + h_x) + 2
+        s = x.spacing**(-2)
+        expected = (0.1*u - 0.6*u._subs(x, x - h_x) + 0.6*u.subs(x, x + h_x)) * s + 2
 
         assert sp.simplify(expected.evalf(_PRECISION) - eq.evaluate.lhs) == 0
 
@@ -200,8 +202,8 @@ def test_staggered_equation(self):
 
         eq_f = Eq(f, f.dx2(weights=weights))
 
-        expected = 'Eq(f(x + h_x/2), 1.0*f(x - h_x/2) - 2.0*f(x + h_x/2)'\
-            ' + 1.0*f(x + 3*h_x/2))'
+        expected = 'Eq(f(x + h_x/2), (1.0*f(x - h_x/2) - 2.0*f(x + h_x/2)'\
+            ' + 1.0*f(x + 3*h_x/2))/h_x**2)'
         assert(str(eq_f.evaluate) == expected)
 
     @pytest.mark.parametrize('stagger', [True, False])
@@ -241,7 +243,8 @@ def test_nested_subs(self):
 
         eq = Eq(p.forward, p.dx2(weights=coeffs0).dy2(weights=coeffs1))
 
-        mul = lambda e: sp.Mul(e, 200, evaluate=False)
+        s = x.spacing**(-2) * y.spacing**(-2)
+        mul = lambda e: sp.Mul(e, 200 * s, evaluate=False)
         term0 = mul(p*100 +
                     p.subs(x, x-hx)*100 +
                     p.subs(x, x+hx)*100)
@@ -270,9 +273,10 @@ def test_compound_subs(self):
         term0 = f*p*100
         term1 = (f*p*100).subs(x, x-hx)
         term2 = (f*p*100).subs(x, x+hx)
+        s = x.spacing**(-2)
 
         # `str` simply because some objects are of type EvalDerivative
-        assert sp.simplify(eq.evaluate.rhs - (term0 + term1 + term2)) == 0
+        assert sp.simplify(eq.evaluate.rhs - (term0 + term1 + term2) * s) == 0
 
     def test_compound_nested_subs(self):
         grid = Grid(shape=(11, 11))
@@ -298,8 +302,9 @@ def test_compound_nested_subs(self):
                     p.subs({x: x-hx, y: y+hy})*100 +
                     p.subs({x: x+hx, y: y+hy})*100, 1)
 
+        s = x.spacing**(-2) * y.spacing**(-2)
         # `str` simply because some objects are of type EvalDerivative
-        assert sp.simplify(eq.evaluate.rhs - (term0 + term1 + term2)) == 0
+        assert sp.simplify(eq.evaluate.rhs - (term0 + term1 + term2) * s) == 0
 
     def test_operators(self):
         grid = Grid(shape=(11, 11))
@@ -324,3 +329,18 @@ def test_operators(self):
         expr3 = laplace(f, w=coeffs0)
         assert expr3 == f.dx2(weights=coeffs0) + f.dy2(weights=coeffs0)
         assert list(expr3.args[0].weights) == coeffs0
+
+    def test_spacing(self):
+        grid = Grid(shape=(11, 11))
+        x, _ = grid.dimensions
+        s = x.spacing
+
+        f = Function(name='f', grid=grid, space_order=2)
+
+        coeffs0 = [100, 100, 100]
+        coeffs1 = [100/s, 100/s, 100/s]
+
+        df = f.dx(weights=coeffs0)
+        df_s = f.dx(weights=coeffs1)
+
+        assert sp.simplify(df_s.evaluate - df.evaluate) == 0