diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index 9c56a1aa7d..1e395e1a38 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -1111,6 +1111,58 @@ def dimension(self): weights = Array.initvalue + @cached_property + def symmetry(self): + """ + The sign relating mirrored coefficients: `1` for symmetric weights, + `-1` for antisymmetric weights, or `None` otherwise. + + An odd antisymmetric sequence must have a zero center coefficient. + """ + weights = self.weights + if not weights: + return None + + middle = len(weights) // 2 + pairs = tuple(zip(weights[:middle], + reversed(weights[middle + len(weights) % 2:]), + strict=True)) + + if all(left.equals(right) is True for left, right in pairs): + return sympy.S.One + if all(left.equals(-right) is True for left, right in pairs) and \ + (len(weights) % 2 == 0 or weights[middle].equals(0) is True): + return sympy.S.NegativeOne + + return None + + @cached_property + def coefficient_span(self): + """ + Paired weights and their offset within the original stencil. + + Symmetric or antisymmetric weights retain their full span, including + symmetric zero padding. Otherwise, trim known zeros at both ends and + return rebuilt weights together with their starting offset. The + resulting weights may still have no symmetry. + """ + if self.symmetry is not None: + return self, 0 + + weights = self.weights + first, last = 0, len(weights) + while first < last and weights[first].is_zero is True: + first += 1 + while last > first and weights[last - 1].is_zero is True: + last -= 1 + + if first == 0 and last == len(weights): + return self, 0 + + d = self.dimension + dimension = d._rebuild(_min=d._min + first, _max=d._min + last - 1) + return self._rebuild(dimensions=dimension, initvalue=weights[first:last]), first + def _xreplace(self, rule): if self in rule: return rule[self], True diff --git a/devito/passes/clusters/blocking.py b/devito/passes/clusters/blocking.py index fea502fbfb..0e6a491abd 100644 --- a/devito/passes/clusters/blocking.py +++ b/devito/passes/clusters/blocking.py @@ -541,7 +541,6 @@ def schedule(self, dims, clusters): else: umt = self.umt - unbound = self.unbound and umt.is_multi umt.iter() diff --git a/devito/symbolics/unevaluation.py b/devito/symbolics/unevaluation.py index cc60854e7f..9fc6e702f3 100644 --- a/devito/symbolics/unevaluation.py +++ b/devito/symbolics/unevaluation.py @@ -1,6 +1,6 @@ import sympy -__all__ = ['Add', 'Mul', 'Pow'] +__all__ = ['Add', 'Mod', 'Mul', 'Pow'] class UnevaluableMixin: @@ -13,6 +13,10 @@ class Add(sympy.Add, UnevaluableMixin): __new__ = UnevaluableMixin.__new__ +class Mod(sympy.Mod, UnevaluableMixin): + __new__ = UnevaluableMixin.__new__ + + class Mul(sympy.Mul, UnevaluableMixin): __new__ = UnevaluableMixin.__new__ diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index fa1645ddff..b5169efbcd 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -941,6 +941,54 @@ def test_unexpand_space_interp_w_saved_timefunc(self): assert_structure(op, ['t,x,y,z', 't,x,y,z,i1', 't,x,y,z,i1,i0']) +class TestWeights: + + @pytest.mark.parametrize('coeffs,expected', [ + ((1, 2, -2, -1), -1), + ((1, 2, 0, -2, -1), -1), + ((1, 2, 2, 1), 1), + ((1, 2, 3, 2, 1), 1), + ((1, 2, 3, -2, -1), None), + ((0, 0, 0), 1), + ((0,), 1), + ((3,), 1), + ((Symbol('a'), 0, -Symbol('a')), -1), + ((Symbol('a'), Symbol('b'), Symbol('a')), 1), + ((Symbol('a'), Symbol('b')), None), + ]) + def test_symmetry(self, coeffs, expected): + i = StencilDimension('i', 0, len(coeffs) - 1) + weights = Weights(name='w', dimensions=i, initvalue=coeffs) + + assert weights.symmetry == expected + assert weights._rebuild().symmetry == expected + + @pytest.mark.parametrize('coeffs,expected,offset,sign', [ + ((0, 1, 2, -2, -1), (1, 2, -2, -1), 1, -1), + ((1, 2, -2, -1, 0), (1, 2, -2, -1), 0, -1), + ((0, 0, 1, 2, 2, 1, 0), (1, 2, 2, 1), 2, 1), + ((0, 1, 2, -2, -1, 0), (0, 1, 2, -2, -1, 0), 0, -1), + ((0, 1, 2, 3, 0), (1, 2, 3), 1, None), + ((1, 2, 3), (1, 2, 3), 0, None), + ((0, 0, 0), (0, 0, 0), 0, 1), + ((0, Symbol('a'), -Symbol('a')), (Symbol('a'), -Symbol('a')), 1, -1), + ]) + def test_coefficient_span(self, coeffs, expected, offset, sign): + i = StencilDimension('i', -3, len(coeffs) - 4) + weights = Weights(name='w', dimensions=i, initvalue=coeffs) + span, start = weights.coefficient_span + + assert span.weights == expected + assert span.symmetry == sign + assert start == offset + assert span.dimension._min == i._min + offset + assert span.dimension._size == len(expected) + assert weights.coefficient_span[0] is span + assert weights._rebuild().coefficient_span == (span, start) + if coeffs == expected: + assert span is weights + + class TestTwoStageEvaluation: def test_exceptions(self): diff --git a/tests/test_gpu_openacc.py b/tests/test_gpu_openacc.py index 2dba097684..eddae9fb1b 100644 --- a/tests/test_gpu_openacc.py +++ b/tests/test_gpu_openacc.py @@ -144,11 +144,10 @@ def test_multiple_tile_sizes(self, par_tile): assert trees[3][1].pragmas[0].ccode.value ==\ f'acc parallel loop {sclause} present(src,src_gp,src_wx,src_wy,src_wz,u)' - def test_short_multi_tile_keeps_outer_dim_blocked(self): + def test_short_multi_tile_keeps_outer_dim_on_device(self): """ - A multi `par-tile` entry shorter than the nest it lands on must not cost - the outermost Dimension its BlockDimension: on a device, dropping it - would leave `x` iterated outside the offloaded nest. + A short multi `par-tile` leaves the outermost Dimension unblocked. + Its full loop must stay inside the OpenACC region. """ grid = Grid(shape=(8, 8, 8)) @@ -166,9 +165,21 @@ def test_short_multi_tile_keeps_outer_dim_blocked(self): 'advanced', {'par-tile': par_tile, 'blocklevels': 1, 'blockinner': True})) - bns, _ = assert_blocking(op, {'x0_blk0', 'x1_blk0'}) + bns, _ = assert_blocking(op, {'x0_blk0', 'y1_blk0'}) - expected = ((4, 4, 32), (4, 4, 16)) + trees = retrieve_iteration_tree(op) + assert len(trees) == 2 + assert all(tree[0].dim is grid.time_dim for tree in trees) + assert trees[0][1].pragmas[0].ccode.value ==\ + 'acc parallel loop tile(32,4,4) present(u)' + + x = grid.dimensions[0] + assert trees[1][1].dim is x + assert trees[1][1].limits == (x.symbolic_min, x.symbolic_max, 1) + assert trees[1][1].pragmas[0].ccode.value ==\ + 'acc parallel loop tile(16,4,4) present(u,v)' + + expected = ((4, 4, 32), (4, 16)) for root, v in zip(bns.values(), expected, strict=True): iters = FindNodes(Iteration).visit(root) iters = [i for i in iters if i.dim.is_Block and i.dim._depth == 1] diff --git a/tests/test_pickle.py b/tests/test_pickle.py index d29dbb15b5..aea96df5e9 100644 --- a/tests/test_pickle.py +++ b/tests/test_pickle.py @@ -22,6 +22,7 @@ CallFromPointer, Cast, DefFunction, FieldFromPointer, IntDiv, ListInitializer, SizeOf, indexify, pow_to_mul ) +from devito.symbolics.unevaluation import Mod from devito.tools import EnrichedTuple from devito.types import ( Array, ComponentAccess, CustomDimension, DefaultDimension, DeviceID, FIndexed, @@ -720,6 +721,15 @@ def test_local_sum(self, pickle): tuple(map(str, reduction.conditionals.values())) assert rebuilt.free_symbols.isdisjoint(rebuilt.bound_symbols) + @pytest.mark.parametrize('args', [(0, 3), (1, 3), (3, 3), (Symbol('i'), 3)]) + def test_unevaluated_mod(self, pickle, args): + expr = Mod(*args) + + rebuilt = pickle.loads(pickle.dumps(expr)) + + assert rebuilt == expr + assert rebuilt.args == args + class TestAdvanced: