Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

api: more robust processing of custom weight with spacing #2489

Merged
merged 2 commits into from
Nov 17, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
misc: make builtin grid unit spacing
  • Loading branch information
mloubout committed Nov 17, 2024
commit 6f00de218a40cca1fbc5268bf4b531cf4dd0643c
4 changes: 3 additions & 1 deletion devito/builtins/initializers.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
14 changes: 5 additions & 9 deletions devito/finite_differences/finite_difference.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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
mloubout marked this conversation as resolved.
Show resolved Hide resolved

return deriv
12 changes: 7 additions & 5 deletions devito/finite_differences/tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion examples/seismic/tutorials/07_DRP_schemes.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
}
],
Expand Down
Loading