From 77edeaf3ba05a7dd66c7c2d5622dc1df5b0d6b68 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Wed, 11 Mar 2026 12:40:35 +0000 Subject: [PATCH 1/6] compiler: Patch degenerating_dimensions to catch SteppingDimensions over size 1 buffer --- devito/ir/support/basic.py | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py index 3bad983426..731f3bc386 100644 --- a/devito/ir/support/basic.py +++ b/devito/ir/support/basic.py @@ -385,6 +385,9 @@ def distance(self, other): # it's provable that they never intersect if sai and sit == oit and disjoint_test(self[n], other[n], sai, sit): return Vector(S.ImaginaryUnit) + + # if 'p_bck_xx,[t, x + 9, y + 16, z + 16' in str(self): + # from IPython import embed; embed() # Compute the distance along the current IterationInterval if self.function._mem_shared: @@ -392,7 +395,9 @@ def distance(self, other): # objects falls back to zero, as any other value would be # nonsensical ret.append(S.Zero) - elif degenerating_dimensions(sai, oai): + # FIXME: This should have returned True + # elif degenerating_indices(sai, oai): + elif degenerating_indices(self[n], other[n], self.function): # Special case: `sai` and `oai` may be different symbolic objects # but they can be proved to systematically generate the same value ret.append(S.Zero) @@ -1425,17 +1430,38 @@ def disjoint_test(e0, e1, d, it): return not bool(i0.intersect(i1)) -def degenerating_dimensions(d0, d1): +def degenerating_indices(i0, i1, function): """ - True if `d0` and `d1` are Dimensions that are possibly symbolically + True if `i0` and `i1` are indices that are possibly symbolically different, but they can be proved to systematically degenerate to the same value, False otherwise. """ # Case 1: ModuloDimensions of size 1 try: - if d0.is_Modulo and d1.is_Modulo and d0.modulo == d1.modulo == 1: + if i0.is_Modulo and i1.is_Modulo and i0.modulo == i1.modulo == 1: return True except AttributeError: pass + # Case 2: SteppingDimension corresponding to buffer of size 1 + # if str(i1) == "t - 1" and str(i0) == 't': + # from IPython import embed; embed() + + # Extract dimension from both IndexAccessFunctions -> d0, d1 + try: + d0 = i0.d + except AttributeError: + d0 = i0 + + try: + d1 = i1.d + except AttributeError: + d1 = i1 + + with suppress(AttributeError): + if d0 is d1 and d0.is_Stepping: + # Buffer is size 1 + if function._size_domain[d0] == 1: + return True + return False From d503c557bf5385e53d0636447376089673844181 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Wed, 11 Mar 2026 14:56:11 +0000 Subject: [PATCH 2/6] tests: Add a test for the Buffer(1) issue --- tests/test_fission.py | 68 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 67 insertions(+), 1 deletion(-) diff --git a/tests/test_fission.py b/tests/test_fission.py index 3c6263e2e1..d27278d7a3 100644 --- a/tests/test_fission.py +++ b/tests/test_fission.py @@ -2,8 +2,11 @@ from conftest import assert_structure from devito import ( - Eq, Function, Grid, Inc, Operator, SubDimension, SubDomain, TimeFunction, solve + Buffer, Eq, Function, Grid, Inc, Operator, SubDimension, SubDomain, TimeFunction, + solve ) +from devito.ir.iet import retrieve_iteration_tree +from devito.ir.support.properties import PARALLEL def test_issue_1725(): @@ -126,3 +129,66 @@ def test_issue_1921(): op1.apply(time_m=1, time_M=5, g=g1) assert np.all(g.data == g1.data) + + +def test_buffer1_fissioning(): + """ + Tests an edge case whereby inability to spot the equivalence of + `f.forward`/`backward` and `f` when using `Buffer(1)` would cause + data dependance analysis to incorrectly determine that `x` and `y` + loops should be sequential. + """ + so = 2 + + class SD0(SubDomain): + name = 'sd0' + + def __init__(self, so, **kwargs): + self.so = so + super().__init__(**kwargs) + + def define(self, dimensions): + map_d = {d: d for d in dimensions} + map_d[dimensions[-1]] = ('left', self.so) + return map_d + + grid = Grid(shape=(101, 101, 101)) + z = grid.dimensions[-1] + + sd0 = SD0(so, grid=grid) + iz = sd0.dimensions[-1] + + f0 = TimeFunction(name='f0', grid=grid, + space_order=so, save=Buffer(1)) + f1 = TimeFunction(name='f1', grid=grid, + space_order=so, save=Buffer(1)) + + f2 = TimeFunction(name='f2', grid=grid, + space_order=so, save=Buffer(1)) + f3 = TimeFunction(name='f3', grid=grid, + space_order=so, save=Buffer(1)) + + # Free-surface like conditions plus derivative creates a dependence + eq0 = Eq(f0.backward.subs(z, iz-so), -f0.backward.subs(z, so-iz)) + eq1 = Eq(f1.backward.subs(z, iz-so), -f1.backward.subs(z, so-iz)) + + eq2 = Eq(f2, f0.dx + f2) + eq3 = Eq(f3, f1.dy + f3) + + eqs = [eq0, eq1, eq2, eq3] + + op = Operator(eqs) + + trees = retrieve_iteration_tree(op) + + # If the equivalence of f0.backward and f0 are not spotted by the + # compiler, data dependence analysis will determine x and y loops + # to be sequential. + for tree in trees: + for i in range(1, 3): + assert PARALLEL in tree[i].properties + + # If parallelisation failed, then there will also be no blocking on + # x and y + assert 'x0_blk0_size' in str(op.parameters) + assert 'y0_blk0_size' in str(op.parameters) From db5175bdcddace100b7cb8fc5190dc09545b07c6 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Wed, 11 Mar 2026 15:02:45 +0000 Subject: [PATCH 3/6] misc: Linting and tidy-up --- devito/ir/support/basic.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py index 731f3bc386..2f7ab430d5 100644 --- a/devito/ir/support/basic.py +++ b/devito/ir/support/basic.py @@ -385,9 +385,6 @@ def distance(self, other): # it's provable that they never intersect if sai and sit == oit and disjoint_test(self[n], other[n], sai, sit): return Vector(S.ImaginaryUnit) - - # if 'p_bck_xx,[t, x + 9, y + 16, z + 16' in str(self): - # from IPython import embed; embed() # Compute the distance along the current IterationInterval if self.function._mem_shared: @@ -1444,24 +1441,18 @@ def degenerating_indices(i0, i1, function): pass # Case 2: SteppingDimension corresponding to buffer of size 1 - # if str(i1) == "t - 1" and str(i0) == 't': - # from IPython import embed; embed() - # Extract dimension from both IndexAccessFunctions -> d0, d1 try: d0 = i0.d except AttributeError: d0 = i0 - try: d1 = i1.d except AttributeError: d1 = i1 with suppress(AttributeError): - if d0 is d1 and d0.is_Stepping: - # Buffer is size 1 - if function._size_domain[d0] == 1: - return True + if d0 is d1 and d0.is_Stepping and function._size_domain[d0] == 1: + return True return False From 2b47ebd3e8df29b58dd45377e0ea75cd00eb330d Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Wed, 11 Mar 2026 15:32:01 +0000 Subject: [PATCH 4/6] misc: Tidy-up --- devito/ir/support/basic.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py index 2f7ab430d5..8357b1b05f 100644 --- a/devito/ir/support/basic.py +++ b/devito/ir/support/basic.py @@ -392,8 +392,6 @@ def distance(self, other): # objects falls back to zero, as any other value would be # nonsensical ret.append(S.Zero) - # FIXME: This should have returned True - # elif degenerating_indices(sai, oai): elif degenerating_indices(self[n], other[n], self.function): # Special case: `sai` and `oai` may be different symbolic objects # but they can be proved to systematically generate the same value From 48e0c0a37283089f79be7979c71c554e3ae118ef Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Mon, 16 Mar 2026 12:45:19 +0000 Subject: [PATCH 5/6] tests: Fix up some halo tests with Buffer --- tests/test_mpi.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/tests/test_mpi.py b/tests/test_mpi.py index ba3636c682..24ac679812 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -1887,10 +1887,12 @@ def test_enforce_haloupdate_if_unwritten_function(self, mode): @pytest.mark.parallel(mode=1) def test_haloupdate_buffer1(self, mode): grid = Grid(shape=(4, 4)) - x, y = grid.dimensions - u = TimeFunction(name='u', grid=grid, time_order=1, save=Buffer(1)) - v = TimeFunction(name='v', grid=grid, time_order=1, save=Buffer(1)) + # With order 2 (1) forward derivatives, the loops can (and will) be fused, + # removing the need for a halo update so space_order=4 is used to ensure + # derivatives are centred resulting in parallel loops and halo updates + u = TimeFunction(name='u', grid=grid, time_order=1, space_order=4, save=Buffer(1)) + v = TimeFunction(name='v', grid=grid, time_order=1, space_order=4, save=Buffer(1)) eqns = [Eq(u.forward, div(v) + 1.), Eq(v.forward, div(u.forward) + 1.)] @@ -1909,22 +1911,22 @@ def test_haloupdate_buffer1(self, mode): @pytest.mark.parallel(mode=1) @pytest.mark.parametrize('sz,fwd,expr,exp0,exp1,args', [ (1, True, 'rec.interpolate(v2)', 3, 2, ('v1', 'v2')), - (1, True, 'Eq(v3.forward, v2.laplace + 1)', 1, 1, ('v2',)), - (1, True, 'Eq(v3.forward, v2.forward.laplace + 1)', 3, 2, ('v1', 'v2',)), - (2, True, 'Eq(v3.forward, v2.forward.laplace + 1)', 3, 2, ('v1', 'v2',)), + (1, True, 'Eq(v3.forward, v2.laplace + 1)', 3, 2, ('v1', 'v2')), + (1, True, 'Eq(v3.forward, v2.forward.laplace + 1)', 3, 2, ('v1', 'v2')), + (2, True, 'Eq(v3.forward, v2.forward.laplace + 1)', 3, 2, ('v1', 'v2')), (1, False, 'rec.interpolate(v2)', 3, 2, ('v1', 'v2')), - (1, False, 'Eq(v3.backward, v2.laplace + 1)', 1, 1, ('v2',)), - (1, False, 'Eq(v3.backward, v2.backward.laplace + 1)', 3, 2, ('v1', 'v2',)), - (2, False, 'Eq(v3.backward, v2.backward.laplace + 1)', 3, 2, ('v1', 'v2',)), + (1, False, 'Eq(v3.backward, v2.laplace + 1)', 3, 2, ('v1', 'v2')), + (1, False, 'Eq(v3.backward, v2.backward.laplace + 1)', 3, 3, ('v2', 'v1', 'v2')), + (2, False, 'Eq(v3.backward, v2.backward.laplace + 1)', 3, 3, ('v2', 'v1', 'v2')), ]) def test_haloupdate_buffer_cases(self, sz, fwd, expr, exp0, exp1, args, mode): grid = Grid((65, 65, 65), topology=('*', 1, '*')) - v1 = TimeFunction(name='v1', grid=grid, space_order=2, time_order=1, + v1 = TimeFunction(name='v1', grid=grid, space_order=4, time_order=1, save=Buffer(1)) - v2 = TimeFunction(name='v2', grid=grid, space_order=2, time_order=1, + v2 = TimeFunction(name='v2', grid=grid, space_order=4, time_order=1, save=Buffer(1)) - v3 = TimeFunction(name='v3', grid=grid, space_order=2, time_order=1, # noqa + v3 = TimeFunction(name='v3', grid=grid, space_order=4, time_order=1, # noqa save=Buffer(1)) rec = SparseTimeFunction(name='rec', grid=grid, nt=500, npoint=65) # noqa From 4da6f326780c143370d18af8f1aa8f1a9b787d2e Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Mon, 16 Mar 2026 14:37:16 +0000 Subject: [PATCH 6/6] tests: Introduce assert_structure to fission test --- tests/test_fission.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_fission.py b/tests/test_fission.py index d27278d7a3..fd54890510 100644 --- a/tests/test_fission.py +++ b/tests/test_fission.py @@ -192,3 +192,7 @@ def define(self, dimensions): # x and y assert 'x0_blk0_size' in str(op.parameters) assert 'y0_blk0_size' in str(op.parameters) + + # Two loop nests: free-surface-like and update-like + assert_structure(op, ['t,x,y,z', 't,x0_blk0,y0_blk0,x,y,z'], + 't,x,y,z,x0_blk0,y0_blk0,x,y,z')