From 3c8edd93ac517ae9230f396b55812618a26d97c5 Mon Sep 17 00:00:00 2001 From: mloubout Date: Sat, 16 Nov 2024 23:46:40 -0500 Subject: [PATCH 1/2] api: more robust processing of custom weight with spacing --- .../finite_differences/finite_difference.py | 9 ++++- tests/test_symbolic_coefficients.py | 34 +++++++++++++++---- 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index ad98f1818e..237f7438d1 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -85,7 +85,7 @@ def cross_derivative(expr, dims, fd_order, deriv_order, x0=None, side=None, **kw return expr -@check_input +# @check_input def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None, coefficients='taylor', expand=True, weights=None, side=None): """ @@ -164,6 +164,11 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici # Enforce fixed precision FD coefficients to avoid variations in results weights = [sympify(w).evalf(_PRECISION) for w in weights] + # Adimensional weight from custom coeffs need to multiply by h^order + scale = None + if wdim is None and not weights[0].has(dim.spacing): + scale = dim.spacing**(-deriv_order) + # Transpose the FD, if necessary if matvec == transpose: weights = weights[::-1] @@ -208,4 +213,6 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici deriv = EvalDerivative(*terms, base=expr) + if scale is not None: + deriv = scale * deriv return deriv diff --git a/tests/test_symbolic_coefficients.py b/tests/test_symbolic_coefficients.py index 5bebb1f240..872995ff68 100644 --- a/tests/test_symbolic_coefficients.py +++ b/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 From 6f00de218a40cca1fbc5268bf4b531cf4dd0643c Mon Sep 17 00:00:00 2001 From: mloubout Date: Sun, 17 Nov 2024 00:04:31 -0500 Subject: [PATCH 2/2] misc: make builtin grid unit spacing --- devito/builtins/initializers.py | 4 +++- devito/finite_differences/finite_difference.py | 14 +++++--------- devito/finite_differences/tools.py | 12 +++++++----- examples/seismic/tutorials/07_DRP_schemes.ipynb | 2 +- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/devito/builtins/initializers.py b/devito/builtins/initializers.py index aa7869e3f7..dbe271eac1 100644 --- a/devito/builtins/initializers.py +++ b/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) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index 237f7438d1..45f4060715 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -85,7 +85,7 @@ def cross_derivative(expr, dims, fd_order, deriv_order, x0=None, side=None, **kw return expr -# @check_input +@check_input def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None, coefficients='taylor', expand=True, weights=None, side=None): """ @@ -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 @@ -164,11 +164,6 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici # Enforce fixed precision FD coefficients to avoid variations in results weights = [sympify(w).evalf(_PRECISION) for w in weights] - # Adimensional weight from custom coeffs need to multiply by h^order - scale = None - if wdim is None and not weights[0].has(dim.spacing): - scale = dim.spacing**(-deriv_order) - # Transpose the FD, if necessary if matvec == transpose: weights = weights[::-1] @@ -213,6 +208,7 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici deriv = EvalDerivative(*terms, base=expr) - if scale is not None: - deriv = scale * deriv + if scale: + deriv = dim.spacing**(-deriv_order) * deriv + return deriv diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py index 8cb6ac608b..d153d082bd 100644 --- a/devito/finite_differences/tools.py +++ b/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 diff --git a/examples/seismic/tutorials/07_DRP_schemes.ipynb b/examples/seismic/tutorials/07_DRP_schemes.ipynb index e1f9b1f4d9..ceddb18547 100644 --- a/examples/seismic/tutorials/07_DRP_schemes.ipynb +++ b/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" ] } ],