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" ] } ],