From 7c3051a0453fd683b51bba409788009dbf615ca5 Mon Sep 17 00:00:00 2001 From: rhodrin Date: Wed, 24 Mar 2021 21:53:38 +0000 Subject: [PATCH] fd: Start overhauling custom coefficients code. --- devito/finite_differences/__init__.py | 1 - devito/finite_differences/coefficients.py | 291 --------------- devito/finite_differences/derivative.py | 26 +- .../finite_differences/finite_difference.py | 38 +- devito/finite_differences/tools.py | 14 - devito/types/dense.py | 18 - devito/types/equation.py | 1 - tests/test_symbolic_coefficients.py | 347 ------------------ 8 files changed, 41 insertions(+), 695 deletions(-) delete mode 100644 devito/finite_differences/coefficients.py delete mode 100644 tests/test_symbolic_coefficients.py diff --git a/devito/finite_differences/__init__.py b/devito/finite_differences/__init__.py index 8369da30a3..d0a38dc353 100644 --- a/devito/finite_differences/__init__.py +++ b/devito/finite_differences/__init__.py @@ -2,5 +2,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 4f7e48cf3a..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(0.1*u(t, x, y) - 0.6*u(t, x - h_x, y) + 0.6*u(t, x + h_x, y) \ -- u(t, x, y)/dt + u(t + dt, x, y)/dt, 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 6d62212cd1..44d4536b3b 100644 --- a/devito/finite_differences/derivative.py +++ b/devito/finite_differences/derivative.py @@ -43,6 +43,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 -------- @@ -84,7 +86,7 @@ class Derivative(sympy.Derivative, Differentiable): evaluation are `x0`, `fd_order` and `side`. """ - _state = ('expr', 'dims', 'side', 'fd_order', 'transpose', '_subs', 'x0') + _state = ('expr', 'dims', 'side', 'fd_order', 'transpose', '_subs', 'x0', 'coeffs') _fd_priority = 3 def __new__(cls, expr, *dims, **kwargs): @@ -107,6 +109,7 @@ def __new__(cls, expr, *dims, **kwargs): obj._transpose = kwargs.get("transpose", direct) obj._subs = 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 @@ -171,18 +174,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 {})) @@ -192,12 +198,12 @@ 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): _kwargs = {'deriv_order': self.deriv_order, 'fd_order': self.fd_order, 'side': self.side, 'transpose': self.transpose, 'subs': self._subs, - 'x0': self.x0, 'preprocessed': True} + 'x0': self.x0, 'coeffs': self.coeffs, 'preprocessed': True} expr = kwargs.pop('expr', self.expr) _kwargs.update(**kwargs) return Derivative(expr, *self.dims, **_kwargs) @@ -250,6 +256,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 @@ -340,13 +350,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 a472ae6e02..1539a1c762 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -1,9 +1,9 @@ from sympy import sympify from devito.finite_differences.differentiable import Add -from devito.finite_differences.tools import (numeric_weights, symbolic_weights, left, +from devito.finite_differences.tools import (numeric_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): """ First-order derivative of a given expression. @@ -38,6 +37,8 @@ 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. Returns ------- @@ -88,8 +89,9 @@ def first_derivative(expr, dim, fd_order=None, side=centered, matvec=direct, 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) + # NOTE: fix up - add function case + if coeffs: + c = list(coeffs) else: c = numeric_weights(1, ind, dim) @@ -97,7 +99,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 +116,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 +157,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 +176,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 +216,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): """ Arbitrary-order derivative of a given expression. @@ -239,6 +243,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 +253,16 @@ 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) # 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 ba3095f17b..945074f692 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']) @@ -196,19 +191,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 4b86a59ec9..2b380f483f 100644 --- a/devito/types/equation.py +++ b/devito/types/equation.py @@ -5,7 +5,6 @@ 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_symbolic_coefficients.py b/tests/test_symbolic_coefficients.py deleted file mode 100644 index fa8a6d2617..0000000000 --- a/tests/test_symbolic_coefficients.py +++ /dev/null @@ -1,347 +0,0 @@ -import numpy as np -import sympy as sp -import pytest - -from devito import (Grid, Function, TimeFunction, Eq, Coefficient, Substitutions, - Dimension, solve, Operator, NODE) -from devito.finite_differences import Differentiable -from devito.tools import as_tuple - -_PRECISION = 9 - - -class TestSC(object): - """ - Class for testing symbolic coefficients functionality - """ - - @pytest.mark.parametrize('order', [1, 2, 6]) - @pytest.mark.parametrize('stagger', [True, False]) - def test_default_rules(self, order, stagger): - """ - Test that the default replacement rules return the same - as standard FD. - """ - grid = Grid(shape=(20, 20)) - if stagger: - staggered = grid.dimensions[0] - else: - staggered = None - u0 = TimeFunction(name='u', grid=grid, time_order=order, space_order=order, - staggered=staggered) - u1 = TimeFunction(name='u', grid=grid, time_order=order, space_order=order, - staggered=staggered, coefficients='symbolic') - eq0 = Eq(-u0.dx+u0.dt) - eq1 = Eq(u1.dt-u1.dx) - assert(eq0.evalf(_PRECISION).__repr__() == eq1.evalf(_PRECISION).__repr__()) - - @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, coefficients='symbolic') - x = grid.dimensions - - order = dorder - dim = x[dim] - weights = np.array(weights) - - coeffs = Coefficient(order, u, dim, weights) - - eq = Eq(eval(expr), coefficients=Substitutions(coeffs)) - 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)) - - def test_coefficients_w_xreplace(self): - """Test custom coefficients with an xreplace before they are applied""" - grid = Grid(shape=(4, 4)) - u = Function(name='u', grid=grid, space_order=2, coefficients='symbolic') - x = grid.dimensions[0] - - dorder = 1 - weights = np.array([-0.6, 0.1, 0.6]) - - coeffs = Coefficient(dorder, u, x, weights) - - c = sp.Symbol('c') - - eq = Eq(u.dx+c, coefficients=Substitutions(coeffs)) - eq = eq.xreplace({c: 2}) - - expected = '0.1*u(x, y) - 0.6*u(x - h_x, y) + 0.6*u(x + h_x, y) + 2' - - assert expected == str(eq.evaluate.lhs) - - @pytest.mark.parametrize('order', [1, 2, 6, 8]) - @pytest.mark.parametrize('extent', [1., 10., 100.]) - @pytest.mark.parametrize('conf', [{'l': 'NODE', 'r1': 'x', 'r2': None}, - {'l': 'NODE', 'r1': 'y', 'r2': None}, - {'l': 'NODE', 'r1': '(x, y)', 'r2': None}, - {'l': 'x', 'r1': 'NODE', 'r2': None}, - {'l': 'y', 'r1': 'NODE', 'r2': None}, - {'l': '(x, y)', 'r1': 'NODE', 'r2': None}, - {'l': 'NODE', 'r1': 'x', 'r2': 'y'}]) - def test_default_rules_vs_standard(self, order, extent, conf): - """ - Test that equations containing default symbolic coefficients evaluate to - the same expressions as standard coefficients for the same function. - """ - def function_setup(name, grid, order, stagger): - x, y = grid.dimensions - if stagger == 'NODE': - staggered = NODE - elif stagger == 'x': - staggered = x - elif stagger == 'y': - staggered = y - elif stagger == '(x, y)': - staggered = (x, y) - else: - raise ValueError("Invalid stagger in configuration") - - f_std = Function(name=name, grid=grid, space_order=order, - staggered=staggered) - f_sym = Function(name=name, grid=grid, space_order=order, - staggered=staggered, coefficients='symbolic') - - return f_std, f_sym - - def get_eq(u, a, b, conf): - if conf['l'] == 'x' or conf['r1'] == 'x': - a_deriv = a.dx - elif conf['l'] == 'y' or conf['r1'] == 'y': - a_deriv = a.dy - elif conf['l'] == '(x, y)' or conf['r1'] == '(x, y)': - a_deriv = a.dx + a.dy - else: - raise ValueError("Invalid configuration") - - if conf['r2'] == 'y': - b_deriv = b.dy - elif conf['r2'] == '(x, y)': - b_deriv = b.dx + b.dy - elif conf['r2'] is None: - b_deriv = 0. - else: - raise ValueError("Invalid configuration") - - return Eq(u, a_deriv + b_deriv) - - grid = Grid(shape=(11, 11), extent=(extent, extent)) - - # Set up functions as specified - u_std, u_sym = function_setup('u', grid, order, conf['l']) - a_std, a_sym = function_setup('a', grid, order, conf['r1']) - a_std.data[::2, ::2] = 1. - a_sym.data[::2, ::2] = 1. - if conf['r2'] is not None: - b_std, b_sym = function_setup('b', grid, order, conf['r2']) - b_std.data[::2, ::2] = 1. - b_sym.data[::2, ::2] = 1. - else: - b_std, b_sym = 0., 0. - - eq_std = get_eq(u_std, a_std, b_std, conf) - eq_sym = get_eq(u_sym, a_sym, b_sym, conf) - - evaluated_std = eq_std.evaluate.evalf(_PRECISION) - evaluated_sym = eq_sym.evaluate.evalf(_PRECISION) - - assert str(evaluated_std) == str(evaluated_sym) - - Operator(eq_std)() - Operator(eq_sym)() - - assert np.all(np.isclose(u_std.data - u_sym.data, 0.0, atol=1e-5, rtol=0)) - - @pytest.mark.parametrize('so, expected', [(1, '-f(x)/h_x + f(x + h_x)/h_x'), - (4, '-9*f(x)/(8*h_x) + f(x - h_x)' - '/(24*h_x) + 9*f(x + h_x)/(8*h_x)' - ' - f(x + 2*h_x)/(24*h_x)'), - (6, '-75*f(x)/(64*h_x)' - ' - 3*f(x - 2*h_x)/(640*h_x)' - ' + 25*f(x - h_x)/(384*h_x)' - ' + 75*f(x + h_x)/(64*h_x)' - ' - 25*f(x + 2*h_x)/(384*h_x)' - ' + 3*f(x + 3*h_x)/(640*h_x)')]) - def test_default_rules_vs_string(self, so, expected): - """ - Test that default_rules generates correct symbolic expressions when used - with staggered grids. - """ - grid = Grid(shape=(11,), extent=(10.,)) - x = grid.dimensions[0] - f = Function(name='f', grid=grid, space_order=so, staggered=NODE, - coefficients='symbolic') - g = Function(name='g', grid=grid, space_order=so, staggered=x, - coefficients='symbolic') - eq = Eq(g, f.dx) - assert str(eq.evaluate.rhs) == expected - - @pytest.mark.parametrize('so', [1, 4, 6]) - @pytest.mark.parametrize('offset', [1, -1]) - def test_default_rules_deriv_offset(self, so, offset): - """ - Test that default_rules generates correct derivatives when derivatives - are evaluated at offset x0. - """ - grid = Grid(shape=(11,), extent=(10.,)) - x = grid.dimensions[0] - h_x = x.spacing - f_std = Function(name='f', grid=grid, space_order=so) - f_sym = Function(name='f', grid=grid, space_order=so, coefficients='symbolic') - g = Function(name='g', grid=grid, space_order=so) - - eval_std = str(Eq(g, f_std.dx(x0=x+offset*h_x/2)).evaluate.evalf(_PRECISION)) - eval_sym = str(Eq(g, f_sym.dx(x0=x+offset*h_x/2)).evaluate.evalf(_PRECISION)) - assert eval_std == eval_sym - - @pytest.mark.parametrize('order', [2, 4, 6]) - def test_staggered_array(self, order): - """Test custom coefficients provided as an array on a staggered grid""" - grid = Grid(shape=(11,), extent=(10.,)) - x = grid.dimensions[0] - - f = Function(name='f', grid=grid, space_order=order, - coefficients='symbolic') - g = Function(name='g', grid=grid, space_order=order, - coefficients='symbolic', staggered=x) - f.data[::2] = 1 - g.data[::2] = 1 - - weights = np.ones(order+1)/grid.spacing[0]**2 - coeffs_f = Coefficient(2, f, x, weights) - coeffs_g = Coefficient(2, g, x, weights) - - eq_f = Eq(f, f.dx2, coefficients=Substitutions(coeffs_f)) - eq_g = Eq(g, g.dx2, coefficients=Substitutions(coeffs_g)) - - Operator([eq_f, eq_g])() - - assert np.allclose(f.data, g.data, atol=1e-7) - - @pytest.mark.parametrize('order', [2, 4, 6]) - def test_staggered_function(self, order): - """Test custom function coefficients on a staggered grid""" - grid = Grid(shape=(11,), extent=(10.,)) - x = grid.dimensions[0] - - f = Function(name='f', grid=grid, space_order=order, - coefficients='symbolic') - g = Function(name='g', grid=grid, space_order=order, - coefficients='symbolic', staggered=x) - f.data[::2] = 1 - g.data[::2] = 1 - - s = Dimension(name='s') - ncoeffs = order+1 - - wshape = grid.shape + (ncoeffs,) - wdims = grid.dimensions + (s,) - - w = Function(name='w', dimensions=wdims, shape=wshape) - w.data[:] = 1.0/grid.spacing[0]**2 - - coeffs_f = Coefficient(2, f, x, w) - coeffs_g = Coefficient(2, g, x, w) - - eq_f = Eq(f, f.dx2, coefficients=Substitutions(coeffs_f)) - eq_g = Eq(g, g.dx2, coefficients=Substitutions(coeffs_g)) - - Operator([eq_f, eq_g])() - - assert np.allclose(f.data, g.data, atol=1e-7) - - def test_staggered_equation(self): - """ - Check that expressions with substitutions are consistent with - those without - """ - grid = Grid(shape=(11,), extent=(10.,)) - x = grid.dimensions[0] - - f = Function(name='f', grid=grid, space_order=2, - coefficients='symbolic', staggered=x) - - weights = np.array([1, -2, 1])/grid.spacing[0]**2 - coeffs_f = Coefficient(2, f, x, weights) - - eq_f = Eq(f, 1.0*f.dx2, coefficients=Substitutions(coeffs_f)) - - 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))' - assert(str(eq_f.evaluate) == expected) - - @pytest.mark.parametrize('stagger', [True, False]) - def test_with_timefunction(self, stagger): - """Check compatibility of custom coefficients and TimeFunctions""" - grid = Grid(shape=(11,), extent=(10.,)) - x = grid.dimensions[0] - if stagger: - staggered = x - else: - staggered = None - - f = TimeFunction(name='f', grid=grid, space_order=2, staggered=staggered) - g = TimeFunction(name='g', grid=grid, space_order=2, staggered=staggered, - coefficients='symbolic') - - f.data[:, ::2] = 1 - g.data[:, ::2] = 1 - - weights = np.array([-1, 2, -1])/grid.spacing[0]**2 - coeffs = Coefficient(2, g, x, weights) - - eq_f = Eq(f.forward, f.dx2) - eq_g = Eq(g.forward, g.dx2, coefficients=Substitutions(coeffs)) - - Operator([eq_f, eq_g])(t_m=0, t_M=1) - - assert np.allclose(f.data[-1], -g.data[-1], atol=1e-7)