From 96f0d51e0c0fde22fe161068388940628083e718 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Wed, 26 Aug 2026 16:25:13 +0000 Subject: [PATCH 1/7] compiler: Tweak preprocess to avoid erroneous matching of subdomains when inserting halo exchanges --- devito/ir/stree/algorithms.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/devito/ir/stree/algorithms.py b/devito/ir/stree/algorithms.py index fb313f7640..f9828e17e8 100644 --- a/devito/ir/stree/algorithms.py +++ b/devito/ir/stree/algorithms.py @@ -206,6 +206,26 @@ def preprocess(clusters, options=None, **kwargs): diff = dims - distributed_aindices intersection = dims & distributed_aindices + # FIXME: I think this might be borked with MPI? Might erroneously skip + # a halo exchange which is needed for a cluster + # Doesn't throw any errors, but I need to logically check this + + # The HaloScheme from a previous cluster cannot be reused in the case + # that its `distributed_aindices` contain a SubDimension which is not + # found in the current Cluster being inspected but which shares the + # same root Dimension as a SubDimension found in the current cluster. + # As such, check that SubDimensions in `distributed_aindices` match + # those in `dims` where they share a root Dimension. If not, then + # skip this entry in queue. + d_by_root = {d.root: d for d in dims if d.is_Sub} + conflict = any( + e.is_Sub and d_by_root.get(e.root, e) is not e + for e in distributed_aindices + ) + + if conflict: + continue + if all(c1.guards.get(d) == c.guards.get(d) for d in diff) and \ len(intersection) > 0: found.append(c1) From 039d1f7756d0b8816b022b7ca59b2e4493deeded Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Fri, 28 Aug 2026 09:33:28 +0000 Subject: [PATCH 2/7] misc: Minor comment tweaks --- devito/ir/stree/algorithms.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/devito/ir/stree/algorithms.py b/devito/ir/stree/algorithms.py index f9828e17e8..a3c933ffec 100644 --- a/devito/ir/stree/algorithms.py +++ b/devito/ir/stree/algorithms.py @@ -206,9 +206,10 @@ def preprocess(clusters, options=None, **kwargs): diff = dims - distributed_aindices intersection = dims & distributed_aindices - # FIXME: I think this might be borked with MPI? Might erroneously skip + # TODO: Double-check this isn't borked with MPI? Might erroneously skip # a halo exchange which is needed for a cluster - # Doesn't throw any errors, but I need to logically check this + # Doesn't throw any errors with the basic test, but I need to logically + # check this # The HaloScheme from a previous cluster cannot be reused in the case # that its `distributed_aindices` contain a SubDimension which is not From 9cdfb2dfec34792e190321d176c0d01ac733ab20 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Fri, 28 Aug 2026 10:44:39 +0000 Subject: [PATCH 3/7] compiler: Tweak to include all dimensions, not just subdimensions --- devito/ir/stree/algorithms.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/devito/ir/stree/algorithms.py b/devito/ir/stree/algorithms.py index a3c933ffec..a25fc31812 100644 --- a/devito/ir/stree/algorithms.py +++ b/devito/ir/stree/algorithms.py @@ -212,15 +212,15 @@ def preprocess(clusters, options=None, **kwargs): # check this # The HaloScheme from a previous cluster cannot be reused in the case - # that its `distributed_aindices` contain a SubDimension which is not + # that its `distributed_aindices` contain a Dimension which is not # found in the current Cluster being inspected but which shares the - # same root Dimension as a SubDimension found in the current cluster. - # As such, check that SubDimensions in `distributed_aindices` match + # same `root` as a Dimension found in the current cluster. + # As such, check that Dimensions in `distributed_aindices` match # those in `dims` where they share a root Dimension. If not, then # skip this entry in queue. - d_by_root = {d.root: d for d in dims if d.is_Sub} + d_by_root = {d.root: d for d in dims} conflict = any( - e.is_Sub and d_by_root.get(e.root, e) is not e + d_by_root.get(e.root, e) is not e for e in distributed_aindices ) From 4c87825fe7de1c8f690d247a6873d0b44d4f4966 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Fri, 28 Aug 2026 16:38:45 +0000 Subject: [PATCH 4/7] compiler: Fix up regression in interpolation tests with functions on subdomains --- devito/ir/stree/algorithms.py | 30 ++++++++++++++++-------------- tests/test_mpi.py | 6 ++++++ 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/devito/ir/stree/algorithms.py b/devito/ir/stree/algorithms.py index a25fc31812..febebf5f20 100644 --- a/devito/ir/stree/algorithms.py +++ b/devito/ir/stree/algorithms.py @@ -1,3 +1,4 @@ +from collections import defaultdict from itertools import groupby from anytree import findall @@ -206,21 +207,22 @@ def preprocess(clusters, options=None, **kwargs): diff = dims - distributed_aindices intersection = dims & distributed_aindices - # TODO: Double-check this isn't borked with MPI? Might erroneously skip - # a halo exchange which is needed for a cluster - # Doesn't throw any errors with the basic test, but I need to logically - # check this - - # The HaloScheme from a previous cluster cannot be reused in the case - # that its `distributed_aindices` contain a Dimension which is not - # found in the current Cluster being inspected but which shares the - # same `root` as a Dimension found in the current cluster. - # As such, check that Dimensions in `distributed_aindices` match - # those in `dims` where they share a root Dimension. If not, then - # skip this entry in queue. - d_by_root = {d.root: d for d in dims} + # A non-empty `intersection` doesn't guarantee `c1` belongs to + # `c`: two SubDomains sharing an axis's (side, thickness) alias + # to the same cached SubDimension, so the overlap may be one + # incidental axis while the rest belongs to an unrelated + # SubDomain. Require instead that, for every root shared + # between `dims` and `distributed_aindices`, the latter's + # Dimension is exactly one of `dims`'s own for that root (a + # root missing from `dims` is still fine, e.g. `t, f` + # triggering a halo for `t, x, y, z, f`). `dims` may hold + # several Dimensions per root (e.g. a point Dimension plus a + # derived radius CustomDimension), hence a set per root below. + d_by_root = defaultdict(set) + for d in dims: + d_by_root[d.root].add(d) conflict = any( - d_by_root.get(e.root, e) is not e + e.root in d_by_root and e not in d_by_root[e.root] for e in distributed_aindices ) diff --git a/tests/test_mpi.py b/tests/test_mpi.py index d5ad76323a..0f57859b80 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -1670,6 +1670,12 @@ def test_process_but_avoid_haloupdate_along_replicated(self, mode): assert len(calls) == 1 assert calls[0].arguments[0] is u + # TODO: Insert a test based off cpml-like reproducer 06 (including the alternative specs) + # TODO: This test should both check that the stree/iet doesn't have any halos inserted in inner iterations + # TODO: There is possibly already a function in this file or conftest which does this + # TODO: It should also check the operator compiles and runs (although it does not need to check resulting + # values as they are not physically meaningful). + @pytest.mark.parallel(mode=1) def test_conditional_dimension(self, mode): """ From 30ee739eb8259e7a40970d87acccf66e48e48baf Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Tue, 1 Sep 2026 09:33:16 +0000 Subject: [PATCH 5/7] tests: WIP tests and refactoring --- devito/ir/stree/algorithms.py | 39 ++++++----- tests/test_ir.py | 62 ++++++++++++++++- tests/test_mpi.py | 125 ++++++++++++++++++++++++++++++++-- 3 files changed, 202 insertions(+), 24 deletions(-) diff --git a/devito/ir/stree/algorithms.py b/devito/ir/stree/algorithms.py index febebf5f20..b8591db853 100644 --- a/devito/ir/stree/algorithms.py +++ b/devito/ir/stree/algorithms.py @@ -207,26 +207,14 @@ def preprocess(clusters, options=None, **kwargs): diff = dims - distributed_aindices intersection = dims & distributed_aindices + # TODO: Can check intersection length here and short-circuit early + # A non-empty `intersection` doesn't guarantee `c1` belongs to # `c`: two SubDomains sharing an axis's (side, thickness) alias # to the same cached SubDimension, so the overlap may be one # incidental axis while the rest belongs to an unrelated - # SubDomain. Require instead that, for every root shared - # between `dims` and `distributed_aindices`, the latter's - # Dimension is exactly one of `dims`'s own for that root (a - # root missing from `dims` is still fine, e.g. `t, f` - # triggering a halo for `t, x, y, z, f`). `dims` may hold - # several Dimensions per root (e.g. a point Dimension plus a - # derived radius CustomDimension), hence a set per root below. - d_by_root = defaultdict(set) - for d in dims: - d_by_root[d.root].add(d) - conflict = any( - e.root in d_by_root and e not in d_by_root[e.root] - for e in distributed_aindices - ) - - if conflict: + # SubDomain -- see `is_halo_scheme_conflicting`. + if is_halo_scheme_conflicting(dims, distributed_aindices): continue if all(c1.guards.get(d) == c.guards.get(d) for d in diff) and \ @@ -266,6 +254,25 @@ def preprocess(clusters, options=None, **kwargs): return processed +def is_halo_scheme_conflicting(dims, distributed_aindices): + """ + True if `distributed_aindices` cannot safely be attached to a Cluster + whose block-promoted itintervals are `dims`, because some Dimension in + `distributed_aindices` shares a root with a Dimension in `dims` without + being identical to it (a root missing from `dims` is still fine, e.g. + `t, f` triggering a halo for `t, x, y, z, f`). `dims` may hold several + Dimensions per root at once (e.g. a point Dimension plus a derived + radius CustomDimension), hence a set per root below. + """ + d_by_root = defaultdict(set) + for d in dims: + d_by_root[d.root].add(d) + return any( + e.root in d_by_root and e not in d_by_root[e.root] + for e in distributed_aindices + ) + + def reuse_partial_subtree(c0, c1, d=None): return c0.guards.get(d) == c1.guards.get(d) diff --git a/tests/test_ir.py b/tests/test_ir.py index 9f26ecb8c5..bc667b4a91 100644 --- a/tests/test_ir.py +++ b/tests/test_ir.py @@ -4,8 +4,8 @@ from conftest import EVAL, skipif # noqa from devito import ( # noqa - Constant, Dimension, Eq, Function, Grid, Inc, Operator, SubDimension, TimeFunction, - switchconfig + Constant, CustomDimension, Dimension, Eq, Function, Grid, Inc, Operator, + SubDimension, TimeFunction, switchconfig ) from devito.ir.cgen import ccode from devito.ir.clusters import Cluster, ClusterGroup @@ -13,6 +13,7 @@ from devito.ir.equations.algorithms import dimension_sort from devito.ir.iet import FindNodes, Iteration from devito.ir.stree import stree_build +from devito.ir.stree.algorithms import is_halo_scheme_conflicting from devito.ir.support.basic import ( AFFINE, IRREGULAR, REGULAR, IterationInstance, Scope, TimedAccess, Vector, mocksym0, mocksym1 @@ -1188,6 +1189,63 @@ def test_from_clusters_mixed_dtypes(self): assert len([i for i in stree.visit() if i.is_Iteration]) == 1 +class TestSubdomainHaloMatching: + + """ + Tests for `devito.ir.stree.algorithms.is_halo_scheme_conflicting`, the + logic that decides whether a "wild" HaloTouch Cluster's HaloScheme may be + attached to a given computational Cluster during `preprocess()`. + """ + + def test_cross_subdomain_aliasing_conflict(self): + """ + Two different SubDomains sharing one axis's (side, thickness) alias + to the same cached SubDimension on that axis; a mismatch on another + shared-root axis must still be flagged as conflicting. + """ + x, y = Dimension('x'), Dimension('y') + ix_shared = SubDimension.left('ix', x, 3) + iy_a = SubDimension.right('iy', y, 3) # domain A's own y + iy_b = SubDimension.left('iy', y, 3) # domain B's y -- different! + + dims = {ix_shared, iy_a} + distributed_aindices = {ix_shared, iy_b} + + assert is_halo_scheme_conflicting(dims, distributed_aindices) + + def test_bare_dimension_vs_subdimension_conflict(self): + """ + A bare (unrestricted) Dimension in `dims` colliding with an + unrelated SubDomain's SubDimension for the same root. + """ + x = Dimension('x') + xi = SubDimension.left('xi', x, 3) + + assert is_halo_scheme_conflicting({x}, {xi}) + + def test_no_conflict_when_root_missing_from_dims(self): + """ + A root present in `distributed_aindices` but absent from `dims` + entirely (e.g. an outer `t, f` Cluster triggering a halo for an + inner `t, x, y, z, f` Cluster) is not a conflict. + """ + x, f = Dimension('x'), Dimension('f') + + assert not is_halo_scheme_conflicting({x}, {x, f}) + + def test_multiple_dimensions_per_root_not_collapsed(self): + """ + `dims` may legitimately hold several Dimensions sharing one root at + once (e.g. a sparse-function point Dimension plus a CustomDimension + derived from it for an interpolation radius); this must not be + collapsed to a single arbitrary match. + """ + p = Dimension('p') + r = CustomDimension(name='r', symbolic_min=0, symbolic_max=1, parent=p) + + assert not is_halo_scheme_conflicting({p, r}, {p}) + + class TestClusterGroup: def test_eq_hash_include_ispace(self): diff --git a/tests/test_mpi.py b/tests/test_mpi.py index 0f57859b80..bc46892549 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -14,7 +14,8 @@ from devito.arch.compiler import OneapiCompiler from devito.data import LEFT, RIGHT from devito.ir.iet import ( - Call, Conditional, FindNodes, FindSymbols, Iteration, retrieve_iteration_tree + Call, Conditional, FindNodes, FindSymbols, HaloSpot, Iteration, + retrieve_iteration_tree ) from devito.ir.support.space import Backward, Forward from devito.mpi import MPI @@ -1107,6 +1108,115 @@ def check_halo_exchanges(op, exp0, exp1): return calls, tloop +def check_cpml_no_misplaced_halo(specs): + """ + Build an Operator modelled on a CPML-style absorbing-boundary formulation: + several overlapping-thickness SubDomains, each carrying a pair of auxiliary + VectorTimeFunctions updated via a "diagonal" derivative pattern (component + `i` derived along dimension `i`). This is the minimal pattern that exposed + a bug in `devito.ir.stree.algorithms.preprocess`, whereby the HaloScheme of + one SubDomain could be erroneously attached to a *different* SubDomain's + Cluster whenever the two happened to share a single axis's SubDimension + (SubDimensions are cached by `(name, parent, thickness, local)`, so two + unrelated SubDomains restricting one axis identically alias to the same + object). + + The Operator is built and run to completion; the values are not physically + meaningful and are not checked. The regression check is purely structural: + no halo-exchange node may end up nested *inside* a blocking (`is_Incr`) + Iteration. + """ + class SpeccedDomain(SubDomain): + def __init__(self, name, thickness, spec, **kwargs): + self.name = name + self.thickness = thickness + self.spec = spec + super().__init__(**kwargs) + + def define(self, dimensions): + retval = {} + for d, s in zip(dimensions, self.spec, strict=True): + if s == 'none': + # Unrestricted -- bare Dimension, not a SubDimension + retval[d] = d + elif s == 'middle': + retval[d] = (s, self.thickness, self.thickness) + else: + retval[d] = (s, self.thickness) + return retval + + so, to, nb = 4, 2, 3 + + grid = Grid(shape=(21, 21, 21), extent=(1., 1., 1.)) + domains = {spec: SpeccedDomain(f"x{spec[0]}_y{spec[1]}_z{spec[2]}", nb, spec, + grid=grid) + for spec in specs} + + p = TimeFunction(name='p', grid=grid, space_order=so, time_order=to) + + psi_eqs, zeta_eqs, p_eqs = [], [], [] + for v in domains.values(): + # Fields live on the base grid, not `v` -- `subdomain=v` is passed + # explicitly on each Eq instead. + psi = VectorTimeFunction(name=f"psi_{v.name}", grid=grid, space_order=so, + time_order=to, staggered=(None, None, None)) + zeta = VectorTimeFunction(name=f"zeta_{v.name}", grid=grid, space_order=so, + time_order=to, staggered=(None, None, None)) + psi_eqs.append(Eq(psi, 1, subdomain=v)) + + # "Diagonal" derivative pattern -- component `i` of zeta is the + # derivative of component `i` of psi along dimension `i`. This + # specific pattern is required to trigger the bug; grad()/div()/a + # plain vector add alone do not. + zeta_diag = VectorTimeFunction([getattr(psi[i], f"d{d.name}") + for i, d in enumerate(grid.dimensions)]) + zeta_eqs.append(Eq(zeta, zeta_diag, subdomain=v)) + p_eqs.append(Eq(p.forward, psi.div(), subdomain=v)) + + op = Operator(psi_eqs + zeta_eqs + p_eqs) + op.apply(time_M=1) + + # No misplaced halo exchanges: every halo-exchange node must sit above + # (never inside) any blocking Iteration. + incr_iterations = [i for i in FindNodes(Iteration).visit(op) if i.dim.is_Incr] + assert incr_iterations, "no blocking Iterations found -- check would be vacuous" + + if configuration['mpi']: + # HaloSpots have already been lowered into concrete Calls by mpiize() + halo_types = (HaloUpdateCall, HaloUpdateList) + else: + # HaloSpots remain in the IET as transparent (no-op) wrappers + halo_types = (HaloSpot,) + + for i in incr_iterations: + assert len(FindNodes(halo_types).visit(i)) == 0 + + +CPML_SPECS = [ + pytest.param( + [('left', 'right', 'left'), + ('left', 'right', 'middle'), + ('left', 'right', 'right'), + ('middle', 'left', 'right')], + id='full-restriction' + ), + pytest.param( + [('left', 'right', 'left'), + ('none', 'none', 'middle'), + ('none', 'none', 'right'), + ('middle', 'left', 'right')], + id='partial-none-xy' + ), + pytest.param( + [('left', 'right', 'left'), + ('left', 'right', 'none'), + ('left', 'right', 'right'), + ('middle', 'left', 'right')], + id='partial-none-z' + ), +] + + class TestCodeGeneration: @pytest.mark.parallel(mode=1) @@ -1670,11 +1780,14 @@ def test_process_but_avoid_haloupdate_along_replicated(self, mode): assert len(calls) == 1 assert calls[0].arguments[0] is u - # TODO: Insert a test based off cpml-like reproducer 06 (including the alternative specs) - # TODO: This test should both check that the stree/iet doesn't have any halos inserted in inner iterations - # TODO: There is possibly already a function in this file or conftest which does this - # TODO: It should also check the operator compiles and runs (although it does not need to check resulting - # values as they are not physically meaningful). + @pytest.mark.parametrize('specs', CPML_SPECS) + def test_cpml_no_misplaced_halo(self, specs): + check_cpml_no_misplaced_halo(specs) + + @pytest.mark.parametrize('specs', CPML_SPECS) + @pytest.mark.parallel(mode=1) + def test_cpml_no_misplaced_halo_mpi(self, specs, mode): + check_cpml_no_misplaced_halo(specs) @pytest.mark.parallel(mode=1) def test_conditional_dimension(self, mode): From d88604f0fc0d4f95ad12d1364e97fe8271174e07 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Tue, 1 Sep 2026 15:35:47 +0100 Subject: [PATCH 6/7] tests: Force blocking on MPI test to ensure it also works on CPU --- tests/test_mpi.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/tests/test_mpi.py b/tests/test_mpi.py index bc46892549..98e111ded1 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -1149,7 +1149,7 @@ def define(self, dimensions): grid = Grid(shape=(21, 21, 21), extent=(1., 1., 1.)) domains = {spec: SpeccedDomain(f"x{spec[0]}_y{spec[1]}_z{spec[2]}", nb, spec, - grid=grid) + grid=grid) for spec in specs} p = TimeFunction(name='p', grid=grid, space_order=so, time_order=to) @@ -1159,9 +1159,9 @@ def define(self, dimensions): # Fields live on the base grid, not `v` -- `subdomain=v` is passed # explicitly on each Eq instead. psi = VectorTimeFunction(name=f"psi_{v.name}", grid=grid, space_order=so, - time_order=to, staggered=(None, None, None)) + time_order=to, staggered=(None, None, None)) zeta = VectorTimeFunction(name=f"zeta_{v.name}", grid=grid, space_order=so, - time_order=to, staggered=(None, None, None)) + time_order=to, staggered=(None, None, None)) psi_eqs.append(Eq(psi, 1, subdomain=v)) # "Diagonal" derivative pattern -- component `i` of zeta is the @@ -1169,24 +1169,23 @@ def define(self, dimensions): # specific pattern is required to trigger the bug; grad()/div()/a # plain vector add alone do not. zeta_diag = VectorTimeFunction([getattr(psi[i], f"d{d.name}") - for i, d in enumerate(grid.dimensions)]) + for i, d in enumerate(grid.dimensions)]) zeta_eqs.append(Eq(zeta, zeta_diag, subdomain=v)) p_eqs.append(Eq(p.forward, psi.div(), subdomain=v)) - op = Operator(psi_eqs + zeta_eqs + p_eqs) + op = Operator(psi_eqs + zeta_eqs + p_eqs, + opt=('advanced', {'blockrelax': 'device-aware'})) op.apply(time_M=1) # No misplaced halo exchanges: every halo-exchange node must sit above # (never inside) any blocking Iteration. incr_iterations = [i for i in FindNodes(Iteration).visit(op) if i.dim.is_Incr] - assert incr_iterations, "no blocking Iterations found -- check would be vacuous" - - if configuration['mpi']: - # HaloSpots have already been lowered into concrete Calls by mpiize() - halo_types = (HaloUpdateCall, HaloUpdateList) - else: - # HaloSpots remain in the IET as transparent (no-op) wrappers - halo_types = (HaloSpot,) + assert incr_iterations, "no blocking Iterations found" + + # HaloSpots have already been lowered into concrete Calls by mpiize() in the MPI + # case, but HaloSpots remain in the IET as transparent (no-op) wrappers in the + # serial case. + halo_types = (HaloUpdateCall, HaloUpdateList) if configuration['mpi'] else (HaloSpot,) for i in incr_iterations: assert len(FindNodes(halo_types).visit(i)) == 0 From 0935d2c8f893c322a4386b1aa84a2c769492e08e Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Tue, 1 Sep 2026 16:09:25 +0100 Subject: [PATCH 7/7] tests: Tighten up MPI test slightly --- tests/test_mpi.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_mpi.py b/tests/test_mpi.py index 98e111ded1..443a804f34 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -1187,6 +1187,9 @@ def define(self, dimensions): # serial case. halo_types = (HaloUpdateCall, HaloUpdateList) if configuration['mpi'] else (HaloSpot,) + # Assert that HaloSpots etc are present, just not located within inner loops + assert len(FindNodes(halo_types).visit(op)) > 0 + for i in incr_iterations: assert len(FindNodes(halo_types).visit(i)) == 0