From 5e60d919780872e3b4811b72f9e22c0492ca95b2 Mon Sep 17 00:00:00 2001 From: mloubout Date: Thu, 12 Oct 2023 11:19:53 -0400 Subject: [PATCH] api: enforce interpolation radius to be smaller than any input space order --- devito/operations/interpolators.py | 17 +++++++++++- devito/types/sparse.py | 4 +++ tests/test_dle.py | 2 +- tests/test_gpu_common.py | 2 +- tests/test_interpolation.py | 44 ++++++++++++++++-------------- tests/test_mpi.py | 6 ++-- tests/test_operator.py | 10 +++---- 7 files changed, 53 insertions(+), 32 deletions(-) diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index 3d2dcb7466..c932802f4f 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -1,11 +1,12 @@ from abc import ABC, abstractmethod +from functools import wraps import sympy from cached_property import cached_property from devito.finite_differences.differentiable import Mul from devito.finite_differences.elementary import floor -from devito.symbolics import retrieve_function_carriers, INT +from devito.symbolics import retrieve_function_carriers, retrieve_functions, INT from devito.tools import as_tuple, flatten from devito.types import (ConditionalDimension, Eq, Inc, Evaluable, Symbol, CustomDimension) @@ -14,6 +15,18 @@ __all__ = ['LinearInterpolator', 'PrecomputedInterpolator'] +def check_radius(func): + @wraps(func) + def wrapper(interp, *args, **kwargs): + r = interp.sfunction.r + funcs = set().union(*[retrieve_functions(a) for a in args]) + so = min({f.space_order for f in funcs if not f.is_SparseFunction} or {r}) + if so < r: + raise ValueError("Space order %d smaller than interpolation r %d" % (so, r)) + return func(interp, *args, **kwargs) + return wrapper + + class UnevaluatedSparseOperation(sympy.Expr, Evaluable): """ @@ -209,6 +222,7 @@ def _interp_idx(self, variables, implicit_dims=None): return idx_subs, temps + @check_radius def interpolate(self, expr, increment=False, self_subs={}, implicit_dims=None): """ Generate equations interpolating an arbitrary expression into ``self``. @@ -226,6 +240,7 @@ def interpolate(self, expr, increment=False, self_subs={}, implicit_dims=None): """ return Interpolation(expr, increment, implicit_dims, self_subs, self) + @check_radius def inject(self, field, expr, implicit_dims=None): """ Generate equations injecting an arbitrary expression into a field. diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 90816fd7ad..f268467dda 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -1008,6 +1008,8 @@ class PrecomputedSparseFunction(AbstractSparseFunction): uses `*args` to (re-)create the Dimension arguments of the symbolic object. """ + is_SparseFunction = True + _sub_functions = ('gridpoints', 'coordinates', 'interpolation_coeffs') __rkwargs__ = (AbstractSparseFunction.__rkwargs__ + @@ -1173,6 +1175,8 @@ class PrecomputedSparseTimeFunction(AbstractSparseTimeFunction, uses ``*args`` to (re-)create the Dimension arguments of the symbolic object. """ + is_SparseTimeFunction = True + __rkwargs__ = tuple(filter_ordered(AbstractSparseTimeFunction.__rkwargs__ + PrecomputedSparseFunction.__rkwargs__)) diff --git a/tests/test_dle.py b/tests/test_dle.py index 03c0b533c9..8e37947a36 100644 --- a/tests/test_dle.py +++ b/tests/test_dle.py @@ -709,7 +709,7 @@ def test_scheduling(self): """ grid = Grid(shape=(11, 11)) - u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=0) + u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=1) sf1 = SparseTimeFunction(name='s', grid=grid, npoint=1, nt=5) eqns = [Eq(u.forward, u + 1)] diff --git a/tests/test_gpu_common.py b/tests/test_gpu_common.py index 071005464d..3bf641a120 100644 --- a/tests/test_gpu_common.py +++ b/tests/test_gpu_common.py @@ -1424,7 +1424,7 @@ def test_empty_arrays(self): """ grid = Grid(shape=(4, 4), extent=(3.0, 3.0)) - f = TimeFunction(name='f', grid=grid, space_order=0) + f = TimeFunction(name='f', grid=grid, space_order=1) f.data[:] = 1. sf1 = SparseTimeFunction(name='sf1', grid=grid, npoint=0, nt=10) sf2 = SparseTimeFunction(name='sf2', grid=grid, npoint=0, nt=10, diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 97d86c1759..efde867402 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -14,19 +14,19 @@ import scipy.sparse -def unit_box(name='a', shape=(11, 11), grid=None): +def unit_box(name='a', shape=(11, 11), grid=None, space_order=1): """Create a field with value 0. to 1. in each dimension""" grid = grid or Grid(shape=shape) - a = Function(name=name, grid=grid) + a = Function(name=name, grid=grid, space_order=space_order) dims = tuple([np.linspace(0., 1., d) for d in shape]) a.data[:] = np.meshgrid(*dims)[1] return a -def unit_box_time(name='a', shape=(11, 11)): +def unit_box_time(name='a', shape=(11, 11), space_order=1): """Create a field with value 0. to 1. in each dimension""" grid = Grid(shape=shape) - a = TimeFunction(name=name, grid=grid, time_order=1) + a = TimeFunction(name=name, grid=grid, time_order=1, space_order=space_order) dims = tuple([np.linspace(0., 1., d) for d in shape]) a.data[0, :] = np.meshgrid(*dims)[1] a.data[1, :] = np.meshgrid(*dims)[1] @@ -117,16 +117,15 @@ def test_precomputed_interpolation(r): origin = (0, 0) grid = Grid(shape=shape, origin=origin) - r = 2 # Constant for linear interpolation - # because we interpolate across 2 neighbouring points in each dimension def init(data): + # This is data with halo so need to shift to match the m.data expectations for i in range(data.shape[0]): for j in range(data.shape[1]): - data[i, j] = sin(grid.spacing[0]*i) + sin(grid.spacing[1]*j) + data[i, j] = sin(grid.spacing[0]*(i-r)) + sin(grid.spacing[1]*(j-r)) return data - m = Function(name='m', grid=grid, initializer=init, space_order=0) + m = Function(name='m', grid=grid, initializer=init, space_order=r) gridpoints, interpolation_coeffs = precompute_linear_interpolation(points, grid, origin, @@ -154,10 +153,8 @@ def test_precomputed_interpolation_time(r): origin = (0, 0) grid = Grid(shape=shape, origin=origin) - r = 2 # Constant for linear interpolation - # because we interpolate across 2 neighbouring points in each dimension - u = TimeFunction(name='u', grid=grid, space_order=0, save=5) + u = TimeFunction(name='u', grid=grid, space_order=r, save=5) for it in range(5): u.data[it, :] = it @@ -190,11 +187,7 @@ def test_precomputed_injection(r): origin = (0, 0) result = 0.25 - # Constant for linear interpolation - # because we interpolate across 2 neighbouring points in each dimension - r = 2 - - m = unit_box(shape=shape) + m = unit_box(shape=shape, space_order=r) m.data[:] = 0. gridpoints, interpolation_coeffs = precompute_linear_interpolation(coords, @@ -228,11 +221,7 @@ def test_precomputed_injection_time(r): result = 0.25 nt = 20 - # Constant for linear interpolation - # because we interpolate across 2 neighbouring points in each dimension - r = 2 - - m = unit_box_time(shape=shape) + m = unit_box_time(shape=shape, space_order=r) m.data[:] = 0. gridpoints, interpolation_coeffs = precompute_linear_interpolation(coords, @@ -761,3 +750,16 @@ def test_inject_function(): for i in [0, 1, 3, 4]: for j in [0, 1, 3, 4]: assert u.data[1, i, j] == 0 + + +def test_interpolation_radius(): + nt = 11 + + grid = Grid(shape=(5, 5)) + u = TimeFunction(name="u", grid=grid, space_order=0) + src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1) + try: + src.interpolate(u) + assert False + except ValueError: + assert True diff --git a/tests/test_mpi.py b/tests/test_mpi.py index ab7092ba1a..46b8e18fa5 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -1501,7 +1501,7 @@ def test_injection_wodup(self): """ grid = Grid(shape=(4, 4), extent=(3.0, 3.0)) - f = Function(name='f', grid=grid, space_order=0) + f = Function(name='f', grid=grid, space_order=1) f.data[:] = 0. coords = np.array([(0.5, 0.5), (0.5, 2.5), (2.5, 0.5), (2.5, 2.5)]) sf = SparseFunction(name='sf', grid=grid, npoint=len(coords), coordinates=coords) @@ -1536,7 +1536,7 @@ def test_injection_wodup_wtime(self): grid = Grid(shape=(4, 4), extent=(3.0, 3.0)) save = 3 - f = TimeFunction(name='f', grid=grid, save=save, space_order=0) + f = TimeFunction(name='f', grid=grid, save=save, space_order=1) f.data[:] = 0. coords = np.array([(0.5, 0.5), (0.5, 2.5), (2.5, 0.5), (2.5, 2.5)]) sf = SparseTimeFunction(name='sf', grid=grid, nt=save, @@ -1611,7 +1611,7 @@ def test_injection_dup(self): def test_interpolation_wodup(self): grid = Grid(shape=(4, 4), extent=(3.0, 3.0)) - f = Function(name='f', grid=grid, space_order=0) + f = Function(name='f', grid=grid, space_order=1) f.data[:] = 4. coords = [(0.5, 0.5), (0.5, 2.5), (2.5, 0.5), (2.5, 2.5)] sf = SparseFunction(name='sf', grid=grid, npoint=len(coords), coordinates=coords) diff --git a/tests/test_operator.py b/tests/test_operator.py index b2188d007c..e685762a77 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -521,7 +521,7 @@ def test_sparsefunction_inject(self): Test injection of a SparseFunction into a Function """ grid = Grid(shape=(11, 11)) - u = Function(name='u', grid=grid, space_order=0) + u = Function(name='u', grid=grid, space_order=1) sf1 = SparseFunction(name='s', grid=grid, npoint=1) op = Operator(sf1.inject(u, expr=sf1)) @@ -542,7 +542,7 @@ def test_sparsefunction_interp(self): Test interpolation of a SparseFunction from a Function """ grid = Grid(shape=(11, 11)) - u = Function(name='u', grid=grid, space_order=0) + u = Function(name='u', grid=grid, space_order=1) sf1 = SparseFunction(name='s', grid=grid, npoint=1) op = Operator(sf1.interpolate(u)) @@ -563,7 +563,7 @@ def test_sparsetimefunction_interp(self): Test injection of a SparseTimeFunction into a TimeFunction """ grid = Grid(shape=(11, 11)) - u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=0) + u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=1) sf1 = SparseTimeFunction(name='s', grid=grid, npoint=1, nt=5) op = Operator(sf1.interpolate(u)) @@ -586,7 +586,7 @@ def test_sparsetimefunction_inject(self): Test injection of a SparseTimeFunction from a TimeFunction """ grid = Grid(shape=(11, 11)) - u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=0) + u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=1) sf1 = SparseTimeFunction(name='s', grid=grid, npoint=1, nt=5) op = Operator(sf1.inject(u, expr=3*sf1)) @@ -611,7 +611,7 @@ def test_sparsetimefunction_inject_dt(self): Test injection of the time deivative of a SparseTimeFunction into a TimeFunction """ grid = Grid(shape=(11, 11)) - u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=0) + u = TimeFunction(name='u', grid=grid, time_order=2, save=5, space_order=1) sf1 = SparseTimeFunction(name='s', grid=grid, npoint=1, nt=5, time_order=2)