Source code for devito.finite_differences.finite_difference

from sympy import finite_diff_weights

from devito.finite_differences.tools import (symbolic_weights, left, right,
                                             generate_indices, centered, check_input,
                                             check_symbolic, direct, transpose)

__all__ = ['first_derivative', 'second_derivative', 'cross_derivative',
           'generic_derivative', 'left', 'right', 'centered', 'transpose',
           'generate_indices']

# Number of digits for FD coefficients to avoid roundup errors and non-deterministic
# code generation
_PRECISION = 9


[docs]@check_input @check_symbolic def first_derivative(expr, dim, fd_order=None, side=centered, matvec=direct, symbolic=False, x0=None): """ First-order derivative of a given expression. Parameters ---------- expr : expr-like Expression for which the first-order derivative is produced. dim : Dimension The Dimension w.r.t. which to differentiate. fd_order : int, optional Coefficient discretization order. Note: this impacts the width of the resulting stencil. Defaults to ``expr.space_order`` side : Side, optional Side of the finite difference location, centered (at x), left (at x - 1) or right (at x +1). Defaults to ``centered``. matvec : Transpose, optional Forward (matvec=direct) or transpose (matvec=transpose) mode of the finite difference. Defaults to ``direct``. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. Returns ------- expr-like First-order derivative of ``expr``. Examples -------- >>> from devito import Function, Grid, first_derivative, transpose >>> grid = Grid(shape=(4, 4)) >>> x, _ = grid.dimensions >>> f = Function(name='f', grid=grid) >>> g = Function(name='g', grid=grid) >>> first_derivative(f*g, dim=x) -f(x, y)*g(x, y)/h_x + f(x + h_x, y)*g(x + h_x, y)/h_x Semantically, this is equivalent to >>> (f*g).dx Derivative(f(x, y)*g(x, y), x) The only difference is that in the latter case derivatives remain unevaluated. The expanded form is obtained via ``evaluate`` >>> (f*g).dx.evaluate -f(x, y)*g(x, y)/h_x + f(x + h_x, y)*g(x + h_x, y)/h_x For the adjoint mode of the first derivative, pass ``matvec=transpose`` >>> g = Function(name='g', grid=grid) >>> first_derivative(f*g, dim=x, matvec=transpose) -f(x, y)*g(x, y)/h_x + f(x - h_x, y)*g(x - h_x, y)/h_x This is also accessible via the .T shortcut >>> (f*g).dx.T.evaluate -f(x, y)*g(x, y)/h_x + f(x - h_x, y)*g(x - h_x, y)/h_x Finally the x0 argument allows to choose the origin of the finite-difference >>> first_derivative(f, dim=x, x0={x: 1}) -f(1, y)/h_x + f(h_x + 1, y)/h_x """ side = side order = fd_order or expr.space_order # Stencil positions for non-symmetric cross-derivatives with symmetric averaging ind = generate_indices(expr, dim, order, side=side, x0=x0)[0] # Finite difference weights from Taylor approximation with these positions if symbolic: c = symbolic_weights(expr, 1, ind, dim) else: c = finite_diff_weights(1, ind, dim)[-1][-1] return indices_weights_to_fd(expr, dim, ind, c, matvec=matvec.val)
[docs]@check_input @check_symbolic def second_derivative(expr, dim, fd_order, **kwargs): """ Second-order derivative of a given expression. Parameters ---------- expr : expr-like Expression for which the derivative is produced. dim : Dimension The Dimension w.r.t. which to differentiate. fd_order : int Coefficient discretization order. Note: this impacts the width of the resulting stencil. stagger : Side, optional Shift of the finite-difference approximation. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. Returns ------- expr-like Second-order derivative of ``expr``. Examples -------- >>> from devito import Function, Grid, second_derivative >>> grid = Grid(shape=(4, 4)) >>> x, _ = grid.dimensions >>> f = Function(name='f', grid=grid, space_order=2) >>> g = Function(name='g', grid=grid, space_order=2) >>> second_derivative(f*g, dim=x, fd_order=2) -2.0*f(x, y)*g(x, y)/h_x**2 + f(x - h_x, y)*g(x - h_x, y)/h_x**2 +\ f(x + h_x, y)*g(x + h_x, y)/h_x**2 Semantically, this is equivalent to >>> (f*g).dx2 Derivative(f(x, y)*g(x, y), (x, 2)) The only difference is that in the latter case derivatives remain unevaluated. The expanded form is obtained via ``evaluate`` >>> (f*g).dx2.evaluate -2.0*f(x, y)*g(x, y)/h_x**2 + f(x - h_x, y)*g(x - h_x, y)/h_x**2 +\ f(x + h_x, y)*g(x + h_x, y)/h_x**2 Finally the x0 argument allows to choose the origin of the finite-difference >>> second_derivative(f, dim=x, fd_order=2, x0={x: 1}) -2.0*f(1, y)/h_x**2 + f(1 - h_x, y)/h_x**2 + f(h_x + 1, y)/h_x**2 """ return generic_derivative(expr, dim, fd_order, 2, **kwargs)
[docs]@check_input @check_symbolic def cross_derivative(expr, dims, fd_order, deriv_order, **kwargs): """ Arbitrary-order cross derivative of a given expression. Parameters ---------- expr : expr-like Expression for which the cross derivative is produced. dims : tuple of Dimension Dimensions w.r.t. which to differentiate. fd_order : tuple of ints Coefficient discretization order. Note: this impacts the width of the resulting stencil. deriv_order : tuple of ints Derivative order, e.g. 2 for a second-order derivative. stagger : tuple of Side, optional Shift of the finite-difference approximation. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. Returns ------- expr-like Cross-derivative of ``expr``. Examples -------- >>> from devito import Function, Grid, second_derivative >>> grid = Grid(shape=(4, 4)) >>> x, y = grid.dimensions >>> f = Function(name='f', grid=grid, space_order=2) >>> g = Function(name='g', grid=grid, space_order=2) >>> cross_derivative(f*g, dims=(x, y), fd_order=(2, 2), deriv_order=(1, 1)) -(-f(x, y)*g(x, y)/h_x + f(x + h_x, y)*g(x + h_x, y)/h_x)/h_y +\ (-f(x, y + h_y)*g(x, y + h_y)/h_x + f(x + h_x, y + h_y)*g(x + h_x, y + h_y)/h_x)/h_y Semantically, this is equivalent to >>> (f*g).dxdy Derivative(f(x, y)*g(x, y), x, y) The only difference is that in the latter case derivatives remain unevaluated. The expanded form is obtained via ``evaluate`` >>> (f*g).dxdy.evaluate -(-f(x, y)*g(x, y)/h_x + f(x + h_x, y)*g(x + h_x, y)/h_x)/h_y +\ (-f(x, y + h_y)*g(x, y + h_y)/h_x + f(x + h_x, y + h_y)*g(x + h_x, y + h_y)/h_x)/h_y Finally the x0 argument allows to choose the origin of the finite-difference >>> cross_derivative(f*g, dims=(x, y), fd_order=(2, 2), deriv_order=(1, 1), \ x0={x: 1, y: 2}) -(-f(1, 2)*g(1, 2)/h_x + f(h_x + 1, 2)*g(h_x + 1, 2)/h_x)/h_y +\ (-f(1, h_y + 2)*g(1, h_y + 2)/h_x + f(h_x + 1, h_y + 2)*g(h_x + 1, h_y + 2)/h_x)/h_y """ x0 = kwargs.get('x0', {}) for d, fd, dim in zip(deriv_order, fd_order, dims): expr = generic_derivative(expr, dim=dim, fd_order=fd, deriv_order=d, x0=x0) return expr
[docs]@check_input @check_symbolic def generic_derivative(expr, dim, fd_order, deriv_order, symbolic=False, matvec=direct, x0=None): """ Arbitrary-order derivative of a given expression. Parameters ---------- expr : expr-like Expression for which the derivative is produced. dim : Dimension The Dimension w.r.t. which to differentiate. fd_order : int Coefficient discretization order. Note: this impacts the width of the resulting stencil. deriv_order : int Derivative order, e.g. 2 for a second-order derivative. stagger : Side, optional Shift of the finite-difference approximation. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. Returns ------- expr-like ``deriv-order`` derivative of ``expr``. """ # First order derivative with 2nd order FD is highly non-recommended so taking # first order fd that is a lot better if deriv_order == 1 and fd_order == 2 and not symbolic: fd_order = 1 # Stencil positions indices, x0 = generate_indices(expr, dim, fd_order, x0=x0) # Finite difference weights from Taylor approximation with these positions if symbolic: c = symbolic_weights(expr, deriv_order, indices, x0) else: c = finite_diff_weights(deriv_order, indices, x0)[-1][-1] return indices_weights_to_fd(expr, dim, indices, c, matvec=matvec.val)
def indices_weights_to_fd(expr, dim, inds, weights, matvec=1): """Expression from lists of indices and weights.""" diff = dim.spacing deriv = 0 all_dims = tuple(set((expr.indices_ref[dim],) + tuple(expr.indices_ref[dim] for i in expr.dimensions if i.root is dim))) d0 = ([d for d in expr.dimensions if d.root is dim] or [dim])[0] # Loop through weights for i, c in zip(inds, weights): try: iloc = i.xreplace({dim: d0, diff: matvec*diff}) except AttributeError: iloc = i subs = dict((d, iloc) for d in all_dims) deriv += expr.subs(subs) * c return deriv.evalf(_PRECISION)