diff --git a/devito/finite_differences/__init__.py b/devito/finite_differences/__init__.py index 4a4c491572..16ba9993a3 100644 --- a/devito/finite_differences/__init__.py +++ b/devito/finite_differences/__init__.py @@ -3,5 +3,4 @@ from .finite_difference import * # noqa from .derivative import * # noqa from .tools import generate_fd_shortcuts # noqa -from .coefficients import * # noqa from .operators import * # noqa diff --git a/devito/finite_differences/coefficients.py b/devito/finite_differences/coefficients.py deleted file mode 100644 index 9320e16740..0000000000 --- a/devito/finite_differences/coefficients.py +++ /dev/null @@ -1,291 +0,0 @@ -import sympy -import numpy as np -from cached_property import cached_property - -from devito.finite_differences import generate_indices -from devito.tools import filter_ordered, as_tuple -from devito.symbolics.search import retrieve_dimensions - -__all__ = ['Coefficient', 'Substitutions', 'default_rules'] - - -class Coefficient(object): - """ - Prepare custom coefficients to pass to a Substitutions object. - - Parameters - ---------- - deriv_order : int - The order of the derivative being taken. - function : Function - The function for which the supplied coefficients - will be used. - dimension : Dimension - The dimension with respect to which the - derivative is being taken. - weights : np.ndarray - The set of finite difference weights - intended to be used in place of the standard - weights (obtained from a Taylor expansion). - - Examples - -------- - >>> import numpy as np - >>> from devito import Grid, Function, Coefficient - >>> grid = Grid(shape=(4, 4)) - >>> u = Function(name='u', grid=grid, space_order=2, coefficients='symbolic') - >>> x, y = grid.dimensions - - Now define some partial d/dx FD coefficients of the Function u: - - >>> u_x_coeffs = Coefficient(1, u, x, np.array([-0.6, 0.1, 0.6])) - - And some partial d^2/dy^2 FD coefficients: - - >>> u_y2_coeffs = Coefficient(2, u, y, np.array([0.0, 0.0, 0.0])) - - """ - - def __init__(self, deriv_order, function, dimension, weights): - - self._check_input(deriv_order, function, dimension, weights) - - # Ensure the given set of weights is the correct length - try: - wl = weights.shape[-1]-1 - except AttributeError: - wl = len(weights)-1 - if dimension.is_Time: - if wl != function.time_order: - raise ValueError("Number of FD weights provided does not " - "match the functions space_order") - elif dimension.is_Space: - if wl != function.space_order: - raise ValueError("Number of FD weights provided does not " - "match the functions space_order") - - self._deriv_order = deriv_order - self._function = function - self._dimension = dimension - self._weights = weights - - @property - def deriv_order(self): - """The derivative order.""" - return self._deriv_order - - @property - def function(self): - """The function to which the coefficients belong.""" - return self._function - - @property - def dimension(self): - """The dimension to which the coefficients will be applied.""" - return self._dimension - - @property - def index(self): - """ - The dimension to which the coefficients will be applied plus the offset - in that dimension if the function is staggered. - """ - return self._function.indices_ref[self._dimension] - - @property - def weights(self): - """The set of weights.""" - return self._weights - - def _check_input(self, deriv_order, function, dimension, weights): - if not isinstance(deriv_order, int): - raise TypeError("Derivative order must be an integer") - try: - if not function.is_Function: - raise TypeError("Object is not of type Function") - except AttributeError: - raise TypeError("Object is not of type Function") - try: - if not dimension.is_Dimension: - raise TypeError("Coefficients must be attached to a valid dimension") - except AttributeError: - raise TypeError("Coefficients must be attached to a valid dimension") - try: - weights.is_Function is True - except AttributeError: - if not isinstance(weights, np.ndarray): - raise TypeError("Weights must be of type np.ndarray or a Devito Function") - return - - -class Substitutions(object): - """ - Devito class to convert Coefficient objects into replacent rules - to be applied when constructing a Devito Eq. - - Examples - -------- - >>> from devito import Grid, TimeFunction, Coefficient - >>> grid = Grid(shape=(4, 4)) - >>> u = TimeFunction(name='u', grid=grid, space_order=2, coefficients='symbolic') - >>> x, y = grid.dimensions - - Now define some partial d/dx FD coefficients of the Function u: - - >>> u_x_coeffs = Coefficient(1, u, x, np.array([-0.6, 0.1, 0.6])) - - And now create our Substitutions object to pass to equation: - - >>> from devito import Substitutions - >>> subs = Substitutions(u_x_coeffs) - - Now create a Devito equation and pass to it 'subs' - - >>> from devito import Eq - >>> eq = Eq(u.dt+u.dx, coefficients=subs) - - When evaluated, the derivatives will use the custom coefficients. We can - check that by - - >>> eq.evaluate - 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) - - Notes - ----- - If a Function is declared with 'symbolic' coefficients and no - replacement rules for any derivative appearing in a Devito equation, - the coefficients will revert to those of the 'default' Taylor - expansion. - """ - - def __init__(self, *args): - - if any(not isinstance(arg, Coefficient) for arg in args): - raise TypeError("Non Coefficient object within input") - - self._coefficients = args - self._function_list = self.function_list - self._rules = self.rules - - @property - def coefficients(self): - """The Coefficient objects passed.""" - return self._coefficients - - @cached_property - def function_list(self): - return filter_ordered((i.function for i in self.coefficients), lambda i: i.name) - - @cached_property - def rules(self): - - def generate_subs(i): - - deriv_order = i.deriv_order - function = i.function - dim = i.dimension - index = i.index - weights = i.weights - - if isinstance(weights, np.ndarray): - fd_order = len(weights)-1 - else: - fd_order = weights.shape[-1]-1 - - subs = {} - - indices, x0 = generate_indices(function, dim, fd_order, side=None) - - # NOTE: This implementation currently assumes that indices are ordered - # according to their position in the FD stencil. This may not be the - # case in all schemes and should be changed such that the weights are - # passed as a dictionary of the form {pos: w} (or something similar). - if isinstance(weights, np.ndarray): - for j in range(len(weights)): - subs.update({function._coeff_symbol - (indices[j], deriv_order, - function, index): weights[j]}) - else: - shape = weights.shape - x = weights.dimensions - for j in range(shape[-1]): - idx = list(x) - idx[-1] = j - subs.update({function._coeff_symbol - (indices[j], deriv_order, function, index): - weights[as_tuple(idx)]}) - - return subs - - # Figure out when symbolic coefficients can be replaced - # with user provided coefficients and, if possible, generate - # replacement rules - rules = {} - for i in self.coefficients: - rules.update(generate_subs(i)) - - return rules - - -def default_rules(obj, functions): - - def generate_subs(deriv_order, function, index): - dim = retrieve_dimensions(index)[0] - - if dim.is_Time: - fd_order = function.time_order - elif dim.is_Space: - fd_order = function.space_order - else: - # Shouldn't arrive here - raise TypeError("Dimension type not recognised") - - subs = {} - - mapper = {dim: index} - - indices, x0 = generate_indices(function, dim, - fd_order, side=None, x0=mapper) - - coeffs = sympy.finite_diff_weights(deriv_order, indices, x0)[-1][-1] - - for j in range(len(coeffs)): - subs.update({function._coeff_symbol - (indices[j], deriv_order, function, index): coeffs[j]}) - - return subs - - # Determine which 'rules' are missing - sym = get_sym(functions) - terms = obj.find(sym) - args_present = filter_ordered(term.args[1:] for term in terms) - - subs = obj.substitutions - if subs: - args_provided = [(i.deriv_order, i.function, i.index) - for i in subs.coefficients] - else: - args_provided = [] - - # NOTE: Do we want to throw a warning if the same arg has - # been provided twice? - args_provided = list(set(args_provided)) - not_provided = [i for i in args_present if i not in frozenset(args_provided)] - - rules = {} - for i in not_provided: - rules = {**rules, **generate_subs(*i)} - - return rules - - -def get_sym(functions): - for f in functions: - try: - sym = f._coeff_symbol - return sym - except AttributeError: - pass - # Shouldn't arrive here - raise TypeError("Failed to retreive symbol") diff --git a/devito/finite_differences/derivative.py b/devito/finite_differences/derivative.py index 9c308da0eb..d7d0c7fe62 100644 --- a/devito/finite_differences/derivative.py +++ b/devito/finite_differences/derivative.py @@ -44,6 +44,8 @@ class Derivative(sympy.Derivative, Differentiable): x0 : dict, optional Origin (where the finite-difference is evaluated at) for the finite-difference scheme, e.g. {x: x, y: y + h_y/2}. + coeffs : iterable or Function, optional + Set of custom finite-difference weights. Examples -------- @@ -85,7 +87,7 @@ class Derivative(sympy.Derivative, Differentiable): evaluation are `x0`, `fd_order` and `side`. """ - _state = ('expr', 'dims', 'side', 'fd_order', 'transpose', '_ppsubs', 'x0') + _state = ('expr', 'dims', 'side', 'fd_order', 'transpose', '_ppsubs', 'x0', 'coeffs') _fd_priority = 3 def __new__(cls, expr, *dims, **kwargs): @@ -108,6 +110,7 @@ def __new__(cls, expr, *dims, **kwargs): obj._transpose = kwargs.get("transpose", direct) obj._ppsubs = as_tuple(frozendict(i) for i in kwargs.get("subs", [])) obj._x0 = frozendict(kwargs.get('x0', {})) + obj._coeffs = kwargs.get('coeffs', None) return obj @classmethod @@ -172,18 +175,21 @@ def _process_kwargs(cls, expr, *dims, **kwargs): for s in filter_ordered(new_dims)] return new_dims, orders, fd_orders, variable_count - def __call__(self, x0=None, fd_order=None, side=None): + def __call__(self, x0=None, coeffs=None, fd_order=None, side=None): if self.ndims == 1: _fd_order = fd_order or self._fd_order _side = side or self._side new_x0 = {self.dims[0]: x0} if x0 is not None else self.x0 - return self._new_from_self(fd_order=_fd_order, side=_side, x0=new_x0) + new_coeffs = coeffs + return self._new_from_self(fd_order=_fd_order, side=_side, x0=new_x0, + coeffs=new_coeffs) if side is not None: raise TypeError("Side only supported for first order single" "Dimension derivative such as `.dxl` or .dx(side=left)") # Cross derivative _x0 = dict(self._x0) + _coeffs = coeffs _fd_order = dict(self.fd_order._getters) try: _fd_order.update(**(fd_order or {})) @@ -193,13 +199,13 @@ def __call__(self, x0=None, fd_order=None, side=None): except AttributeError: raise TypeError("Multi-dimensional Derivative, input expected as a dict") - return self._new_from_self(fd_order=_fd_order, x0=_x0) + return self._new_from_self(fd_order=_fd_order, x0=_x0, coeffs=_coeffs) def _new_from_self(self, **kwargs): expr = kwargs.pop('expr', self.expr) _kwargs = {'deriv_order': self.deriv_order, 'fd_order': self.fd_order, 'side': self.side, 'transpose': self.transpose, 'subs': self._ppsubs, - 'x0': self.x0, 'preprocessed': True} + 'x0': self.x0, 'coeffs': self.coeffs, 'preprocessed': True} _kwargs.update(**kwargs) return Derivative(expr, *self.dims, **_kwargs) @@ -254,6 +260,10 @@ def fd_order(self): def deriv_order(self): return self._deriv_order + @property + def coeffs(self): + return self._coeffs + @property def side(self): return self._side @@ -341,13 +351,13 @@ def _eval_fd(self, expr): if self.side is not None and self.deriv_order == 1: res = first_derivative(expr, self.dims[0], self.fd_order, side=self.side, matvec=self.transpose, - x0=self.x0) + x0=self.x0, coeffs=self.coeffs) elif len(self.dims) > 1: res = cross_derivative(expr, self.dims, self.fd_order, self.deriv_order, - matvec=self.transpose, x0=self.x0) + matvec=self.transpose, x0=self.x0, coeffs=self.coeffs) else: res = generic_derivative(expr, *self.dims, self.fd_order, self.deriv_order, - matvec=self.transpose, x0=self.x0) + matvec=self.transpose, x0=self.x0, coeffs=self.coeffs) # Step 3: Evaluate remaining part of expression res = res.evaluate diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index b731d5c338..75d4141c40 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -3,7 +3,7 @@ from devito.finite_differences.differentiable import EvalDerivative from devito.finite_differences.tools import (numeric_weights, symbolic_weights, left, right, generate_indices, centered, direct, - transpose, check_input, check_symbolic) + transpose, check_input) __all__ = ['first_derivative', 'second_derivative', 'cross_derivative', 'generic_derivative', 'left', 'right', 'centered', 'transpose', @@ -15,9 +15,8 @@ @check_input -@check_symbolic def first_derivative(expr, dim, fd_order=None, side=centered, matvec=direct, - symbolic=False, x0=None): + x0=None, coeffs=None, indices=None): """ First-order derivative of a given expression. @@ -38,6 +37,10 @@ def first_derivative(expr, dim, fd_order=None, side=centered, matvec=direct, finite difference. Defaults to ``direct``. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. + coeffs : iterable or Function, optional + Set of custom finite-difference weights. + indices : iterable, optional + The set of indices at which to apply coeffs. Returns ------- @@ -85,11 +88,13 @@ def first_derivative(expr, dim, fd_order=None, side=centered, matvec=direct, 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] + ind = list(indices) if indices else 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) + # NOTE: fix up - add function case + if coeffs: + #TODO: if indices not none add consistency check + c = list(coeffs) else: c = numeric_weights(1, ind, dim) @@ -97,7 +102,6 @@ def first_derivative(expr, dim, fd_order=None, side=centered, matvec=direct, @check_input -@check_symbolic def second_derivative(expr, dim, fd_order, **kwargs): """ Second-order derivative of a given expression. @@ -115,6 +119,8 @@ def second_derivative(expr, dim, fd_order, **kwargs): Shift of the finite-difference approximation. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. + coeffs : iterable or Function, optional + Set of custom finite-difference weights. Returns ------- @@ -154,7 +160,6 @@ def second_derivative(expr, dim, fd_order, **kwargs): @check_input -@check_symbolic def cross_derivative(expr, dims, fd_order, deriv_order, **kwargs): """ Arbitrary-order cross derivative of a given expression. @@ -174,6 +179,8 @@ def cross_derivative(expr, dims, fd_order, deriv_order, **kwargs): Shift of the finite-difference approximation. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. + coeffs : iterable or Function, optional + Set of custom finite-difference weights. Returns ------- @@ -212,15 +219,15 @@ def cross_derivative(expr, dims, fd_order, deriv_order, **kwargs): """ 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) + expr = generic_derivative(expr, dim=dim, fd_order=fd, deriv_order=d, x0=x0, + coeffs=coeffs) return expr @check_input -@check_symbolic -def generic_derivative(expr, dim, fd_order, deriv_order, symbolic=False, - matvec=direct, x0=None): +def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, + x0=None, coeffs=None, indices=None): """ Arbitrary-order derivative of a given expression. @@ -239,6 +246,8 @@ def generic_derivative(expr, dim, fd_order, deriv_order, symbolic=False, Shift of the finite-difference approximation. x0 : dict, optional Origin of the finite-difference scheme as a map dim: origin_dim. + coeffs : iterable or Function, optional + Set of custom finite-difference weights. Returns ------- @@ -247,14 +256,19 @@ def generic_derivative(expr, dim, fd_order, deriv_order, symbolic=False, """ # 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: + if deriv_order == 1 and fd_order == 2 and not coeffs: fd_order = 1 # Stencil positions - indices, x0 = generate_indices(expr, dim, fd_order, x0=x0) + if indices: + indices = list(indices) + else: + 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) + # NOTE: Fix me. Add function case + # TODO: Add coeffs sanity check. Create decotrator for this + if coeffs: + c = list(coeffs) else: c = numeric_weights(deriv_order, indices, x0) diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py index e97436c482..bfe87f8c16 100644 --- a/devito/finite_differences/tools.py +++ b/devito/finite_differences/tools.py @@ -54,20 +54,6 @@ def wrapper(expr, *args, **kwargs): return wrapper -def check_symbolic(func): - @wraps(func) - def wrapper(expr, *args, **kwargs): - if expr._uses_symbolic_coefficients: - expr_dict = expr.as_coefficients_dict() - if any(len(expr_dict) > 1 for item in expr_dict): - raise NotImplementedError("Applying the chain rule to functions " - "with symbolic coefficients is not currently " - "supported") - kwargs['symbolic'] = expr._uses_symbolic_coefficients - return func(expr, *args, **kwargs) - return wrapper - - def dim_with_order(dims, orders): """ Create all possible derivative order for each dims diff --git a/devito/types/dense.py b/devito/types/dense.py index 7dc773a6da..37a61dcd2a 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -67,11 +67,6 @@ def __init_finalize__(self, *args, **kwargs): # There may or may not be a `Grid` attached to the DiscreteFunction self._grid = kwargs.get('grid') - # Symbolic (finite difference) coefficients - self._coefficients = kwargs.get('coefficients', 'standard') - if self._coefficients not in ('standard', 'symbolic'): - raise ValueError("coefficients must be `standard` or `symbolic`") - # Data-related properties and data initialization self._data = None self._first_touch = kwargs.get('first_touch', configuration['first-touch']) @@ -197,19 +192,6 @@ def grid(self): def staggered(self): return self._staggered - @property - def coefficients(self): - """Form of the coefficients of the function.""" - return self._coefficients - - @cached_property - def _coeff_symbol(self): - if self.coefficients == 'symbolic': - return sympy.Function('W') - else: - raise ValueError("Function was not declared with symbolic " - "coefficients.") - @cached_property def shape(self): """ diff --git a/devito/types/equation.py b/devito/types/equation.py index 8b927c0882..d6eeb7182c 100644 --- a/devito/types/equation.py +++ b/devito/types/equation.py @@ -4,7 +4,7 @@ from cached_property import cached_property -from devito.finite_differences import default_rules +from devito.logger import warning from devito.tools import as_tuple from devito.types.lazy import Evaluable diff --git a/tests/test_custom_coefficients.py b/tests/test_custom_coefficients.py new file mode 100644 index 0000000000..ab6865283b --- /dev/null +++ b/tests/test_custom_coefficients.py @@ -0,0 +1,73 @@ +import numpy as np +import sympy as sp +import pytest + +from devito import (Grid, Function, TimeFunction, Eq, + Dimension, solve, Operator, NODE) +from devito.finite_differences import Differentiable +from devito.tools import as_tuple + +_PRECISION = 9 + + +class TestSC(object): + """ + Class for testing custom coefficients functionality + """ + + @pytest.mark.parametrize('expr, sorder, dorder, dim, weights, expected', [ + ('u.dx', 2, 1, 0, [-0.6, 0.1, 0.6], + '0.1*u(x, y) - 0.6*u(x - h_x, y) + 0.6*u(x + h_x, y)'), + ('u.dy2', 3, 2, 1, [0.121, -0.223, 1.648, -2.904], + '1.648*u(x, y) + 0.121*u(x, y - 2*h_y) - 0.223*u(x, y - h_y) \ +- 2.904*u(x, y + h_y)')]) + def test_coefficients(self, expr, sorder, dorder, dim, weights, expected): + """Test that custom coefficients return the expected result""" + grid = Grid(shape=(10, 10)) + u = Function(name='u', grid=grid, space_order=sorder) + + eq = Eq(eval(expr)(coeffs=weights)) + assert isinstance(eq.lhs, Differentiable) + assert expected == str(eq.evaluate.lhs) + + #def test_function_coefficients(self): + #"""Test that custom function coefficients return the expected result""" + #so = 2 + #grid = Grid(shape=(4, 4)) + #f0 = TimeFunction(name='f0', grid=grid, space_order=so, coefficients='symbolic') + #f1 = TimeFunction(name='f1', grid=grid, space_order=so) + #x, y = grid.dimensions + + #s = Dimension(name='s') + #ncoeffs = so+1 + + #wshape = list(grid.shape) + #wshape.append(ncoeffs) + #wshape = as_tuple(wshape) + + #wdims = list(grid.dimensions) + #wdims.append(s) + #wdims = as_tuple(wdims) + + #w = Function(name='w', dimensions=wdims, shape=wshape) + #w.data[:, :, 0] = 0.0 + #w.data[:, :, 1] = -1.0/grid.spacing[0] + #w.data[:, :, 2] = 1.0/grid.spacing[0] + + #f_x_coeffs = Coefficient(1, f0, x, w) + + #subs = Substitutions(f_x_coeffs) + + #eq0 = Eq(f0.dt + f0.dx, 1, coefficients=subs) + #eq1 = Eq(f1.dt + f1.dx, 1) + + #stencil0 = solve(eq0.evaluate, f0.forward) + #stencil1 = solve(eq1.evaluate, f1.forward) + + #op0 = Operator(Eq(f0.forward, stencil0)) + #op1 = Operator(Eq(f1.forward, stencil1)) + + #op0(time_m=0, time_M=5, dt=1.0) + #op1(time_m=0, time_M=5, dt=1.0) + + #assert np.all(np.isclose(f0.data[:] - f1.data[:], 0.0, atol=1e-5, rtol=0))