From 3c0f3ea183b7d7c05ccfe969346a5b67be2e2765 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Thu, 24 Sep 2026 13:35:30 +0100
Subject: [PATCH 1/6] mpi: Add raise_mpi utility for collective exceptions

---
 devito/mpi/distributed.py | 16 ++++++++++++++++
 tests/test_mpi.py         | 29 ++++++++++++++++++++++++++++-
 2 files changed, 44 insertions(+), 1 deletion(-)

diff --git a/devito/mpi/distributed.py b/devito/mpi/distributed.py
index 35596c39d64..2c3c4ad3641 100644
--- a/devito/mpi/distributed.py
+++ b/devito/mpi/distributed.py
@@ -70,6 +70,7 @@ def __getattr__(self, name):
     'SubDistributor',
     'devito_mpi_finalize',
     'devito_mpi_init',
+    'mpi_raise',
 ]
 
 
@@ -101,6 +102,21 @@ def devito_mpi_finalize():
         MPI.Finalize()
 
 
+def mpi_raise(error, exception=ValueError, comm=None):
+    """
+    Raise `exception` with the first non-None error message in rank order.
+
+    All ranks in `comm` must call this routine, including those with no local
+    error (`error=None`). This prevents a rank-local exception from stranding
+    peers in subsequent MPI calls. With no communicator or `MPI.COMM_NULL`,
+    only the local error is checked.
+    """
+    if comm is not None and comm is not MPI.COMM_NULL:
+        error = next((i for i in comm.allgather(error) if i is not None), None)
+    if error is not None:
+        raise exception(error)
+
+
 class AbstractDistributor(ABC):
 
     """
diff --git a/tests/test_mpi.py b/tests/test_mpi.py
index 6081fa627eb..ac3408d0c41 100644
--- a/tests/test_mpi.py
+++ b/tests/test_mpi.py
@@ -13,13 +13,14 @@
 )
 from devito.arch.compiler import OneapiCompiler
 from devito.data import LEFT, RIGHT
+from devito.exceptions import InvalidArgument
 from devito.ir import Cluster, Interval, IterationSpace
 from devito.ir.clusters.algorithms import check_halo_writes
 from devito.ir.iet import (
     Call, Conditional, FindNodes, FindSymbols, Iteration, retrieve_iteration_tree
 )
 from devito.ir.support.space import Backward, Forward
-from devito.mpi import MPI
+from devito.mpi import MPI, mpi_raise
 from devito.mpi.distributed import CustomTopology
 from devito.mpi.routines import ComputeCall, HaloUpdateCall, HaloUpdateList, MPICall
 from devito.tools import Bunch
@@ -28,6 +29,32 @@
 from examples.seismic.acoustic import acoustic_setup
 
 
+class TestMPIUtils:
+
+    @pytest.mark.parametrize('comm', [None, MPI.COMM_NULL], ids=['none', 'null'])
+    @pytest.mark.parametrize('error', [None, '', 'invalid argument'])
+    @pytest.mark.parametrize('exception', [ValueError, InvalidArgument])
+    def test_raise_serial(self, comm, error, exception):
+        if error is None:
+            mpi_raise(error, exception, comm=comm)
+        else:
+            with pytest.raises(exception, match=f'^{error}$'):
+                mpi_raise(error, exception, comm=comm)
+
+    @pytest.mark.parametrize('failing_ranks', [(), (0,), (1,), (0, 1)])
+    @pytest.mark.parallel(mode=[(2, 'basic')])
+    def test_raise_collective(self, failing_ranks, mode):
+        grid = Grid(shape=(16, 16))
+        rank = grid.distributor.myrank
+        error = f'rank {rank}' if rank in failing_ranks else None
+
+        if failing_ranks:
+            with pytest.raises(InvalidArgument, match=f'^rank {failing_ranks[0]}$'):
+                mpi_raise(error, InvalidArgument, comm=grid.comm)
+        else:
+            mpi_raise(error, InvalidArgument, comm=grid.comm)
+
+
 class TestDistributor:
 
     @pytest.mark.parallel(mode=[2, 4])

From 68f52ce6efc077c8d061b0296f9ed504138d8876 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Wed, 23 Sep 2026 10:24:05 +0100
Subject: [PATCH 2/6] compiler: Revamp SubDimension DDA

---
 devito/ir/support/basic.py         | 208 +++++++++++++-----------
 devito/operator/operator.py        |   2 +-
 devito/passes/iet/orchestration.py |  10 +-
 devito/types/dimension.py          |  42 ++++-
 devito/types/grid.py               |  23 ++-
 tests/test_ir.py                   | 252 +++++++++++++++++++++++++++--
 tests/test_operator.py             |   6 +-
 tests/test_subdomains.py           | 183 ++++++++++++++++++++-
 8 files changed, 609 insertions(+), 117 deletions(-)

diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py
index 518e1d0f442..034dbddfa63 100644
--- a/devito/ir/support/basic.py
+++ b/devito/ir/support/basic.py
@@ -7,7 +7,7 @@
 from sympy import Expr, S
 
 from devito.ir.support.space import Backward, null_ispace
-from devito.ir.support.utils import AccessMode, extrema
+from devito.ir.support.utils import AccessMode, erange, extrema
 from devito.ir.support.vector import LabeledVector, Vector
 from devito.symbolics import (
     compare_ops, q_affine, q_comp_acc, q_constant, retrieve_indexed, search
@@ -358,6 +358,9 @@ def distance(self, other, logical=False):
             # E.g., `uv(x).x` and `uv(x).y` -- not a real dependence!
             return Vector(S.ImaginaryUnit)
 
+        if disjoint_subdims(self, other):
+            return Vector(S.ImaginaryUnit)
+
         ret = []
         for sit, oit in zip(self.itintervals, other.itintervals, strict=False):
             n = len(ret)
@@ -369,20 +372,14 @@ def distance(self, other, logical=False):
                 # E.g., `self=R<f,[x]>` and `self.itintervals=(x, i)`
                 break
 
-            # If over SubDimensions, check disjointness
-            test = disjoint_subdims(self[n], other[n], sai, oai, sit, oit)
-            if test == DISJOINT:
-                return Vector(S.ImaginaryUnit)
-            elif test == MAYBE_OVERLAP:
-                ret.append(S.Infinity)
-                continue
-
             try:
                 if not (sit == oit and sai.root is oai.root):
                     # E.g., `self=R<f,[x + 2]>` and `other=W<f,[i + 1]>`
                     # E.g., `self=R<f,[x]>`, `other=W<f,[x + 1]>`,
                     #       `self.itintervals=(x<0>,)`, `other.itintervals=(x<1>,)`
-                    return vinf(ret)
+                    # Keep looking: a later axis may prove disjointness
+                    ret.append(S.Infinity)
+                    continue
             except AttributeError:
                 # E.g., `self=R<f,[cy]>` and `self.itintervals=(y,)` => `sai=None`
                 pass
@@ -1152,7 +1149,29 @@ def reads_smart_gen(self, f):
         """
         Generate all read accesses to a given function.
 
-        StencilDimensions, if any, are replaced with their extrema.
+        StencilDimensions, if any, are replaced with:
+
+            * in presence of SubDimensions: the range of points they span;
+            * in all other cases: just their extrema, since it suffices to
+              capture all possible dependencies.
+
+        The reason SubDimensions must be treated specially -- with a full set
+        of TimedAccess objects getting generated -- is to handle the special
+        case of slabs thinner than the stencil’s reach. For example, consider
+        the following scenario:
+
+        * A SubDimension with just two points, 10 and 11;
+        * One equation writes `F[10]` and `F[11]`;
+        * Another equation runs over the same SubDimension reading the stencil
+          `F[x-4] ... F[x+4]`.
+
+        If we examine only the two extreme stencil offsets:
+
+        * `F[x-4]` reads points 6–7: no overlap.
+        * `F[x+4]` reads points 14–15: no overlap.
+
+        But interior offsets certainly overlap -- for instance, `F[x-1]` reads
+        9–10, which includes the producer’s point 10.
 
         Notes
         -----
@@ -1163,9 +1182,13 @@ def reads_smart_gen(self, f):
         be found. For example, a DiscreteFunction would never appear among
         the iteration symbols.
         """
+        uses_subdims = lambda i: any(d.is_Sub for d in i.ispace.dimensions)
+
         if isinstance(f, (Function, Temp, TempArray, TBArray)):
             for i in self.getreads(f):
-                for j in extrema(i.access):
+                expand = erange if uses_subdims(i) else extrema
+
+                for j in expand(i.access):
                     yield TimedAccess(j, i.mode, i.timestamp, i.ispace)
 
         else:
@@ -1581,90 +1604,95 @@ def skippable_interval(d, ispace, it):
     return d is None or (d in ispace and not d._defines & it.dim._defines)
 
 
-# Possible return values for `disjoint_subdims`
-INAPPLICABLE = 0
-DISJOINT = 1
-MAYBE_OVERLAP = 2
-
-
-def disjoint_subdims(e0, e1, d0, d1, it0, it1):
+def disjoint_subdims(a0, a1):
     """
-    Determine whether two accesses span distinct pieces of the same
-    SubDimension decomposition.
-
-    Consider a root Dimension `x` with bounds `x_m` and `x_M`. A valid
-    left/middle/right decomposition with thicknesses `L` and `R` is::
-
-        xl = [x_m,         x_m + L - 1]
-        xm = [x_m + L,     x_M - R]
-        xr = [x_M - R + 1, x_M]
-
-    These intervals are pairwise disjoint. Replacing `xl`, `xm`, or `xr`
-    with `x` in an affine access removes the choice of partition piece while
-    retaining the relative access. If two such normalized accesses have zero
-    distance, they apply the same affine map to disjoint intervals and therefore
-    cannot refer to the same data point. The apparent dependence is imaginary.
-
-    For example, `f[xl]` and `f[xm]` normalize to `f[x]` and `f[x]`;
-    they are independent. The same holds for `f[xl + 1]` and `f[xm + 1]`
-    when their iteration intervals have equal offsets. By contrast, `f[xl]`
-    and `f[xm - 1]` normalize to different accesses, and the latter may reach
-    into the left piece, so they must be treated conservatively.
-
-    This proof requires distinct pieces of the same root, compatible declared
-    thicknesses, affine accesses, and iteration intervals with equal offsets and
-    directions. Runtime bounds are assumed to preserve the declared partition.
-    Return DISJOINT if disjointness is proven, and MAYBE_OVERLAP if the
-    intervals are aligned SubDimensions but are not proven disjoint. In
-    particular, two declarations of the same left, right, or middle piece
-    overlap along this Dimension. MAYBE_OVERLAP lets the caller record an
-    infinite distance and inspect later Dimensions, which may still prove the
-    multidimensional accesses disjoint. Return INAPPLICABLE if this test does not
-    apply, so that the general distance analysis can classify the dependence.
+    Determine whether two TimedAccesses touch disjoint SubDimension regions
+    of the same Function.
+
+    Compare symbolic accessed bounds, including shifts and stencil points.
+    Block intervals are promoted to their logical SubDimensions. Bounds and
+    thicknesses remain symbolic: MPI decomposition and runtime overrides can
+    change their values independently of the defaults.
+
+    For example, `xl = [m, m + L - 1]` and `xm = [m + L, M - R]` are
+    disjoint when they share the symbol `L`, whatever its runtime value.
+    Equal default thicknesses alone do not establish that relationship.
+
+    Opposite left/right slabs are assumed to form a valid partition: their
+    thicknesses satisfy `L + R <= N`. For translated stencil accesses, the
+    interior must also accommodate their combined inward reach. For example,
+    a pointwise left write and a right read at offset -4 require four interior
+    points. Runtime space_order checks cover explicit middle SubDimensions,
+    not arbitrary left/right pairs; no concrete domain size or thickness is
+    used here.
+
+    Match data axes independently of the iteration nests. Return True if any
+    axis proves separation, False otherwise. Accesses over the same interval
+    use the general distance analysis.
     """
-    try:
-        # E.g., `f[xl]` over `(xl,)` and `f[xm]` over `(xm,)` need this
-        # special test, while accesses over the same `(xl,)` should use general
-        # distance analysis, so we can return immediately in such a case
-        if not (d0.is_Sub and
-                d1.is_Sub and
-                d0.root is d1.root and
-                it0.dim.root is d0.root and
-                it1.dim.root is d1.root and
+    for e0, e1, d0, d1 in zip(a0, a1, a0.aindices, a1.aindices, strict=False):
+        it0 = a0.intervals[d0]
+        it1 = a1.intervals[d1]
+        if it0.is_Null or it1.is_Null:
+            continue
+
+        it0 = it0.promote(lambda d: d.is_Incr)
+        it1 = it1.promote(lambda d: d.is_Incr)
+        if not (it0.dim.is_Sub and
+                it1.dim.is_Sub and
+                it0.dim.root is it1.dim.root and
                 it0 != it1):
-            return INAPPLICABLE
-    except AttributeError:
-        return INAPPLICABLE
-
-    if (d0.is_left and d1.is_middle) or \
-       (d0.is_middle and d1.is_left):
-        is_partition = d0.ltkn.value == d1.ltkn.value
-    elif (d0.is_middle and d1.is_right) or \
-         (d0.is_right and d1.is_middle):
-        is_partition = d0.rtkn.value == d1.rtkn.value
-    elif d0.is_left and d1.is_right:
-        is_partition = d0.ltkn.value is not None and d1.rtkn.value is not None
-    elif d0.is_right and d1.is_left:
-        is_partition = d0.rtkn.value is not None and d1.ltkn.value is not None
-    else:
-        is_partition = False
-
-    if not is_partition:
-        return MAYBE_OVERLAP
-
-    if not q_affine(e0, d0) or not q_affine(e1, d1):
-        return MAYBE_OVERLAP
+            continue
 
-    if it0.offsets != it1.offsets or it0.direction is not it1.direction:
-        return MAYBE_OVERLAP
+        bounds = []
+        for e, d, it in ((e0, d0, it0), (e1, d1, it1)):
+            if not q_affine(e, d):
+                break
 
-    e0 = e0._subs(d0, d0.root)
-    e1 = e1._subs(d1, d1.root)
+            lower, upper = [], []
+            for v in erange(e):
+                slope = v.diff(d)
+                if slope.is_nonnegative:
+                    m, M = it.symbolic_min, it.symbolic_max
+                elif slope.is_nonpositive:
+                    M, m = it.symbolic_min, it.symbolic_max
+                else:
+                    break
+                lower.append(v._subs(d, m))
+                upper.append(v._subs(d, M))
+            else:
+                bounds.append((sympy.Min(*lower), sympy.Max(*upper)))
+
+        if len(bounds) == 2:
+            (m0, M0), (m1, M1) = bounds
+            mapper = {}
+
+            dl, dr = (it0.dim, it1.dim) if it0.dim.is_left else (it1.dim, it0.dim)
+            dlp, drp = dl.parent, dr.parent
+
+            if dl.is_left and dr.is_right and dlp is drp:
+                # A valid partition satisfies L + R <= N, where N is the parent
+                # extent; an explicit middle SubDomain checks this at construction.
+                # Further, for stencils, we require that:
+                # `N - L - R >= the combined inward reach`
+                # so accesses from opposite slabs cannot meet. Explicit middle
+                # SubDimensions check for at least space_order interior points
+                # at *op.apply time*, accounting for runtime overrides.  Without
+                # an explicit middle, the gap assumption is unchecked
+                gap = sympy.Dummy(nonnegative=True)
+                if e0.diff(d0) == e1.diff(d1) == 1:
+                    M, m = (M0, m1) if it0.dim.is_left else (M1, m0)
+                    reach = (M - dl.symbolic_max - m + dr.symbolic_min).expand()
+                    if is_integer(reach):
+                        gap += max(0, reach)
+
+                mapper[dlp.symbolic_max] = dlp.symbolic_min + dl.ltkn + dr.rtkn + gap - 1
+
+            if (M0 - m1).subs(mapper).is_negative or \
+               (M1 - m0).subs(mapper).is_negative:
+                return True
 
-    if e0 - e1 == 0:
-        return DISJOINT
-    else:
-        return MAYBE_OVERLAP
+    return False
 
 
 def disjoint_test(e0, e1, d, it):
diff --git a/devito/operator/operator.py b/devito/operator/operator.py
index c97870f075a..875eb584a7f 100644
--- a/devito/operator/operator.py
+++ b/devito/operator/operator.py
@@ -716,7 +716,7 @@ def _prepare_arguments(self, autotune=None, estimate_memory=False, **kwargs):
             except AttributeError:
                 pass
             if d.is_Derived:
-                d._arg_check(args)
+                d._arg_check(args, **kwargs)
 
         # Turn arguments into a format suitable for the generated code
         # E.g., instead of NumPy arrays for Functions, the generated code expects
diff --git a/devito/passes/iet/orchestration.py b/devito/passes/iet/orchestration.py
index f9969801617..443a6b78ef7 100644
--- a/devito/passes/iet/orchestration.py
+++ b/devito/passes/iet/orchestration.py
@@ -125,9 +125,9 @@ def _make_syncarray(self, iet, sync_ops, layer):
     def _make_prefetchupdate(self, iet, sync_ops, layer, wrap=True):
         return self._make_async_task(prefetchupdate, iet, sync_ops, layer, wrap)
 
-    @iet_pass
-    def process(self, iet):
-        callbacks = {
+    @property
+    def _callbacks(self):
+        return {
             WaitLock: self._make_waitlock,
             WithLock: self._make_withlock,
             SyncArray: self._make_syncarray,
@@ -139,13 +139,15 @@ def process(self, iet):
             AsyncCallable: self._make_async_callable
         }
 
+    @iet_pass
+    def process(self, iet):
         # Collect the compatible asynchronous task groups, if any
         task_groups = TaskGroups()
         if self.npthreads:
             CollectTasks(task_groups).visit(iet)
 
         # Lower the SyncSpots in a single bottom-up traversal, atomically lowering
-        lowerer = LowerSyncSpots(callbacks, task_groups, self.sregistry)
+        lowerer = LowerSyncSpots(self._callbacks, task_groups, self.sregistry)
         iet = lowerer.visit(iet)
 
         return iet, {'efuncs': lowerer.efuncs}
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index fec9fb78a33..c131cd39d06 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -818,6 +818,46 @@ def _arg_values(self, interval, grid=None, **kwargs):
         # themselves
         return {}
 
+    def _arg_check(self, args, *_args, **kwargs):
+        # These modules depend on Dimension, so importing them above would cycle
+        from devito.mpi import mpi_raise  # noqa: PLC0415
+        from devito.symbolics import subs_op_args  # noqa: PLC0415
+
+        if not self.is_middle:
+            return
+
+        # Function._arg_check visits original axes (e.g. x_ltkn), whereas `args`
+        # contains the concretized thicknesses (x_ltkn0, ...). The Operator checks
+        # the matching concrete SubDimensions separately in its dimension loop
+        if self not in args.op.dimensions:
+            return
+
+        d = self.root
+        if args.grid is not None and args.grid.is_distributed(d):
+            # Check the global runtime region: its MPI-local slices may be empty
+            # or smaller than space_order even for a non-degenerate global interior
+            size = args.grid.size_map[d].glb
+            values = {**args,
+                      d.min_name: kwargs.get(d.min_name, 0),
+                      d.max_name: kwargs.get(d.max_name, kwargs.get(d.name, size - 1)),
+                      **{t.name: kwargs.get(t.name, t.value) for t in self.thickness}}
+        else:
+            values = args
+        size = int(subs_op_args(self.symbolic_size, values))
+
+        # Runtime overrides do not change the compiled stencil order
+        items = [f.space_order for f in args.op.input if f.is_DiscreteFunction]
+        space_order = max(items, default=0)
+
+        if size < space_order:
+            error = (f"Expected at least {space_order} interior points along "
+                     f"`{self.parent}` (space_order), but runtime arguments leave {size}")
+        else:
+            error = None
+
+        comm = args.comm if args.options['mpi'] else None
+        mpi_raise(error, InvalidArgument, comm=comm)
+
 
 class MultiSubDimension(AbstractSubDimension):
 
@@ -1436,7 +1476,7 @@ def _arg_values(self, interval, grid=None, args=None, **kwargs):
                 # Avoid OOB (will end up here only in case of tiny iteration spaces)
                 return {name: 1}
 
-    def _arg_check(self, args, *_args):
+    def _arg_check(self, args, *_args, **kwargs):
         try:
             name = self.step.name
         except AttributeError:
diff --git a/devito/types/grid.py b/devito/types/grid.py
index 5884b0fc30d..ae6df92cb49 100644
--- a/devito/types/grid.py
+++ b/devito/types/grid.py
@@ -641,12 +641,6 @@ def __subdomain_finalize_legacy__(self, grid):
                 try:
                     # Case ('middle', int, int)
                     side, ltkn, rtkn = v
-                    if side != 'middle':
-                        raise ValueError(f"Expected side 'middle', not `{side}`")
-                    sub_dimensions.append(SubDimension.middle(f'i{k.name}',
-                                                              k, ltkn, rtkn))
-                    thickness = s-ltkn-rtkn
-                    sdshape.append(thickness)
                 except ValueError:
                     side, thickness = v
                     constructor = {'left': SubDimension.left,
@@ -663,6 +657,23 @@ def __subdomain_finalize_legacy__(self, grid):
                         ) from None
                     sub_dimensions.append(constructor(f'i{k.name}', k, thickness))
                     sdshape.append(thickness)
+                else:
+                    if side != 'middle':
+                        raise ValueError(f"Expected side 'middle', not `{side}`")
+
+                    # A `middle` region expects `ltkn + rtkn <= s` in the global Grid.
+                    # This ensures that the left and right regions won't overlap
+                    thickness = s-ltkn-rtkn
+                    if thickness < 0:
+                        raise ValueError(
+                            f"SubDomain `{self.name}` has combined thickness "
+                            f"{ltkn + rtkn} along `{k}`, exceeding the Grid size {s}"
+                        )
+
+                    sub_dimensions.append(
+                        SubDimension.middle(f'i{k.name}', k, ltkn, rtkn)
+                    )
+                    sdshape.append(thickness)
 
         self._shape = tuple(sdshape)
         self._dimensions = tuple(sub_dimensions)
diff --git a/tests/test_ir.py b/tests/test_ir.py
index 1757fe29a7f..63a462b1a91 100644
--- a/tests/test_ir.py
+++ b/tests/test_ir.py
@@ -15,8 +15,8 @@
 from devito.ir.iet import FindNodes, Iteration
 from devito.ir.stree import stree_build
 from devito.ir.support.basic import (
-    AFFINE, IRREGULAR, REGULAR, IterationInstance, Scope, TimedAccess, Vector, mocksym0,
-    mocksym1
+    AFFINE, IRREGULAR, REGULAR, IterationInstance, Relation, Scope, TimedAccess, Vector,
+    mocksym0, mocksym1
 )
 from devito.ir.support.guards import GuardOverflow
 from devito.ir.support.space import (
@@ -28,8 +28,8 @@
 from devito.tools import prod
 from devito.tools.data_structures import frozendict
 from devito.types import (
-    Array, Bundle, ConditionalDimension, CriticalRegion, CustomDimension, Jump, Scalar,
-    Symbol
+    Array, BlockDimension, Bundle, ConditionalDimension, CriticalRegion, CustomDimension,
+    Jump, Scalar, StencilDimension, Symbol
 )
 
 
@@ -316,9 +316,11 @@ def test_timed_access_distance(self, x, y, ta_literal):
         assert tcxy_r0.distance(tcx1y1_r1) == (-1, -1)
         assert tcx1y1_r1.distance(tcx1y_r1) == (0, 1)
 
-        # Distance should go to infinity due to mismatching directions
-        assert rev_tcxy_w0.distance(tcx1y_r1) == (S.Infinity,)
-        assert tcx1y_r1.distance(rev_tcxy_w0) == (S.Infinity,)
+        # Mismatching x directions do not prevent computing the y distance
+        assert rev_tcxy_w0.distance(tcx1y_r1) == (S.Infinity, 0)
+        assert tcx1y_r1.distance(rev_tcxy_w0) == (S.Infinity, 0)
+        assert rev_tcxy_w0.distance(tcx1y1_r1) == (S.Infinity, -1)
+        assert tcx1y1_r1.distance(rev_tcxy_w0) == (S.Infinity, 1)
 
         # Distance when both source and since go backwards along the x Dimension
         assert rev_tcxy_w0.distance(rev_tcx1y1_r1) == (1, -1)
@@ -409,6 +411,10 @@ def test_timed_access_distance_subdimensions(self):
         xm_overlap = SubDimension.middle('xm_overlap', x, 4, 4)
         xr_overlap = SubDimension.right('xr_overlap', x, 4)
 
+        # Adjacent pieces share the runtime boundary, not just its default value
+        xm = xm._rebuild(thickness=(xl.ltkn, xr.rtkn))
+        ym = ym._rebuild(thickness=(yl.ltkn, ym.rtkn))
+
         f = Function(name='f', grid=grid)
 
         left = TimedAccess(
@@ -471,10 +477,149 @@ def test_timed_access_distance_subdimensions(self):
         assert middle.distance(middle_overlap) == (S.Infinity, 0)
         assert right.distance(right_overlap) == (S.Infinity, 0)
         assert left.distance(bad) == (S.Infinity, 0)
-        assert left.distance(shifted) == (S.Infinity, 0)
-        assert left.distance(shifted_range) == (S.Infinity, 0)
+        assert left.distance(shifted) == (S.ImaginaryUnit,)
+        assert left.distance(shifted_range) == (S.ImaginaryUnit,)
         assert left_nonlinear.distance(middle_nonlinear) == (S.Infinity, 0)
-        assert left.distance(orthogonal) == (S.Infinity,)
+        assert left.distance(orthogonal) == (S.Infinity, S.Infinity)
+
+    @pytest.mark.parametrize('shared_boundary', [False, True])
+    @pytest.mark.parametrize('offset,independent', [(-5, False), (-4, True), (1, True)])
+    def test_subdimension_stencil_distance(self, shared_boundary, offset, independent):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        xl = SubDimension.left('xl', x, 8)
+        xr = SubDimension.right('xr', x, 20)
+        if shared_boundary:
+            xr = SubDimension.middle('xm', x, 8, 0)
+            xr = xr._rebuild(thickness=(xl.ltkn, xr.rtkn))
+            interval = Interval(xr, 4, 4)
+        else:
+            interval = Interval(xr)
+        h = StencilDimension('h', 0, 8)
+        f = Function(name='f', grid=grid)
+        a = TimedAccess(f[xl], 'W', 0, IterationSpace([Interval(xl)]))
+        b = TimedAccess(f[xr + offset + h], 'R', 1,
+                        IterationSpace([interval]))
+        independent = independent if shared_boundary else True
+        assert (S.ImaginaryUnit in a.distance(b)) is independent
+        assert (S.ImaginaryUnit in b.distance(a)) is independent
+
+    @pytest.mark.parametrize('direction', [Forward, Backward])
+    @pytest.mark.parametrize('blocked', [False, True])
+    @pytest.mark.parametrize('side,thickness,shift,expected', [
+        ('right', 20, 0, (S.ImaginaryUnit,)),
+        ('right', 24, 0, (S.ImaginaryUnit,)),
+        ('right', 20, -8, (S.ImaginaryUnit,)),
+        ('middle', 8, 0, (S.ImaginaryUnit,)),
+        ('middle', 8, -1, (S.Infinity, S.Infinity))
+    ])
+    def test_subdimension_distance_different_nests(self, direction, blocked, side,
+                                                   thickness, shift, expected):
+        grid = Grid(shape=(32, 32))
+        x, y = grid.dimensions
+        yl = SubDimension.left('yl', y, 8)
+        if side == 'middle':
+            yr = SubDimension.middle('ym', y, thickness, 0)
+            yr = yr._rebuild(thickness=(yl.ltkn, yr.rtkn))
+        else:
+            yr = SubDimension.right('yr', y, thickness)
+        i = Dimension(name='i')
+        f = Function(name='f', grid=grid)
+        if blocked:
+            yl = BlockDimension('ylb', yl, yl.symbolic_min, yl.symbolic_max, step=1)
+            yr = BlockDimension('yrb', yr, yr.symbolic_min, yr.symbolic_max, step=1)
+
+        a = TimedAccess(f[x, yl], 'W', 0, IterationSpace([Interval(x), Interval(yl)]))
+        b = TimedAccess(f[x, yr], 'R', 1,
+                        IterationSpace([Interval(i), Interval(x),
+                                        Interval(yr, shift, shift)],
+                                       directions={yr: direction}))
+
+        # The first iteration intervals differ, but y may still prove disjointness
+        assert a.distance(b) == b.distance(a) == expected
+
+    def test_indexedbase_distance(self):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        f = Function(name='f', grid=grid)
+        a = TimedAccess(f.indexed, 'R', 0)
+        b = TimedAccess(f[x], 'W', 1, IterationSpace([Interval(x)]))
+        assert a.distance(b) == b.distance(a) == (S.Infinity,)
+
+    @pytest.mark.parametrize('shared_boundary', [False, True])
+    @pytest.mark.parametrize('symbolic', [False, True])
+    @pytest.mark.parametrize('slope,offset,expected', [
+        (1, 0, S.ImaginaryUnit), (1, -1, S.Infinity),
+        (-1, 0, S.ImaginaryUnit), (-1, 1, S.Infinity),
+        (2, 0, S.ImaginaryUnit), (2, -2, S.Infinity),
+        (1, None, S.Infinity)
+    ])
+    def test_subdimension_affine_bounds(self, shared_boundary, symbolic, slope, offset,
+                                        expected):
+        grid = Grid(shape=(64,))
+        x, = grid.dimensions
+        xl = SubDimension.left('xl', x, 8)
+        xm = SubDimension.middle('xm', x, 8, 40)
+        if shared_boundary:
+            xm = xm._rebuild(thickness=(xl.ltkn, xm.rtkn))
+        else:
+            expected = S.Infinity
+        if symbolic:
+            f = Array(name='f', dimensions=(x,))
+        else:
+            f = Function(name='f', grid=grid)
+        base = 0 if slope > 0 else 63
+        offset = Symbol(name='offset', integer=True) if offset is None else offset
+
+        a = TimedAccess(f[base + slope*xl], 'W', 0,
+                        IterationSpace([Interval(xl)]))
+        b = TimedAccess(f[base + slope*xm + offset], 'R', 1,
+                        IterationSpace([Interval(xm)]))
+        assert a.distance(b) == b.distance(a) == (expected,)
+
+    @pytest.mark.parametrize('index', [
+        lambda d: d % 2,
+        lambda d: Symbol(name='s', integer=True)*d
+    ])
+    @pytest.mark.parametrize('shift,expected', [
+        (0, (S.ImaginaryUnit,)), (-1, (S.Infinity, S.Infinity))
+    ])
+    def test_subdimension_disjoint_later_axis(self, index, shift, expected):
+        grid = Grid(shape=(32, 32))
+        x, y = grid.dimensions
+        xl = SubDimension.left('xl', x, 8)
+        xr = SubDimension.right('xr', x, 20)
+        yl = SubDimension.left('yl', y, 8)
+        ym = SubDimension.middle('ym', y, 8, 0)
+        ym = ym._rebuild(thickness=(yl.ltkn, ym.rtkn))
+        f = Function(name='f', grid=grid)
+
+        a = TimedAccess(f[index(xl), yl], 'W', 0,
+                        IterationSpace([Interval(xl), Interval(yl)]))
+        b = TimedAccess(f[index(xr), ym], 'R', 1,
+                        IterationSpace([Interval(xr), Interval(ym, shift, shift)]))
+
+        # An unresolved x axis must not prevent y from proving disjointness
+        assert a.distance(b) == b.distance(a) == expected
+
+    @pytest.mark.parametrize('slope', [1, -1, 2])
+    @pytest.mark.parametrize('offset,expected', [
+        (0, S.ImaginaryUnit), (1, S.ImaginaryUnit), (-1, S.Infinity)
+    ])
+    def test_opposite_subdimension_bounds(self, slope, offset, expected):
+        x = Dimension(name='x')
+        xl = SubDimension.left('xl', x, 8)
+        xr = SubDimension.right('xr', x, 8)
+        f = Array(name='f', dimensions=(x,))
+
+        a = TimedAccess(f[slope*xl], 'W', 0, IterationSpace([Interval(xl)]))
+        b = TimedAccess(f[slope*(xr + offset)], 'R', 1,
+                        IterationSpace([Interval(xr)]))
+
+        # Translated accesses assume the interior accommodates their inward reach
+        if slope == 1:
+            expected = S.ImaginaryUnit
+        assert a.distance(b) == b.distance(a) == (expected,)
 
 
 class TestSpace:
@@ -1158,6 +1303,93 @@ def test_bundle_components(self):
         dep, = scope.d_flow
         assert dep.function is f
 
+    @pytest.mark.parametrize('symbolic', [False, True])
+    @pytest.mark.parametrize('lower,upper,ndeps', [
+        (0, 31, (8, 32)), (0, 7, (0, 8)), (16, 31, (0, 16)),
+        (7, 8, (1, 2)), (15, 16, (1, 2))
+    ])
+    def test_stencil_contains_producer(self, symbolic, lower, upper, ndeps):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        xl = SubDimension.left('xl', x, 1)
+        xm = SubDimension.middle('xm', x, 8, 16)
+        h = StencilDimension('h', lower, upper)
+        f = Function(name='f', grid=grid)
+        g = Function(name='g', grid=grid)
+        if symbolic:
+            pi, ci = Interval(xm), Interval(xl)
+        else:
+            # Encode actual fixed iteration bounds, not runtime defaults
+            pi = Interval(xm, 8 - xm.symbolic_min, 15 - xm.symbolic_max)
+            ci = Interval(xl, -xl.symbolic_min, -xl.symbolic_max)
+        producer = Cluster(Eq(f[xm], 1), IterationSpace([pi]))
+        consumer = Cluster(Eq(g[xl], f[xl + h]), IterationSpace([ci]))
+        ndeps = ndeps[symbolic]
+
+        # In particular, [0, 31] covers the producer despite both endpoints missing it
+        scope = Scope.from_scopes(producer.scope, consumer.scope)
+        w, = scope.getwrites(f)
+        r, = scope.getreads(f)
+        for relation in (Relation(w, r), Relation(r, w)):
+            assert relation.distance == ((S.Infinity if ndeps else S.ImaginaryUnit),)
+        assert len(scope.d_flow) == ndeps
+        assert len(scope.d_anti) == ndeps
+
+    @pytest.mark.parametrize('lower,upper,flow,anti', [
+        (-2, -1, {(1,), (2,)}, set()),
+        (-2, 2, {(2,)}, {(2,)}),
+        (0, 2, {(0,)}, {(2,)}),
+        (1, 2, set(), {(1,), (2,)})
+    ])
+    def test_stencil_same_domain(self, lower, upper, flow, anti):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        h = StencilDimension('h', lower, upper)
+        f = Function(name='f', grid=grid)
+        g = Function(name='g', grid=grid)
+        ispace = IterationSpace([Interval(x)])
+        producer = Cluster(Eq(f[x], 1), ispace)
+        consumer = Cluster(Eq(g[x], f[x + h]), ispace)
+
+        scope = Scope.from_scopes(producer.scope, consumer.scope)
+        assert {tuple(d.distance) for d in scope.d_flow} == flow
+        assert {tuple(d.distance) for d in scope.d_anti} == anti
+
+    @pytest.mark.parametrize('lower,upper,flow,anti', [
+        (-16, 16, True, True), (-16, 0, True, False), (0, 16, True, True),
+        (-16, -8, False, False), (8, 16, False, False)
+    ])
+    def test_stencil_same_subdimension(self, lower, upper, flow, anti):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        xl = SubDimension.left('xl', x, 8)
+        h = StencilDimension('h', lower, upper)
+        f = Function(name='f', grid=grid)
+        g = Function(name='g', grid=grid)
+        ispace = IterationSpace([Interval(xl)])
+        producer = Cluster(Eq(f[xl], 1), ispace)
+        consumer = Cluster(Eq(g[xl], f[xl + h]), ispace)
+
+        # Even in the same domain, stencil endpoints may lie beyond the producer
+        scope = Scope.from_scopes(producer.scope, consumer.scope)
+        assert bool(scope.d_flow) is flow
+        assert bool(scope.d_anti) is anti
+
+    def test_stencil_multidimensional_distance(self):
+        grid = Grid(shape=(32, 32))
+        x, y = grid.dimensions
+        h = StencilDimension('h', 0, 1)
+        f = Function(name='f', grid=grid)
+        g = Function(name='g', grid=grid)
+        ispace = IterationSpace([Interval(x), Interval(y)])
+        producer = Cluster(Eq(f[x, y], 1), ispace)
+        consumer = Cluster(Eq(g[x, y], f[x + h, y + 2 - h]), ispace)
+
+        scope = Scope.from_scopes(producer.scope, consumer.scope)
+        assert not scope.d_flow
+        assert {tuple(d.distance) for d in scope.d_anti} == {(0, 2), (1, 1)}
+        assert {d.cause for d in scope.d_anti} == {frozenset({x}), frozenset({y})}
+
 
 class TestParallelismAnalysis:
 
diff --git a/tests/test_operator.py b/tests/test_operator.py
index 59cf0013c78..4518b324299 100644
--- a/tests/test_operator.py
+++ b/tests/test_operator.py
@@ -18,9 +18,9 @@
 
 from devito import (  # noqa
     CELL, NODE, Buffer, CondEq, Constant, Dimension, Eq, Function, Ge, Grid, Gt, Inc, Le,
-    Lt, Operator, SpaceDimension, SparseFunction, SparseTimeFunction, TensorFunction,
-    TensorTimeFunction, TimeFunction, VectorFunction, VectorTimeFunction, configuration,
-    dimensions, div, error, exp, grad, sin, switchconfig
+    Lt, Operator, SpaceDimension, SparseFunction, SparseTimeFunction, SubDimension,
+    TensorFunction, TensorTimeFunction, TimeFunction, VectorFunction, VectorTimeFunction,
+    configuration, dimensions, div, error, exp, grad, sin, switchconfig
 )
 from devito.arch.archinfo import Device
 from devito.exceptions import InvalidOperator
diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py
index e6cab673116..0fbb47c5666 100644
--- a/tests/test_subdomains.py
+++ b/tests/test_subdomains.py
@@ -7,9 +7,10 @@
 from conftest import assert_structure, opts_tiling
 from devito import (
     Border, Buffer, ConditionalDimension, Constant, Eq, Function, Grid, Lt, Operator,
-    SparseFunction, SparseTimeFunction, SubDomain, SubDomainSet, TensorFunction,
-    TimeFunction, VectorFunction, solve
+    SparseFunction, SparseTimeFunction, SubDimension, SubDomain, SubDomainSet,
+    TensorFunction, TimeFunction, VectorFunction, solve
 )
+from devito.exceptions import InvalidArgument
 from devito.ir import (
     Expression, FindNodes, FindSymbols, Iteration, SymbolRegistry, retrieve_iteration_tree
 )
@@ -121,6 +122,29 @@ def define(self, dimensions):
         assert s_d1.shape == (4, 2)
         assert s_d2.shape == (3, 7)
 
+    @pytest.mark.parametrize('legacy', [False, True])
+    @pytest.mark.parametrize('size', [15, 16, 17])
+    @pytest.mark.parametrize('spec', [('left', 16), ('middle', 8, 8), ('right', 16)])
+    def test_partition_thickness(self, legacy, size, spec):
+        class Region(SubDomain):
+            name = 'region'
+
+            def define(self, dimensions):
+                x, = dimensions
+                return {x: spec}
+
+        def make_region():
+            if legacy:
+                return Grid(shape=(size,), subdomains=(Region(),)).subdomains['region']
+            return Region(grid=Grid(shape=(size,)))
+
+        if size < 16:
+            with pytest.raises(ValueError, match='thickness'):
+                make_region()
+        else:
+            expected = size - 16 if spec[0] == 'middle' else 16
+            assert make_region().shape == (expected,)
+
     def test_definitions(self):
 
         class sd0(SubDomain):
@@ -1908,3 +1932,158 @@ def define(self, dimensions):
         eq = Eq(g, g + f.dx)
         eqe = eq.evaluate
         assert eqe.rhs == g + f.dx(x0=x).evaluate._subs(x, g.dimensions[1])
+
+
+class TestSubDomainArguments:
+
+    @staticmethod
+    def _make_operator(left_shift=0, right_shift=0, grid=None, middle=False):
+        grid = grid or Grid(shape=(8, 32))
+        y = grid.dimensions[-1]
+
+        yl = SubDimension.left('yl', y, 8)
+        yr = SubDimension.right('yr', y, 8)
+        if middle:
+            yl = yr = SubDimension.middle('ym', y, 8, 8)
+
+        u = TimeFunction(name='u', grid=grid, space_order=8)
+        v = TimeFunction(name='v', grid=grid, space_order=8)
+
+        eqs = [
+            Eq(u.forward.subs(y, y + left_shift), 1).subs(y, yl),
+            Eq(v.forward, u.forward.subs(y, y + right_shift) + 1).subs(y, yr)
+        ]
+
+        op = Operator(eqs, name='subdomain_arguments')
+
+        return op, (u, v)
+
+    @pytest.mark.parametrize('middle', [False, True])
+    @pytest.mark.parametrize('left_shift,right_shift', [
+        (0, 0), (0, -4), (2, -4)
+    ])
+    @pytest.mark.parametrize('margin', [-1, 0, 1])
+    @pytest.mark.parametrize('override', ['function', 'bounds', 'thickness'])
+    def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
+        op, (u, v) = self._make_operator(left_shift, right_shift, middle=middle)
+
+        size = 16 + 8 + margin
+        if override == 'function':
+            grid = Grid(shape=(8, size))
+
+            kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
+                                           space_order=8) for f in (u, v)}
+        elif override == 'bounds':
+            kwargs = {'y_m': 3, 'y_M': 3 + size - 1}
+        else:
+            dl, = [d for d in op.dimensions if d.is_Sub and not d.is_right]
+            kwargs = {dl.ltkn.name: 32 - 8 - 8 - margin}
+
+        if middle and margin < 0:
+            with pytest.raises(InvalidArgument, match='at least 8 interior points'):
+                op.arguments(time_M=0, **kwargs)
+        else:
+            op.arguments(time_M=0, **kwargs)
+
+    @pytest.mark.parametrize('middle', [False, True])
+    @pytest.mark.parametrize('space_order', [4, 8, 12])
+    @pytest.mark.parametrize('margin', [-1, 0, 1])
+    def test_runtime_space_order(self, middle, space_order, margin):
+        """Override metadata does not change the compiled stencil order."""
+        op, fields = self._make_operator(middle=middle)
+
+        required = 8
+        grid = Grid(shape=(8, 16 + required + margin))
+
+        kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
+                                       space_order=space_order) for f in fields}
+
+        if middle and margin < 0:
+            with pytest.raises(InvalidArgument,
+                               match=f'at least {required} interior points'):
+                op.arguments(time_M=0, **kwargs)
+        else:
+            op.arguments(time_M=0, **kwargs)
+
+    @pytest.mark.parametrize('side', ['left', 'right'])
+    def test_empty_slab(self, side):
+        op, _ = self._make_operator(right_shift=-4)
+
+        d, = [d for d in op.dimensions if d.is_Sub and getattr(d, f'is_{side}')]
+        thickness = d.ltkn if side == 'left' else d.rtkn
+
+        # The active slab fills the local domain; the opposite slab is absent
+        op.arguments(time_M=0, y_M=7, **{thickness.name: 0})
+
+    def test_before_autotuning(self):
+        op, _ = self._make_operator(right_shift=-4, middle=True)
+
+        with pytest.raises(InvalidArgument, match='interior points'):
+            op.arguments(time_M=0, y_M=22, autotune=True)
+
+        assert 'autotuning' not in op._state
+
+        op.arguments(time_M=0, y_M=23, autotune=True)
+
+        assert len(op._state['autotuning']) == 1
+
+    @pytest.mark.parametrize('left', [15, 16, 17, 25])
+    @pytest.mark.parallel(mode=[(2, 'basic')])
+    def test_distributed_middle(self, left, mode):
+        grid = Grid(shape=(16, 32), topology=(1, 2))
+
+        op, _ = self._make_operator(grid=grid, middle=True)
+
+        d, = [d for d in op.dimensions if d.is_Sub]
+
+        # Rank 0 has one point or an empty middle (possibly with inverted bounds).
+        # Only the global size determines whether the middle is large enough
+        kwargs = {d.ltkn.name: left, d.rtkn.name: 0}
+
+        if left == 25:
+            with pytest.raises(InvalidArgument, match='at least 8 interior points'):
+                op.arguments(time_M=0, **kwargs)
+        else:
+            op.arguments(time_M=0, **kwargs)
+
+    @pytest.mark.parallel(mode=[(2, 'basic')])
+    def test_collective_rejection(self, mode):
+        grid = Grid(shape=(32, 32), topology=(2, 1))
+
+        op, _ = self._make_operator(right_shift=-4, grid=grid, middle=True)
+
+        dl, = [d for d in op.dimensions if d.is_Sub and d.is_middle]
+        left = 24 if grid.distributor.myrank == 0 else 0
+
+        # Only rank 0 has an insufficient middle; its peer has 24 interior points
+        with pytest.raises(InvalidArgument, match='interior points'):
+            op.arguments(time_M=0, **{dl.ltkn.name: left})
+
+    def test_function_on_subdomain(self):
+        class Interior(SubDomain):
+
+            def define(self, dimensions):
+                x, y = dimensions
+                return {x: x, y: ('middle', 8, 8)}
+
+        grid = Grid(shape=(16, 32))
+
+        f = Function(name='f', grid=Interior(grid=grid), space_order=8)
+        original = f.dimensions[-1]
+
+        eq = Eq(f, f + 1)
+
+        op = Operator(eq, name='subdomain_function_arguments')
+
+        concrete, = [d for d in op.dimensions if d.is_Sub]
+
+        # Function validation visits `original`; Operator validation visits `concrete`
+        assert original not in op.dimensions
+
+        args = op.arguments()
+
+        assert original.ltkn.name not in args
+        assert concrete.ltkn.name in args
+
+        with pytest.raises(InvalidArgument, match='at least 8 interior points'):
+            op.arguments(**{concrete.ltkn.name: 17})

From 91540c5c935b196877e8c8f779d1b42e6f884cbd Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Thu, 24 Sep 2026 15:42:16 +0100
Subject: [PATCH 3/6] tests: Update input grid to honour SubDomain rules

---
 tests/test_builtins.py | 3 ++-
 tests/test_dse.py      | 6 +++---
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/tests/test_builtins.py b/tests/test_builtins.py
index 275ee4bbe24..c3a86ea1cf4 100644
--- a/tests/test_builtins.py
+++ b/tests/test_builtins.py
@@ -147,7 +147,8 @@ def test_gs_1d_int(self, sigma):
     def test_gs_1d_float(self, sigma):
         """Test the Gaussian smoother in 1d on array of float."""
 
-        a = np.array([1.2, 2.7, 3.9, 4.1, 5.2, 6.5, 7.1, 9.3, 11.0])
+        a = np.array([1.2, 2.7, 3.9, 4.1, 5.2, 6.5, 7.1, 9.3, 11.0,
+                      8.7, 6.4, 4.6, 3.1, 2.9, 1.4, 0.8, 2.3])
         sp_smoothed = gaussian_filter(a, sigma=sigma)
         dv_smoothed = gaussian_smooth(a, sigma=sigma)
 
diff --git a/tests/test_dse.py b/tests/test_dse.py
index 2c75f60efcb..a0221244368 100644
--- a/tests/test_dse.py
+++ b/tests/test_dse.py
@@ -586,7 +586,7 @@ def test_full_shape_w_subdims(self, rotate):
         """
         Like `test_full_shape`, but SubDomains (and therefore SubDimensions) are used.
         """
-        grid = Grid(shape=(3, 3, 3))
+        grid = Grid(shape=(5, 5, 5))
         x, y, z = grid.dimensions
         t = grid.stepping_dim
 
@@ -753,12 +753,12 @@ def test_mixed_shapes_v2_w_subdims(self, rotate):
         Analogous `test_mixed_shapes`, but with different sets of aliasing expressions.
         Also, uses SubDimensions.
         """
-        grid = Grid(shape=(3, 3, 3))
+        grid = Grid(shape=(5, 5, 5))
         x, y, z = grid.dimensions
         t = grid.stepping_dim
         d = Dimension(name='d')
 
-        c = Function(name='c', grid=grid, shape=(2, 3), dimensions=(d, z))
+        c = Function(name='c', grid=grid, shape=(2, 5), dimensions=(d, z))
         u = TimeFunction(name='u', grid=grid, space_order=3)
         u1 = TimeFunction(name='u1', grid=grid, space_order=3)
 

From beaa385248dc1267bd79dd7d7acb0ed745b295b7 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Thu, 24 Sep 2026 16:58:58 +0100
Subject: [PATCH 4/6] compiler: Improve SubDimension DDA

---
 devito/ir/support/basic.py | 27 +++++++-----
 devito/types/dimension.py  | 24 +++++++---
 tests/test_dimension.py    |  9 ++--
 tests/test_ir.py           |  8 ++--
 tests/test_subdomains.py   | 89 +++++++++++++++++++++++++++++---------
 5 files changed, 109 insertions(+), 48 deletions(-)

diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py
index 034dbddfa63..127c5e5967a 100644
--- a/devito/ir/support/basic.py
+++ b/devito/ir/support/basic.py
@@ -1610,21 +1610,20 @@ def disjoint_subdims(a0, a1):
     of the same Function.
 
     Compare symbolic accessed bounds, including shifts and stencil points.
-    Block intervals are promoted to their logical SubDimensions. Bounds and
-    thicknesses remain symbolic: MPI decomposition and runtime overrides can
-    change their values independently of the defaults.
+    Block intervals are promoted to their logical SubDimensions. Declared
+    thicknesses determine the global regions: explicit overrides are forbidden,
+    while MPI clips these regions to each rank. Parent bounds and access offsets
+    remain symbolic; only iteration bounds use the declared thicknesses.
 
-    For example, `xl = [m, m + L - 1]` and `xm = [m + L, M - R]` are
-    disjoint when they share the symbol `L`, whatever its runtime value.
-    Equal default thicknesses alone do not establish that relationship.
+    For example, a left slab of thickness 4 ends before a middle excluding 4
+    points, even when the two thickness symbols are distinct.
 
     Opposite left/right slabs are assumed to form a valid partition: their
     thicknesses satisfy `L + R <= N`. For translated stencil accesses, the
     interior must also accommodate their combined inward reach. For example,
     a pointwise left write and a right read at offset -4 require four interior
     points. Runtime space_order checks cover explicit middle SubDimensions,
-    not arbitrary left/right pairs; no concrete domain size or thickness is
-    used here.
+    not arbitrary left/right pairs; no concrete domain size is used here.
 
     Match data axes independently of the iteration nests. Return True if any
     axis proves separation, False otherwise. Accesses over the same interval
@@ -1644,6 +1643,8 @@ def disjoint_subdims(a0, a1):
                 it0 != it1):
             continue
 
+        thicknesses = {t: t.value for it in (it0, it1)
+                       for t in it.dim.thickness if t.value is not None}
         bounds = []
         for e, d, it in ((e0, d0, it0), (e1, d1, it1)):
             if not q_affine(e, d):
@@ -1658,8 +1659,8 @@ def disjoint_subdims(a0, a1):
                     M, m = it.symbolic_min, it.symbolic_max
                 else:
                     break
-                lower.append(v._subs(d, m))
-                upper.append(v._subs(d, M))
+                lower.append(v._subs(d, m.xreplace(thicknesses)))
+                upper.append(v._subs(d, M.xreplace(thicknesses)))
             else:
                 bounds.append((sympy.Min(*lower), sympy.Max(*upper)))
 
@@ -1682,11 +1683,13 @@ def disjoint_subdims(a0, a1):
                 gap = sympy.Dummy(nonnegative=True)
                 if e0.diff(d0) == e1.diff(d1) == 1:
                     M, m = (M0, m1) if it0.dim.is_left else (M1, m0)
-                    reach = (M - dl.symbolic_max - m + dr.symbolic_min).expand()
+                    reach = (M - dl.symbolic_max.xreplace(thicknesses) - m +
+                             dr.symbolic_min.xreplace(thicknesses)).expand()
                     if is_integer(reach):
                         gap += max(0, reach)
 
-                mapper[dlp.symbolic_max] = dlp.symbolic_min + dl.ltkn + dr.rtkn + gap - 1
+                mapper[dlp.symbolic_max] = (dlp.symbolic_min + dl.ltkn.value +
+                                            dr.rtkn.value + gap - 1)
 
             if (M0 - m1).subs(mapper).is_negative or \
                (M1 - m0).subs(mapper).is_negative:
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index c131cd39d06..570d1771211 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -554,7 +554,12 @@ def _arg_check(self, *args, **kwargs):
 # the user
 
 class Thickness(DataSymbol):
-    """A DataSymbol to represent a thickness of a SubDimension"""
+    """
+    A SubDimension thickness, fixed at construction and localized by MPI.
+
+    Explicit runtime overrides are not supported: dependence analysis uses the
+    declared thickness to determine the global region before MPI decomposition.
+    """
 
     __rkwargs__ = DataSymbol.__rkwargs__ + ('root', 'side', 'local', 'value')
 
@@ -589,11 +594,7 @@ def value(self):
         return self._value
 
     def _arg_values(self, grid=None, **kwargs):
-        # Allow override of thickness values to disable BCs
-        # However, arguments from the user are considered global
-        # So overriding the thickness to a nonzero value should not cause
-        # boundaries to exist between ranks where they did not before
-        rtkn = kwargs.get(self.name, self.value)
+        rtkn = self.value
         if grid is not None and grid.is_distributed(self.root):
             # Get local thickness
             if self.local:
@@ -612,6 +613,15 @@ def _arg_values(self, grid=None, **kwargs):
 
         return {self.name: tkn}
 
+    def _arg_check(self, args, *_args, **kwargs):
+        # This module depends on Dimension, so importing it above would cycle
+        from devito.mpi import mpi_raise  # noqa: PLC0415
+
+        error = (f"Cannot override SubDimension thickness `{self.name}`"
+                 if self.name in kwargs else None)
+        comm = args.comm if args.options['mpi'] else None
+        mpi_raise(error, InvalidArgument, comm=comm)
+
 
 class AbstractSubDimension(DerivedDimension):
 
@@ -840,7 +850,7 @@ def _arg_check(self, args, *_args, **kwargs):
             values = {**args,
                       d.min_name: kwargs.get(d.min_name, 0),
                       d.max_name: kwargs.get(d.max_name, kwargs.get(d.name, size - 1)),
-                      **{t.name: kwargs.get(t.name, t.value) for t in self.thickness}}
+                      **{t.name: t.value for t in self.thickness}}
         else:
             values = args
         size = int(subs_op_args(self.symbolic_size, values))
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index 96bec129353..d40d002d1ed 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -794,17 +794,18 @@ def test_expandingbox_like(self, opt):
         grid = Grid(shape=(8, 8))
         x, y = grid.dimensions
 
+        # Declare the widest box; runtime root bounds control the active box
+        xi = SubDimension.middle(name='xi', parent=x, thickness_left=0, thickness_right=0)
+        yi = SubDimension.middle(name='yi', parent=y, thickness_left=0, thickness_right=0)
+
         u = TimeFunction(name='u', grid=grid)
-        xi = SubDimension.middle(name='xi', parent=x, thickness_left=2, thickness_right=2)
-        yi = SubDimension.middle(name='yi', parent=y, thickness_left=2, thickness_right=2)
 
         eqn = Eq(u.forward, u + 1)
         eqn = eqn.subs({x: xi, y: yi})
 
         op = Operator(eqn, opt=opt)
 
-        op.apply(time=3, x_m=2, x_M=5, y_m=2, y_M=5,
-                 x_ltkn0=0, x_rtkn0=0, y_ltkn0=0, y_rtkn0=0)
+        op.apply(time=3, x_m=2, x_M=5, y_m=2, y_M=5)
 
         assert np.all(u.data[0, 2:-2, 2:-2] == 4.)
         assert np.all(u.data[1, 2:-2, 2:-2] == 3.)
diff --git a/tests/test_ir.py b/tests/test_ir.py
index 63a462b1a91..8629c9fa854 100644
--- a/tests/test_ir.py
+++ b/tests/test_ir.py
@@ -562,8 +562,6 @@ def test_subdimension_affine_bounds(self, shared_boundary, symbolic, slope, offs
         xm = SubDimension.middle('xm', x, 8, 40)
         if shared_boundary:
             xm = xm._rebuild(thickness=(xl.ltkn, xm.rtkn))
-        else:
-            expected = S.Infinity
         if symbolic:
             f = Array(name='f', dimensions=(x,))
         else:
@@ -1305,8 +1303,8 @@ def test_bundle_components(self):
 
     @pytest.mark.parametrize('symbolic', [False, True])
     @pytest.mark.parametrize('lower,upper,ndeps', [
-        (0, 31, (8, 32)), (0, 7, (0, 8)), (16, 31, (0, 16)),
-        (7, 8, (1, 2)), (15, 16, (1, 2))
+        (0, 31, (8, 24)), (0, 7, (0, 0)), (16, 31, (0, 16)),
+        (7, 8, (1, 1)), (15, 16, (1, 2))
     ])
     def test_stencil_contains_producer(self, symbolic, lower, upper, ndeps):
         grid = Grid(shape=(32,))
@@ -1317,6 +1315,8 @@ def test_stencil_contains_producer(self, symbolic, lower, upper, ndeps):
         f = Function(name='f', grid=grid)
         g = Function(name='g', grid=grid)
         if symbolic:
+            # Thickness 8 fixes the left boundary; the parent upper bound remains
+            # symbolic, so stencil offsets >= 8 can still touch the producer
             pi, ci = Interval(xm), Interval(xl)
         else:
             # Encode actual fixed iteration bounds, not runtime defaults
diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py
index 0fbb47c5666..c8f7e3686a3 100644
--- a/tests/test_subdomains.py
+++ b/tests/test_subdomains.py
@@ -1937,14 +1937,16 @@ def define(self, dimensions):
 class TestSubDomainArguments:
 
     @staticmethod
-    def _make_operator(left_shift=0, right_shift=0, grid=None, middle=False):
+    def _make_operator(left_shift=0, right_shift=0, grid=None, middle=False,
+                       thickness=(8, 8)):
         grid = grid or Grid(shape=(8, 32))
         y = grid.dimensions[-1]
 
-        yl = SubDimension.left('yl', y, 8)
-        yr = SubDimension.right('yr', y, 8)
+        left, right = thickness
+        yl = SubDimension.left('yl', y, left)
+        yr = SubDimension.right('yr', y, right)
         if middle:
-            yl = yr = SubDimension.middle('ym', y, 8, 8)
+            yl = yr = SubDimension.middle('ym', y, left, right)
 
         u = TimeFunction(name='u', grid=grid, space_order=8)
         v = TimeFunction(name='v', grid=grid, space_order=8)
@@ -1979,7 +1981,11 @@ def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
             dl, = [d for d in op.dimensions if d.is_Sub and not d.is_right]
             kwargs = {dl.ltkn.name: 32 - 8 - 8 - margin}
 
-        if middle and margin < 0:
+        if override == 'thickness':
+            with pytest.raises(InvalidArgument,
+                               match='Cannot override SubDimension thickness'):
+                op.arguments(time_M=0, **kwargs)
+        elif middle and margin < 0:
             with pytest.raises(InvalidArgument, match='at least 8 interior points'):
                 op.arguments(time_M=0, **kwargs)
         else:
@@ -2007,13 +2013,11 @@ def test_runtime_space_order(self, middle, space_order, margin):
 
     @pytest.mark.parametrize('side', ['left', 'right'])
     def test_empty_slab(self, side):
-        op, _ = self._make_operator(right_shift=-4)
-
-        d, = [d for d in op.dimensions if d.is_Sub and getattr(d, f'is_{side}')]
-        thickness = d.ltkn if side == 'left' else d.rtkn
+        thickness = (0, 8) if side == 'left' else (8, 0)
+        op, _ = self._make_operator(right_shift=-4, thickness=thickness)
 
         # The active slab fills the local domain; the opposite slab is absent
-        op.arguments(time_M=0, y_M=7, **{thickness.name: 0})
+        op.arguments(time_M=0, y_M=7)
 
     def test_before_autotuning(self):
         op, _ = self._make_operator(right_shift=-4, middle=True)
@@ -2032,19 +2036,15 @@ def test_before_autotuning(self):
     def test_distributed_middle(self, left, mode):
         grid = Grid(shape=(16, 32), topology=(1, 2))
 
-        op, _ = self._make_operator(grid=grid, middle=True)
-
-        d, = [d for d in op.dimensions if d.is_Sub]
+        op, _ = self._make_operator(grid=grid, middle=True, thickness=(left, 0))
 
         # Rank 0 has one point or an empty middle (possibly with inverted bounds).
         # Only the global size determines whether the middle is large enough
-        kwargs = {d.ltkn.name: left, d.rtkn.name: 0}
-
         if left == 25:
             with pytest.raises(InvalidArgument, match='at least 8 interior points'):
-                op.arguments(time_M=0, **kwargs)
+                op.arguments(time_M=0)
         else:
-            op.arguments(time_M=0, **kwargs)
+            op.arguments(time_M=0)
 
     @pytest.mark.parallel(mode=[(2, 'basic')])
     def test_collective_rejection(self, mode):
@@ -2052,12 +2052,11 @@ def test_collective_rejection(self, mode):
 
         op, _ = self._make_operator(right_shift=-4, grid=grid, middle=True)
 
-        dl, = [d for d in op.dimensions if d.is_Sub and d.is_middle]
-        left = 24 if grid.distributor.myrank == 0 else 0
+        upper = 22 if grid.distributor.myrank == 0 else 31
 
-        # Only rank 0 has an insufficient middle; its peer has 24 interior points
+        # Only rank 0 has an insufficient middle; its peer has 16 interior points
         with pytest.raises(InvalidArgument, match='interior points'):
-            op.arguments(time_M=0, **{dl.ltkn.name: left})
+            op.arguments(time_M=0, y_M=upper)
 
     def test_function_on_subdomain(self):
         class Interior(SubDomain):
@@ -2086,4 +2085,52 @@ def define(self, dimensions):
         assert concrete.ltkn.name in args
 
         with pytest.raises(InvalidArgument, match='at least 8 interior points'):
+            op.arguments(y_M=22)
+
+        with pytest.raises(InvalidArgument,
+                           match='Cannot override SubDimension thickness'):
             op.arguments(**{concrete.ltkn.name: 17})
+
+    def test_thickness_overrides(self):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        xi = SubDimension.left('xi', x, 4)
+
+        f = Function(name='f', grid=grid, space_order=0)
+
+        eq = Eq(f[xi], 1)
+
+        op = Operator(eq, name='thickness_overrides')
+
+        d, = [d for d in op.dimensions if d.is_Sub]
+
+        with pytest.raises(InvalidArgument,
+                           match='Cannot override SubDimension thickness'):
+            op.apply(**{d.ltkn.name: 2})
+
+        assert np.all(f.data == 0)
+
+    @pytest.mark.parallel(mode=[(2, 'basic')])
+    def test_collective_thickness_rejection(self, mode):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        xi = SubDimension.middle('xi', x, 4, 4)
+
+        f = Function(name='f', grid=grid, space_order=0)
+
+        eq = Eq(f[xi], 1)
+
+        op = Operator(eq, name='collective_thickness_rejection')
+
+        d, = [d for d in op.dimensions if d.is_Sub]
+        t = d.ltkn
+
+        # MPI clips the declared thicknesses without any explicit overrides
+        args = op.arguments()
+        assert set(grid.distributor.comm.allgather(args[t.name])) == {0, 4}
+
+        # Only rank 0 supplies an override; all ranks must reject it
+        kwargs = {t.name: 4} if grid.distributor.myrank == 0 else {}
+        with pytest.raises(InvalidArgument,
+                           match='Cannot override SubDimension thickness'):
+            op.arguments(**kwargs)

From 1c1616f9e221600d7c4f6856cf606e072f3b90ce Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Fri, 25 Sep 2026 08:54:12 +0100
Subject: [PATCH 5/6] compiler: Simplify and strengthten SubDimensions DDA

---
 devito/exceptions.py        |  16 ++++
 devito/ir/support/basic.py  |  36 +++-----
 devito/mpi/distributed.py   |  16 ----
 devito/operator/operator.py |   4 +-
 devito/types/dimension.py   | 147 +++++++++++++++++------------
 tests/test_dimension.py     |  53 ++++++++---
 tests/test_ir.py            |  30 +++++-
 tests/test_mpi.py           |   6 +-
 tests/test_pickle.py        |  18 +++-
 tests/test_subdomains.py    | 178 ++++++++++++++++++++++++++++--------
 10 files changed, 349 insertions(+), 155 deletions(-)

diff --git a/devito/exceptions.py b/devito/exceptions.py
index b15c5fd32fb..281419d39e7 100644
--- a/devito/exceptions.py
+++ b/devito/exceptions.py
@@ -54,3 +54,19 @@ class ExecutionError(DevitoError):
         * Device shared memory or registers (e.g., too many threads per block);
     * etc.
     """
+
+
+def mpi_raise(error, exception=ValueError, comm=None):
+    """
+    Raise `exception` with the first non-None error message in rank order.
+
+    All ranks in `comm` must call this routine, including those with no local
+    error (`error=None`). This prevents a rank-local exception from stranding
+    peers in subsequent MPI calls. With no communicator or `MPI.COMM_NULL`,
+    only the local error is checked.
+    """
+    # A null MPI communicator is false, like None
+    if comm:
+        error = next((i for i in comm.allgather(error) if i is not None), None)
+    if error is not None:
+        raise exception(error)
diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py
index 127c5e5967a..b51949d2862 100644
--- a/devito/ir/support/basic.py
+++ b/devito/ir/support/basic.py
@@ -1157,7 +1157,7 @@ def reads_smart_gen(self, f):
 
         The reason SubDimensions must be treated specially -- with a full set
         of TimedAccess objects getting generated -- is to handle the special
-        case of slabs thinner than the stencil’s reach. For example, consider
+        case of SubDimensions thinner than the stencil’s reach. For example, consider
         the following scenario:
 
         * A SubDimension with just two points, 10 and 11;
@@ -1615,15 +1615,14 @@ def disjoint_subdims(a0, a1):
     while MPI clips these regions to each rank. Parent bounds and access offsets
     remain symbolic; only iteration bounds use the declared thicknesses.
 
-    For example, a left slab of thickness 4 ends before a middle excluding 4
-    points, even when the two thickness symbols are distinct.
+    For example, a left SubDimension of thickness 4 ends before a middle
+    SubDimension excluding 4 points, even when the two thickness symbols are distinct.
 
-    Opposite left/right slabs are assumed to form a valid partition: their
-    thicknesses satisfy `L + R <= N`. For translated stencil accesses, the
-    interior must also accommodate their combined inward reach. For example,
-    a pointwise left write and a right read at offset -4 require four interior
-    points. Runtime space_order checks cover explicit middle SubDimensions,
-    not arbitrary left/right pairs; no concrete domain size is used here.
+    Left/right SubDimensions of the same parent with `overlap=False` satisfy
+    `L + R <= N`, checked against the full global parent extent at runtime.
+    Their regions may be adjacent; stencil accesses extending from one into
+    the other may still induce dependences. If either SubDimension allows
+    overlap, no minimum separation is assumed.
 
     Match data axes independently of the iteration nests. Return True if any
     axis proves separation, False otherwise. Accesses over the same interval
@@ -1671,23 +1670,10 @@ def disjoint_subdims(a0, a1):
             dl, dr = (it0.dim, it1.dim) if it0.dim.is_left else (it1.dim, it0.dim)
             dlp, drp = dl.parent, dr.parent
 
-            if dl.is_left and dr.is_right and dlp is drp:
-                # A valid partition satisfies L + R <= N, where N is the parent
-                # extent; an explicit middle SubDomain checks this at construction.
-                # Further, for stencils, we require that:
-                # `N - L - R >= the combined inward reach`
-                # so accesses from opposite slabs cannot meet. Explicit middle
-                # SubDimensions check for at least space_order interior points
-                # at *op.apply time*, accounting for runtime overrides.  Without
-                # an explicit middle, the gap assumption is unchecked
+            if dl.is_left and dr.is_right and dlp is drp and \
+               not (dl.overlap or dr.overlap):
+                # The runtime partition check guarantees a nonnegative gap only
                 gap = sympy.Dummy(nonnegative=True)
-                if e0.diff(d0) == e1.diff(d1) == 1:
-                    M, m = (M0, m1) if it0.dim.is_left else (M1, m0)
-                    reach = (M - dl.symbolic_max.xreplace(thicknesses) - m +
-                             dr.symbolic_min.xreplace(thicknesses)).expand()
-                    if is_integer(reach):
-                        gap += max(0, reach)
-
                 mapper[dlp.symbolic_max] = (dlp.symbolic_min + dl.ltkn.value +
                                             dr.rtkn.value + gap - 1)
 
diff --git a/devito/mpi/distributed.py b/devito/mpi/distributed.py
index 2c3c4ad3641..35596c39d64 100644
--- a/devito/mpi/distributed.py
+++ b/devito/mpi/distributed.py
@@ -70,7 +70,6 @@ def __getattr__(self, name):
     'SubDistributor',
     'devito_mpi_finalize',
     'devito_mpi_init',
-    'mpi_raise',
 ]
 
 
@@ -102,21 +101,6 @@ def devito_mpi_finalize():
         MPI.Finalize()
 
 
-def mpi_raise(error, exception=ValueError, comm=None):
-    """
-    Raise `exception` with the first non-None error message in rank order.
-
-    All ranks in `comm` must call this routine, including those with no local
-    error (`error=None`). This prevents a rank-local exception from stranding
-    peers in subsequent MPI calls. With no communicator or `MPI.COMM_NULL`,
-    only the local error is checked.
-    """
-    if comm is not None and comm is not MPI.COMM_NULL:
-        error = next((i for i in comm.allgather(error) if i is not None), None)
-    if error is not None:
-        raise exception(error)
-
-
 class AbstractDistributor(ABC):
 
     """
diff --git a/devito/operator/operator.py b/devito/operator/operator.py
index 875eb584a7f..1b4e9957ea1 100644
--- a/devito/operator/operator.py
+++ b/devito/operator/operator.py
@@ -43,7 +43,7 @@
     split, timed_pass, timed_region
 )
 from devito.types import Buffer, Evaluable, device_layer, disk_layer, host_layer
-from devito.types.dimension import Thickness
+from devito.types.dimension import SubDimension, Thickness
 from devito.warnings import warn
 
 __all__ = ['Operator']
@@ -718,6 +718,8 @@ def _prepare_arguments(self, autotune=None, estimate_memory=False, **kwargs):
             if d.is_Derived:
                 d._arg_check(args, **kwargs)
 
+        SubDimension._arg_check_thickness(self.dimensions, args)
+
         # Turn arguments into a format suitable for the generated code
         # E.g., instead of NumPy arrays for Functions, the generated code expects
         # pointers to ctypes.Struct
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index 570d1771211..71b484634d3 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -9,9 +9,9 @@
 
 from devito.data import LEFT, RIGHT
 from devito.deprecations import deprecations
-from devito.exceptions import InvalidArgument
+from devito.exceptions import InvalidArgument, mpi_raise
 from devito.logger import debug
-from devito.tools import Pickable, is_integer, is_number, memoized_meth
+from devito.tools import Pickable, as_mapper, is_integer, is_number, memoized_meth
 from devito.types.args import ArgProvider
 from devito.types.basic import DataSymbol, Scalar, Symbol
 from devito.types.constant import Constant
@@ -561,13 +561,14 @@ class Thickness(DataSymbol):
     declared thickness to determine the global region before MPI decomposition.
     """
 
-    __rkwargs__ = DataSymbol.__rkwargs__ + ('root', 'side', 'local', 'value')
+    __rkwargs__ = DataSymbol.__rkwargs__ + ('root', 'side', 'local', 'value', 'overlap')
 
-    def __new__(cls, *args, root=None, side=None, local=False, **kwargs):
+    def __new__(cls, *args, root=None, side=None, local=False, overlap=False, **kwargs):
         newobj = super().__new__(cls, *args, **kwargs)
         newobj._root = root
         newobj._side = side
         newobj._local = local
+        newobj._overlap = overlap
 
         return newobj
 
@@ -589,6 +590,10 @@ def side(self):
     def local(self):
         return self._local
 
+    @property
+    def overlap(self):
+        return self._overlap
+
     @property
     def value(self):
         return self._value
@@ -614,11 +619,9 @@ def _arg_values(self, grid=None, **kwargs):
         return {self.name: tkn}
 
     def _arg_check(self, args, *_args, **kwargs):
-        # This module depends on Dimension, so importing it above would cycle
-        from devito.mpi import mpi_raise  # noqa: PLC0415
-
         error = (f"Cannot override SubDimension thickness `{self.name}`"
                  if self.name in kwargs else None)
+
         comm = args.comm if args.options['mpi'] else None
         mpi_raise(error, InvalidArgument, comm=comm)
 
@@ -632,20 +635,31 @@ class AbstractSubDimension(DerivedDimension):
     Notes
     -----
     This is just the abstract base class for various types of SubDimensions.
+
+    SubDimensions cannot be nested.
     """
 
     is_AbstractSub = True
 
     __rargs__ = DerivedDimension.__rargs__ + ('thickness',)
-    __rkwargs__ = ()
+    __rkwargs__ = ('overlap',)
 
     _thickness_type = Thickness
 
-    def __init_finalize__(self, name, parent, thickness, **kwargs):
+    def __init_finalize__(self, name, parent, thickness, overlap=False, **kwargs):
         super().__init_finalize__(name, parent)
+        if parent.is_AbstractSub:
+            raise ValueError("Nested SubDimensions are not supported")
+
+        self._overlap = overlap
+
         thickness = thickness or (None, None)
         if any(isinstance(tkn, self._thickness_type) for tkn in thickness):
-            self._thickness = SubDimensionThickness(*thickness)
+            # Preserve identity when unchanged: expressions may share these symbols
+            self._thickness = SubDimensionThickness(*[
+                t if t.overlap == overlap else t._rebuild(overlap=overlap)
+                for t in thickness
+            ])
         else:
             self._thickness = self._symbolic_thickness(thickness=thickness)
 
@@ -657,7 +671,8 @@ def _interval(self):
 
     @memoized_meth
     def _symbolic_thickness(self, **kwargs):
-        kwargs = {'dtype': np.int32, 'is_const': True, 'nonnegative': True}
+        kwargs = {'dtype': np.int32, 'is_const': True, 'nonnegative': True,
+                  'overlap': self.overlap}
 
         names = [f"{self.parent.name}_{s}tkn" for s in ('l', 'r')]
         return SubDimensionThickness(*[Thickness(name=n, **kwargs) for n in names])
@@ -691,6 +706,10 @@ def rtkn(self):
         # Shortcut for the right thickness symbol
         return self.thickness.right
 
+    @property
+    def overlap(self):
+        return self._overlap
+
     def __hash__(self):
         return id(self)
 
@@ -718,6 +737,11 @@ class SubDimension(AbstractSubDimension):
     local : bool
         True if, in case of domain decomposition, the SubDimension is
         guaranteed not to span more than one domain, False otherwise.
+    overlap : bool, optional, default=False
+        Allow a left SubDimension to overlap a right SubDimension of the same
+        parent. When both have this flag unset, their declared thicknesses must sum
+        to at most the full global parent extent. This does not guarantee
+        separation of stencil accesses extending beyond the SubDimensions themselves.
 
     Examples
     --------
@@ -741,6 +765,9 @@ class SubDimension(AbstractSubDimension):
     local (i.e., non-distributed) Dimensions, as they are assumed to fit entirely
     within a single domain. This is the most typical use case (e.g., to set up
     boundary conditions). To drop this assumption, pass ``local=False``.
+
+    Explicit runtime overrides of the root Dimension's bounds are not supported.
+    MPI may still localize the global regions automatically.
     """
 
     is_Sub = True
@@ -749,27 +776,31 @@ class SubDimension(AbstractSubDimension):
 
     _thickness_type = Thickness
 
-    def __init_finalize__(self, name, parent, thickness, local,
-                          **kwargs):
+    def __init_finalize__(self, name, parent, thickness, local, **kwargs):
         self._local = local
-        super().__init_finalize__(name, parent, thickness)
+        super().__init_finalize__(name, parent, thickness, **kwargs)
 
     @classmethod
-    def left(cls, name, parent, thickness, local=True):
-        return cls(name, parent, thickness=(thickness, None), local=local)
+    def left(cls, name, parent, thickness, local=True, overlap=False):
+        return cls(name, parent, thickness=(thickness, None), local=local,
+                   overlap=overlap)
 
     @classmethod
-    def right(cls, name, parent, thickness, local=True):
-        return cls(name, parent, thickness=(None, thickness), local=local)
+    def right(cls, name, parent, thickness, local=True, overlap=False):
+        return cls(name, parent, thickness=(None, thickness), local=local,
+                   overlap=overlap)
 
     @classmethod
-    def middle(cls, name, parent, thickness_left, thickness_right, local=False):
-        return cls(name, parent, thickness=(thickness_left, thickness_right), local=local)
+    def middle(cls, name, parent, thickness_left, thickness_right, local=False,
+               overlap=False):
+        return cls(name, parent, thickness=(thickness_left, thickness_right), local=local,
+                   overlap=overlap)
 
     @memoized_meth
     def _symbolic_thickness(self, thickness=None):
         kwargs = {'dtype': np.int32, 'is_const': True, 'nonnegative': True,
-                  'root': self.root, 'local': self.local}
+                  'root': self.root, 'local': self.local,
+                  'overlap': self.overlap}
 
         names = [f"{self.parent.name}_{s}tkn" for s in ('l', 'r')]
         sides = [LEFT, RIGHT]
@@ -829,41 +860,43 @@ def _arg_values(self, interval, grid=None, **kwargs):
         return {}
 
     def _arg_check(self, args, *_args, **kwargs):
-        # These modules depend on Dimension, so importing them above would cycle
-        from devito.mpi import mpi_raise  # noqa: PLC0415
-        from devito.symbolics import subs_op_args  # noqa: PLC0415
-
-        if not self.is_middle:
-            return
-
-        # Function._arg_check visits original axes (e.g. x_ltkn), whereas `args`
-        # contains the concretized thicknesses (x_ltkn0, ...). The Operator checks
-        # the matching concrete SubDimensions separately in its dimension loop
-        if self not in args.op.dimensions:
-            return
-
         d = self.root
-        if args.grid is not None and args.grid.is_distributed(d):
-            # Check the global runtime region: its MPI-local slices may be empty
-            # or smaller than space_order even for a non-degenerate global interior
-            size = args.grid.size_map[d].glb
-            values = {**args,
-                      d.min_name: kwargs.get(d.min_name, 0),
-                      d.max_name: kwargs.get(d.max_name, kwargs.get(d.name, size - 1)),
-                      **{t.name: t.value for t in self.thickness}}
-        else:
-            values = args
-        size = int(subs_op_args(self.symbolic_size, values))
+        names = [k for k in (d.min_name, d.max_name, d.name) if k in kwargs]
+        error = (f"Cannot override bounds {names} of Dimension `{d}` used by "
+                 f"SubDimension `{self}`" if names else None)
 
-        # Runtime overrides do not change the compiled stencil order
-        items = [f.space_order for f in args.op.input if f.is_DiscreteFunction]
-        space_order = max(items, default=0)
+        comm = args.comm if args.options['mpi'] else None
+        mpi_raise(error, InvalidArgument, comm=comm)
 
-        if size < space_order:
-            error = (f"Expected at least {space_order} interior points along "
-                     f"`{self.parent}` (space_order), but runtime arguments leave {size}")
-        else:
-            error = None
+    @classmethod
+    def _arg_check_thickness(cls, dimensions, args):
+        """
+        Check that non-overlapping left/right SubDimensions fit the global domain.
+        """
+        subdims = [d for d in dimensions if isinstance(d, cls) and not d.overlap]
+
+        grid = args.grid
+        error = None
+
+        for parent, dims in as_mapper(subdims, lambda d: d.parent).items():
+            # Check left/right SubDimensions, not a middle's excluded thicknesses
+            left = max((d for d in dims if d.is_left),
+                       key=lambda d: d.ltkn.value, default=None)
+            right = max((d for d in dims if d.is_right),
+                        key=lambda d: d.rtkn.value, default=None)
+            if left is None or right is None:
+                continue
+
+            root = parent.root
+
+            # Use the full global extent, not a rank's local extent
+            size = grid.size_map[root].glb if grid is not None else args[root.size_name]
+            thickness = left.ltkn.value + right.rtkn.value
+            if thickness > size:
+                error = (f"SubDimensions `{left}` and `{right}` have combined "
+                         f"thickness {thickness} along `{parent}`, exceeding "
+                         f"the runtime extent {size}")
+                break
 
         comm = args.comm if args.options['mpi'] else None
         mpi_raise(error, InvalidArgument, comm=comm)
@@ -877,12 +910,14 @@ class MultiSubDimension(AbstractSubDimension):
 
     is_MultiSub = True
 
-    __rkwargs__ = ('functions', 'bounds_indices', 'implicit_dimension')
+    __rkwargs__ = AbstractSubDimension.__rkwargs__ + (
+        'functions', 'bounds_indices', 'implicit_dimension'
+    )
 
     def __init_finalize__(self, name, parent, thickness, functions=None,
-                          bounds_indices=None, implicit_dimension=None):
+                          bounds_indices=None, implicit_dimension=None, **kwargs):
 
-        super().__init_finalize__(name, parent, thickness)
+        super().__init_finalize__(name, parent, thickness, **kwargs)
         self.functions = functions
         self.bounds_indices = bounds_indices
         self.implicit_dimension = implicit_dimension
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index d40d002d1ed..6bfba0b79dd 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -12,6 +12,7 @@
     SparseTimeFunction, SubDimension, SubDomain, TimeFunction, configuration, dimensions,
     floor, norm, sin, sum, switchconfig
 )
+from devito.exceptions import InvalidArgument
 from devito.ir import SymbolRegistry
 from devito.ir.equations.algorithms import concretize_subdims
 from devito.ir.iet import (
@@ -21,7 +22,7 @@
 from devito.symbolics import INT, IntDiv, indexify, retrieve_functions
 from devito.types import Array, StencilDimension, Symbol
 from devito.types.basic import Scalar
-from devito.types.dimension import AffineIndexAccessFunction, Thickness
+from devito.types.dimension import AffineIndexAccessFunction, MultiSubDimension, Thickness
 from devito.types.misc import Temp
 
 
@@ -787,16 +788,15 @@ def test_arrays_defined_over_subdims(self):
         op()
 
     @pytest.mark.parametrize('opt', opts_tiling)
-    def test_expandingbox_like(self, opt):
+    def test_box_bounds(self, opt):
         """
-        Make sure SubDimensions aren't an obstacle to expanding boxes.
+        SubDimension boxes use declared thicknesses, not runtime root bounds.
         """
         grid = Grid(shape=(8, 8))
         x, y = grid.dimensions
 
-        # Declare the widest box; runtime root bounds control the active box
-        xi = SubDimension.middle(name='xi', parent=x, thickness_left=0, thickness_right=0)
-        yi = SubDimension.middle(name='yi', parent=y, thickness_left=0, thickness_right=0)
+        xi = SubDimension.middle(name='xi', parent=x, thickness_left=2, thickness_right=2)
+        yi = SubDimension.middle(name='yi', parent=y, thickness_left=2, thickness_right=2)
 
         u = TimeFunction(name='u', grid=grid)
 
@@ -805,7 +805,10 @@ def test_expandingbox_like(self, opt):
 
         op = Operator(eqn, opt=opt)
 
-        op.apply(time=3, x_m=2, x_M=5, y_m=2, y_M=5)
+        with pytest.raises(InvalidArgument, match='Cannot override bounds'):
+            op.apply(time=3, x_m=2, x_M=5, y_m=2, y_M=5)
+
+        op.apply(time=3)
 
         assert np.all(u.data[0, 2:-2, 2:-2] == 4.)
         assert np.all(u.data[1, 2:-2, 2:-2] == 3.)
@@ -826,6 +829,18 @@ def test_standalone_thickness(self):
         op(x_m=0)
         assert np.all(f.data == np.array([0, 1, 0, 0, 0]))
 
+    def test_no_nesting(self):
+        x = Dimension('x')
+        xi = SubDimension.middle('xi', x, 1, 1)
+        xm = MultiSubDimension('xm', x, None)
+
+        for parent in (xi, xm):
+            with pytest.raises(ValueError, match='Nested SubDimensions'):
+                SubDimension.left('xl', parent, 1)
+
+            with pytest.raises(ValueError, match='Nested SubDimensions'):
+                MultiSubDimension('xs', parent, None)
+
 
 class TestConditionalDimension:
 
@@ -2242,20 +2257,36 @@ class TestConcretization:
     during compilation.
     """
 
-    def test_correct_thicknesses(self):
+    @pytest.mark.parametrize('overlap', [False, True])
+    def test_correct_thicknesses(self, overlap):
         """
         Check that thicknesses aren't created where they shouldn't be.
         """
         x = Dimension('x')
-        ix0 = SubDimension.left('x', x, 2)
-        ix1 = SubDimension.right('x', x, 2)
-        ix2 = SubDimension.middle('x', x, 2, 2)
+        ix0 = SubDimension.left('x', x, 2, overlap=overlap)
+        ix1 = SubDimension.right('x', x, 2, overlap=overlap)
+        ix2 = SubDimension.middle('x', x, 2, 2, overlap=overlap)
 
         rebuilt = concretize_subdims([ix0, ix1, ix2], sregistry=SymbolRegistry())
 
         assert rebuilt[0].is_left
         assert rebuilt[1].is_right
         assert rebuilt[2].is_middle
+        assert all(d.overlap is overlap for d in rebuilt)
+        assert all(t.overlap is overlap for d in rebuilt for t in d.thickness)
+
+        for d in rebuilt:
+            changed = d._rebuild(overlap=not overlap)
+            assert changed.overlap is not overlap
+            assert all(t.overlap is not overlap for t in changed.thickness)
+
+    def test_shared_thickness(self):
+        x = Dimension('x')
+        xl = SubDimension.left('xl', x, 4)
+
+        d, t = concretize_subdims([xl, xl.ltkn], sregistry=SymbolRegistry())
+
+        assert d.ltkn is t
 
     def test_condition_concretization(self):
         """
diff --git a/tests/test_ir.py b/tests/test_ir.py
index 8629c9fa854..a730b445f89 100644
--- a/tests/test_ir.py
+++ b/tests/test_ir.py
@@ -500,7 +500,7 @@ def test_subdimension_stencil_distance(self, shared_boundary, offset, independen
         a = TimedAccess(f[xl], 'W', 0, IterationSpace([Interval(xl)]))
         b = TimedAccess(f[xr + offset + h], 'R', 1,
                         IterationSpace([interval]))
-        independent = independent if shared_boundary else True
+        independent = independent if shared_boundary else offset >= 0
         assert (S.ImaginaryUnit in a.distance(b)) is independent
         assert (S.ImaginaryUnit in b.distance(a)) is independent
 
@@ -509,7 +509,7 @@ def test_subdimension_stencil_distance(self, shared_boundary, offset, independen
     @pytest.mark.parametrize('side,thickness,shift,expected', [
         ('right', 20, 0, (S.ImaginaryUnit,)),
         ('right', 24, 0, (S.ImaginaryUnit,)),
-        ('right', 20, -8, (S.ImaginaryUnit,)),
+        ('right', 20, -8, (S.Infinity, S.Infinity)),
         ('middle', 8, 0, (S.ImaginaryUnit,)),
         ('middle', 8, -1, (S.Infinity, S.Infinity))
     ])
@@ -614,11 +614,31 @@ def test_opposite_subdimension_bounds(self, slope, offset, expected):
         b = TimedAccess(f[slope*(xr + offset)], 'R', 1,
                         IterationSpace([Interval(xr)]))
 
-        # Translated accesses assume the interior accommodates their inward reach
-        if slope == 1:
-            expected = S.ImaginaryUnit
+        # Left/right SubDimensions may be adjacent, so inward reads can overlap
         assert a.distance(b) == b.distance(a) == (expected,)
 
+    @pytest.mark.parametrize('side', ['left', 'right'])
+    def test_overlap_subdimension_bounds(self, side):
+        grid = Grid(shape=(32, 32))
+        x, y = grid.dimensions
+        xl = SubDimension.left('xl', x, 8, overlap=side == 'left')
+        xr = SubDimension.right('xr', x, 8, overlap=side == 'right')
+        yl = SubDimension.left('yl', y, 4)
+        ym = SubDimension.middle('ym', y, 4, 0)
+
+        f = Function(name='f', grid=grid)
+
+        a = TimedAccess(f[xl, yl], 'W', 0,
+                        IterationSpace([Interval(xl), Interval(yl)]))
+        b = TimedAccess(f[xr, yl], 'R', 1,
+                        IterationSpace([Interval(xr), Interval(yl)]))
+        c = TimedAccess(f[xr, ym], 'R', 1,
+                        IterationSpace([Interval(xr), Interval(ym)]))
+
+        assert a.distance(b) == b.distance(a) == (S.Infinity, 0)
+        # Permitting overlap along x does not manufacture a dependence along y
+        assert a.distance(c) == c.distance(a) == (S.ImaginaryUnit,)
+
 
 class TestSpace:
 
diff --git a/tests/test_mpi.py b/tests/test_mpi.py
index ac3408d0c41..24e46074fb1 100644
--- a/tests/test_mpi.py
+++ b/tests/test_mpi.py
@@ -13,14 +13,14 @@
 )
 from devito.arch.compiler import OneapiCompiler
 from devito.data import LEFT, RIGHT
-from devito.exceptions import InvalidArgument
+from devito.exceptions import InvalidArgument, mpi_raise
 from devito.ir import Cluster, Interval, IterationSpace
 from devito.ir.clusters.algorithms import check_halo_writes
 from devito.ir.iet import (
     Call, Conditional, FindNodes, FindSymbols, Iteration, retrieve_iteration_tree
 )
 from devito.ir.support.space import Backward, Forward
-from devito.mpi import MPI, mpi_raise
+from devito.mpi import MPI
 from devito.mpi.distributed import CustomTopology
 from devito.mpi.routines import ComputeCall, HaloUpdateCall, HaloUpdateList, MPICall
 from devito.tools import Bunch
@@ -613,7 +613,7 @@ def test_local_indices(self, shape, expected, mode):
         )
 
     @pytest.mark.parallel(mode=4)
-    @pytest.mark.parametrize('shape', [(1,), (2, 3), (4, 5, 6)])
+    @pytest.mark.parametrize('shape', [(2,), (2, 3), (4, 5, 6)])
     def test_mpi4py_nodevmpi(self, shape, mode):
 
         with switchconfig(mpi=False):
diff --git a/tests/test_pickle.py b/tests/test_pickle.py
index 954d39cc80b..fca861c1444 100644
--- a/tests/test_pickle.py
+++ b/tests/test_pickle.py
@@ -30,6 +30,7 @@
 from devito.types import Symbol as dSymbol
 from devito.types import TempFunction, ThreadID, Timer
 from devito.types.basic import AbstractSymbol, BoundSymbol
+from devito.types.dimension import MultiSubDimension
 from examples.seismic import (
     AcquisitionGeometry, Receiver, RickerSource, TimeAxis, demo_model
 )
@@ -332,7 +333,7 @@ def test_array(self, pickle):
         assert new_pa.array.name == 'a'
 
     def test_sub_dimension(self, pickle):
-        di = SubDimension.middle('di', Dimension(name='d'), 1, 1)
+        di = SubDimension.middle('di', Dimension(name='d'), 1, 1, overlap=True)
 
         pkl_di = pickle.dumps(di)
         new_di = pickle.loads(pkl_di)
@@ -342,6 +343,21 @@ def test_sub_dimension(self, pickle):
         assert di.parent.name == new_di.parent.name
         assert di._thickness == new_di._thickness
         assert di._interval == new_di._interval
+        assert new_di.overlap
+        assert all(t.overlap for t in new_di.thickness)
+
+    def test_multi_sub_dimension(self, pickle):
+        di = MultiSubDimension('di', Dimension(name='d'), None, overlap=True)
+
+        new_di = pickle.loads(pickle.dumps(di))
+
+        assert new_di.overlap
+        assert all(t.overlap for t in new_di.thickness)
+
+        rebuilt = new_di._rebuild(overlap=False)
+
+        assert not rebuilt.overlap
+        assert all(not t.overlap for t in rebuilt.thickness)
 
     def test_conditional_dimension(self, pickle):
         d = Dimension(name='d')
diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py
index c8f7e3686a3..8f02caf841f 100644
--- a/tests/test_subdomains.py
+++ b/tests/test_subdomains.py
@@ -152,7 +152,8 @@ class sd0(SubDomain):
 
             def define(self, dimensions):
                 x, y = dimensions
-                return {x: ('middle', 2, 2), y: ('right', 10)}
+                return {x: ('middle', 2, 2),
+                        y: SubDimension.right('iy', y, 10, overlap=True)}
 
         class sd1(SubDomain):
             name = 'sd1'
@@ -1965,7 +1966,7 @@ def _make_operator(left_shift=0, right_shift=0, grid=None, middle=False,
         (0, 0), (0, -4), (2, -4)
     ])
     @pytest.mark.parametrize('margin', [-1, 0, 1])
-    @pytest.mark.parametrize('override', ['function', 'bounds', 'thickness'])
+    @pytest.mark.parametrize('override', ['function', 'thickness'])
     def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
         op, (u, v) = self._make_operator(left_shift, right_shift, middle=middle)
 
@@ -1975,8 +1976,6 @@ def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
 
             kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
                                            space_order=8) for f in (u, v)}
-        elif override == 'bounds':
-            kwargs = {'y_m': 3, 'y_M': 3 + size - 1}
         else:
             dl, = [d for d in op.dimensions if d.is_Sub and not d.is_right]
             kwargs = {dl.ltkn.name: 32 - 8 - 8 - margin}
@@ -1985,9 +1984,6 @@ def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
             with pytest.raises(InvalidArgument,
                                match='Cannot override SubDimension thickness'):
                 op.arguments(time_M=0, **kwargs)
-        elif middle and margin < 0:
-            with pytest.raises(InvalidArgument, match='at least 8 interior points'):
-                op.arguments(time_M=0, **kwargs)
         else:
             op.arguments(time_M=0, **kwargs)
 
@@ -1995,7 +1991,7 @@ def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
     @pytest.mark.parametrize('space_order', [4, 8, 12])
     @pytest.mark.parametrize('margin', [-1, 0, 1])
     def test_runtime_space_order(self, middle, space_order, margin):
-        """Override metadata does not change the compiled stencil order."""
+        """Neither compiled nor override space_order constrains region widths."""
         op, fields = self._make_operator(middle=middle)
 
         required = 8
@@ -2004,30 +2000,34 @@ def test_runtime_space_order(self, middle, space_order, margin):
         kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
                                        space_order=space_order) for f in fields}
 
-        if middle and margin < 0:
-            with pytest.raises(InvalidArgument,
-                               match=f'at least {required} interior points'):
-                op.arguments(time_M=0, **kwargs)
-        else:
-            op.arguments(time_M=0, **kwargs)
+        op.arguments(time_M=0, **kwargs)
 
     @pytest.mark.parametrize('side', ['left', 'right'])
-    def test_empty_slab(self, side):
+    def test_empty_subdimension(self, side):
+        grid = Grid(shape=(8, 8))
+
         thickness = (0, 8) if side == 'left' else (8, 0)
-        op, _ = self._make_operator(right_shift=-4, thickness=thickness)
+        op, _ = self._make_operator(right_shift=-4, thickness=thickness, grid=grid)
 
-        # The active slab fills the local domain; the opposite slab is absent
-        op.arguments(time_M=0, y_M=7)
+        # One left/right SubDimension fills the domain; the other is empty
+        op.arguments(time_M=0)
 
     def test_before_autotuning(self):
-        op, _ = self._make_operator(right_shift=-4, middle=True)
+        grid = Grid(shape=(8, 15))
+
+        op, fields = self._make_operator(right_shift=-4, grid=grid)
 
-        with pytest.raises(InvalidArgument, match='interior points'):
-            op.arguments(time_M=0, y_M=22, autotune=True)
+        with pytest.raises(InvalidArgument, match='combined thickness'):
+            op.arguments(time_M=0, autotune=True)
 
         assert 'autotuning' not in op._state
 
-        op.arguments(time_M=0, y_M=23, autotune=True)
+        grid = Grid(shape=(8, 16))
+
+        kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
+                                       space_order=8) for f in fields}
+
+        op.arguments(time_M=0, autotune=True, **kwargs)
 
         assert len(op._state['autotuning']) == 1
 
@@ -2038,25 +2038,19 @@ def test_distributed_middle(self, left, mode):
 
         op, _ = self._make_operator(grid=grid, middle=True, thickness=(left, 0))
 
-        # Rank 0 has one point or an empty middle (possibly with inverted bounds).
-        # Only the global size determines whether the middle is large enough
-        if left == 25:
-            with pytest.raises(InvalidArgument, match='at least 8 interior points'):
-                op.arguments(time_M=0)
-        else:
-            op.arguments(time_M=0)
+        # Empty local middles and global middles smaller than space_order are valid
+        op.arguments(time_M=0)
 
     @pytest.mark.parallel(mode=[(2, 'basic')])
     def test_collective_rejection(self, mode):
-        grid = Grid(shape=(32, 32), topology=(2, 1))
-
-        op, _ = self._make_operator(right_shift=-4, grid=grid, middle=True)
+        grid = Grid(shape=(32, 16), topology=(2, 1))
 
-        upper = 22 if grid.distributor.myrank == 0 else 31
+        right = 9 if grid.distributor.myrank == 0 else 8
+        op, _ = self._make_operator(right_shift=-4, grid=grid, thickness=(8, right))
 
-        # Only rank 0 has an insufficient middle; its peer has 16 interior points
-        with pytest.raises(InvalidArgument, match='interior points'):
-            op.arguments(time_M=0, y_M=upper)
+        # Only rank 0 declares left/right SubDimensions too thick for the global domain
+        with pytest.raises(InvalidArgument, match='combined thickness'):
+            op.arguments(time_M=0)
 
     def test_function_on_subdomain(self):
         class Interior(SubDomain):
@@ -2084,7 +2078,7 @@ def define(self, dimensions):
         assert original.ltkn.name not in args
         assert concrete.ltkn.name in args
 
-        with pytest.raises(InvalidArgument, match='at least 8 interior points'):
+        with pytest.raises(InvalidArgument, match='Cannot override bounds'):
             op.arguments(y_M=22)
 
         with pytest.raises(InvalidArgument,
@@ -2110,6 +2104,30 @@ def test_thickness_overrides(self):
 
         assert np.all(f.data == 0)
 
+    def test_bound_overrides(self):
+        op, fields = self._make_operator()
+
+        for name, value in [('y_m', 1), ('y_M', 30), ('y', 30)]:
+            with pytest.raises(InvalidArgument, match='Cannot override bounds'):
+                op.apply(time_M=0, **{name: value})
+
+        assert all(np.all(f.data == 0) for f in fields)
+
+        # Bounds along axes without SubDimensions remain overridable
+        op.arguments(time_M=0, x_m=1, x_M=6)
+
+    @pytest.mark.parallel(mode=[(2, 'basic')])
+    def test_collective_bound_rejection(self, mode):
+        grid = Grid(shape=(8, 32), topology=(1, 2))
+
+        op, _ = self._make_operator(grid=grid)
+
+        kwargs = {'y_M': 30} if grid.distributor.myrank == 0 else {}
+
+        # One rank supplies an override; all ranks must reject it
+        with pytest.raises(InvalidArgument, match='Cannot override bounds'):
+            op.arguments(time_M=0, **kwargs)
+
     @pytest.mark.parallel(mode=[(2, 'basic')])
     def test_collective_thickness_rejection(self, mode):
         grid = Grid(shape=(32,))
@@ -2134,3 +2152,89 @@ def test_collective_thickness_rejection(self, mode):
         with pytest.raises(InvalidArgument,
                            match='Cannot override SubDimension thickness'):
             op.arguments(**kwargs)
+
+    def test_left_right_partition(self):
+        op, fields = self._make_operator()
+
+        for size in (15, 16):
+            grid = Grid(shape=(8, size))
+
+            kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
+                                           space_order=8) for f in fields}
+
+            if size < 16:
+                with pytest.raises(InvalidArgument, match='combined thickness'):
+                    op.arguments(time_M=0, **kwargs)
+            else:
+                op.arguments(time_M=0, **kwargs)
+
+    def test_largest_left_right_thickness(self):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        xl0 = SubDimension.left('xl0', x, 2)
+        xl1 = SubDimension.left('xl1', x, 8)
+        xr0 = SubDimension.right('xr0', x, 4)
+        xr1 = SubDimension.right('xr1', x, 25)
+        xm = SubDimension.middle('xm', x, 8, 8)
+
+        f = Function(name='f', grid=grid, space_order=0)
+
+        eqs = [Eq(f[d], 1) for d in (xl0, xl1, xr0, xr1, xm)]
+
+        op = Operator(eqs, name='largest_left_right_thickness')
+
+        # A valid middle SubDimension says nothing about the much larger right one
+        with pytest.raises(InvalidArgument, match='combined thickness 33'):
+            op.arguments()
+
+    def test_left_right_pairing(self):
+        grid = Grid(shape=(10, 10))
+        x, y = grid.dimensions
+        xl = SubDimension.left('xl', x, 8)
+        yr = SubDimension.right('yr', y, 8)
+        xm0 = SubDimension.middle('xm0', x, 0, 8)
+        xm1 = SubDimension.middle('xm1', x, 8, 0)
+
+        f = Function(name='f', grid=grid, space_order=0)
+
+        eqs = [Eq(f[xl, yr], 1), Eq(f[xm0, y], 2), Eq(f[xm1, y], 3)]
+
+        op = Operator(eqs, name='left_right_pairing')
+
+        # Do not pair different axes or include a middle's excluded thicknesses
+        op.arguments()
+
+    @pytest.mark.parametrize('side', ['left', 'right'])
+    def test_left_right_overlap(self, side):
+        grid = Grid(shape=(10,))
+        x, = grid.dimensions
+        xl = SubDimension.left('xl', x, 6, overlap=side == 'left')
+        xr = SubDimension.right('xr', x, 6, overlap=side == 'right')
+
+        f = Function(name='f', grid=grid, space_order=0)
+
+        eqs = [Eq(f[xl], 1), Eq(f[xr], 2)]
+
+        op = Operator(eqs, name='left_right_overlap')
+
+        op.apply()
+
+        assert np.all(f.data[:4] == 1)
+        assert np.all(f.data[4:] == 2)
+
+    @pytest.mark.parallel(mode=[(2, 'basic')])
+    def test_distributed_left_right(self, mode):
+        grid = Grid(shape=(8, 32), topology=(1, 2))
+
+        op, fields = self._make_operator(grid=grid, thickness=(12, 12))
+
+        # Left/right SubDimensions fit globally, although their sum exceeds a rank's size
+        op.arguments(time_M=0)
+
+        grid = Grid(shape=(8, 23), topology=(1, 2))
+
+        kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
+                                       space_order=8) for f in fields}
+
+        with pytest.raises(InvalidArgument, match='combined thickness 24'):
+            op.arguments(time_M=0, **kwargs)

From e92b29be9e64a8726c9f7e881393a7a7398bf0e7 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Fri, 25 Sep 2026 16:08:57 +0100
Subject: [PATCH 6/6] compiler: Strengthen SubDimensions DDA

---
 devito/ir/support/basic.py | 21 +++++++-----
 devito/types/dimension.py  | 68 +++++++++++++++++++++----------------
 tests/test_dimension.py    | 20 +++++------
 tests/test_ir.py           | 28 +++++++++++++---
 tests/test_pickle.py       | 18 +++++-----
 tests/test_subdomains.py   | 69 ++++++++++++++++++++++++++++++++------
 6 files changed, 153 insertions(+), 71 deletions(-)

diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py
index b51949d2862..5116fec5130 100644
--- a/devito/ir/support/basic.py
+++ b/devito/ir/support/basic.py
@@ -1618,11 +1618,11 @@ def disjoint_subdims(a0, a1):
     For example, a left SubDimension of thickness 4 ends before a middle
     SubDimension excluding 4 points, even when the two thickness symbols are distinct.
 
-    Left/right SubDimensions of the same parent with `overlap=False` satisfy
-    `L + R <= N`, checked against the full global parent extent at runtime.
-    Their regions may be adjacent; stencil accesses extending from one into
-    the other may still induce dependences. If either SubDimension allows
-    overlap, no minimum separation is assumed.
+    Left/right SubDimensions of the same parent with `separated=True` satisfy
+    `L + R + space_order <= N`, checked against the full global parent extent
+    at `Operator.apply`. Their gap therefore accommodates stencil accesses;
+    larger shifts are still compared explicitly. If either SubDimension has
+    `separated=False`, no minimum separation is assumed.
 
     Match data axes independently of the iteration nests. Return True if any
     axis proves separation, False otherwise. Accesses over the same interval
@@ -1670,12 +1670,15 @@ def disjoint_subdims(a0, a1):
             dl, dr = (it0.dim, it1.dim) if it0.dim.is_left else (it1.dim, it0.dim)
             dlp, drp = dl.parent, dr.parent
 
-            if dl.is_left and dr.is_right and dlp is drp and \
-               not (dl.overlap or dr.overlap):
-                # The runtime partition check guarantees a nonnegative gap only
+            if dl.is_left and dr.is_right and \
+               dl.separated and dr.separated and \
+               dlp is drp:
+                # Runtime validation guarantees N - L - R >= space_order
+                f = a0.function.c0
+                space_order = f.space_order if isinstance(f, Function) else 0
                 gap = sympy.Dummy(nonnegative=True)
                 mapper[dlp.symbolic_max] = (dlp.symbolic_min + dl.ltkn.value +
-                                            dr.rtkn.value + gap - 1)
+                                            dr.rtkn.value + space_order + gap - 1)
 
             if (M0 - m1).subs(mapper).is_negative or \
                (M1 - m0).subs(mapper).is_negative:
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index 71b484634d3..3f0fbfd18c2 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -561,14 +561,14 @@ class Thickness(DataSymbol):
     declared thickness to determine the global region before MPI decomposition.
     """
 
-    __rkwargs__ = DataSymbol.__rkwargs__ + ('root', 'side', 'local', 'value', 'overlap')
+    __rkwargs__ = DataSymbol.__rkwargs__ + ('root', 'side', 'local', 'value', 'separated')
 
-    def __new__(cls, *args, root=None, side=None, local=False, overlap=False, **kwargs):
+    def __new__(cls, *args, root=None, side=None, local=False, separated=True, **kwargs):
         newobj = super().__new__(cls, *args, **kwargs)
         newobj._root = root
         newobj._side = side
         newobj._local = local
-        newobj._overlap = overlap
+        newobj._separated = separated
 
         return newobj
 
@@ -591,8 +591,8 @@ def local(self):
         return self._local
 
     @property
-    def overlap(self):
-        return self._overlap
+    def separated(self):
+        return self._separated
 
     @property
     def value(self):
@@ -642,22 +642,22 @@ class AbstractSubDimension(DerivedDimension):
     is_AbstractSub = True
 
     __rargs__ = DerivedDimension.__rargs__ + ('thickness',)
-    __rkwargs__ = ('overlap',)
+    __rkwargs__ = ('separated',)
 
     _thickness_type = Thickness
 
-    def __init_finalize__(self, name, parent, thickness, overlap=False, **kwargs):
+    def __init_finalize__(self, name, parent, thickness, separated=True, **kwargs):
         super().__init_finalize__(name, parent)
         if parent.is_AbstractSub:
             raise ValueError("Nested SubDimensions are not supported")
 
-        self._overlap = overlap
+        self._separated = separated
 
         thickness = thickness or (None, None)
         if any(isinstance(tkn, self._thickness_type) for tkn in thickness):
             # Preserve identity when unchanged: expressions may share these symbols
             self._thickness = SubDimensionThickness(*[
-                t if t.overlap == overlap else t._rebuild(overlap=overlap)
+                t if t.separated == separated else t._rebuild(separated=separated)
                 for t in thickness
             ])
         else:
@@ -672,7 +672,7 @@ def _interval(self):
     @memoized_meth
     def _symbolic_thickness(self, **kwargs):
         kwargs = {'dtype': np.int32, 'is_const': True, 'nonnegative': True,
-                  'overlap': self.overlap}
+                  'separated': self.separated}
 
         names = [f"{self.parent.name}_{s}tkn" for s in ('l', 'r')]
         return SubDimensionThickness(*[Thickness(name=n, **kwargs) for n in names])
@@ -707,8 +707,8 @@ def rtkn(self):
         return self.thickness.right
 
     @property
-    def overlap(self):
-        return self._overlap
+    def separated(self):
+        return self._separated
 
     def __hash__(self):
         return id(self)
@@ -737,11 +737,14 @@ class SubDimension(AbstractSubDimension):
     local : bool
         True if, in case of domain decomposition, the SubDimension is
         guaranteed not to span more than one domain, False otherwise.
-    overlap : bool, optional, default=False
-        Allow a left SubDimension to overlap a right SubDimension of the same
-        parent. When both have this flag unset, their declared thicknesses must sum
-        to at most the full global parent extent. This does not guarantee
-        separation of stencil accesses extending beyond the SubDimensions themselves.
+    separated : bool, optional, default=True
+        Require a stencil-safe gap between nonempty left/right SubDimensions of
+        the same parent. If both are separated, `N - L - R >= space_order`, where
+        `N` is the global parent extent, `L` and `R` are their thicknesses, and
+        `space_order` is the maximum compiled Function order along that axis.
+        This is checked at `Operator.apply`, including runtime Grid overrides.
+        Set to False to allow touching or overlapping regions and retain
+        conservative dependence analysis. Middle SubDimensions are unaffected.
 
     Examples
     --------
@@ -781,26 +784,26 @@ def __init_finalize__(self, name, parent, thickness, local, **kwargs):
         super().__init_finalize__(name, parent, thickness, **kwargs)
 
     @classmethod
-    def left(cls, name, parent, thickness, local=True, overlap=False):
+    def left(cls, name, parent, thickness, local=True, separated=True):
         return cls(name, parent, thickness=(thickness, None), local=local,
-                   overlap=overlap)
+                   separated=separated)
 
     @classmethod
-    def right(cls, name, parent, thickness, local=True, overlap=False):
+    def right(cls, name, parent, thickness, local=True, separated=True):
         return cls(name, parent, thickness=(None, thickness), local=local,
-                   overlap=overlap)
+                   separated=separated)
 
     @classmethod
     def middle(cls, name, parent, thickness_left, thickness_right, local=False,
-               overlap=False):
+               separated=True):
         return cls(name, parent, thickness=(thickness_left, thickness_right), local=local,
-                   overlap=overlap)
+                   separated=separated)
 
     @memoized_meth
     def _symbolic_thickness(self, thickness=None):
         kwargs = {'dtype': np.int32, 'is_const': True, 'nonnegative': True,
                   'root': self.root, 'local': self.local,
-                  'overlap': self.overlap}
+                  'separated': self.separated}
 
         names = [f"{self.parent.name}_{s}tkn" for s in ('l', 'r')]
         sides = [LEFT, RIGHT]
@@ -871,9 +874,9 @@ def _arg_check(self, args, *_args, **kwargs):
     @classmethod
     def _arg_check_thickness(cls, dimensions, args):
         """
-        Check that non-overlapping left/right SubDimensions fit the global domain.
+        Check that separated left/right SubDimensions leave a stencil-safe gap.
         """
-        subdims = [d for d in dimensions if isinstance(d, cls) and not d.overlap]
+        subdims = [d for d in dimensions if isinstance(d, cls) and d.separated]
 
         grid = args.grid
         error = None
@@ -892,10 +895,17 @@ def _arg_check_thickness(cls, dimensions, args):
             # Use the full global extent, not a rank's local extent
             size = grid.size_map[root].glb if grid is not None else args[root.size_name]
             thickness = left.ltkn.value + right.rtkn.value
-            if thickness > size:
+
+            # Runtime overrides do not change the compiled stencil order
+            orders = [getattr(f, 'space_order', 0) for f in args._op_functions
+                      if root in {d.root for d in f.dimensions}]
+            gap = max(orders, default=0) if left.ltkn.value and right.rtkn.value else 0
+            if thickness + gap > size:
                 error = (f"SubDimensions `{left}` and `{right}` have combined "
-                         f"thickness {thickness} along `{parent}`, exceeding "
-                         f"the runtime extent {size}")
+                         f"thickness {thickness} along `{parent}` and require "
+                         f"a gap of at least {gap} points (space_order), exceeding "
+                         f"the runtime extent {size}; use separated=False to allow "
+                         "closer regions")
                 break
 
         comm = args.comm if args.options['mpi'] else None
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index 6bfba0b79dd..d24db35ee03 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -2257,28 +2257,28 @@ class TestConcretization:
     during compilation.
     """
 
-    @pytest.mark.parametrize('overlap', [False, True])
-    def test_correct_thicknesses(self, overlap):
+    @pytest.mark.parametrize('separated', [False, True])
+    def test_correct_thicknesses(self, separated):
         """
         Check that thicknesses aren't created where they shouldn't be.
         """
         x = Dimension('x')
-        ix0 = SubDimension.left('x', x, 2, overlap=overlap)
-        ix1 = SubDimension.right('x', x, 2, overlap=overlap)
-        ix2 = SubDimension.middle('x', x, 2, 2, overlap=overlap)
+        ix0 = SubDimension.left('x', x, 2, separated=separated)
+        ix1 = SubDimension.right('x', x, 2, separated=separated)
+        ix2 = SubDimension.middle('x', x, 2, 2, separated=separated)
 
         rebuilt = concretize_subdims([ix0, ix1, ix2], sregistry=SymbolRegistry())
 
         assert rebuilt[0].is_left
         assert rebuilt[1].is_right
         assert rebuilt[2].is_middle
-        assert all(d.overlap is overlap for d in rebuilt)
-        assert all(t.overlap is overlap for d in rebuilt for t in d.thickness)
+        assert all(d.separated is separated for d in rebuilt)
+        assert all(t.separated is separated for d in rebuilt for t in d.thickness)
 
         for d in rebuilt:
-            changed = d._rebuild(overlap=not overlap)
-            assert changed.overlap is not overlap
-            assert all(t.overlap is not overlap for t in changed.thickness)
+            changed = d._rebuild(separated=not separated)
+            assert changed.separated is not separated
+            assert all(t.separated is not separated for t in changed.thickness)
 
     def test_shared_thickness(self):
         x = Dimension('x')
diff --git a/tests/test_ir.py b/tests/test_ir.py
index a730b445f89..bf9bab86bd6 100644
--- a/tests/test_ir.py
+++ b/tests/test_ir.py
@@ -614,15 +614,15 @@ def test_opposite_subdimension_bounds(self, slope, offset, expected):
         b = TimedAccess(f[slope*(xr + offset)], 'R', 1,
                         IterationSpace([Interval(xr)]))
 
-        # Left/right SubDimensions may be adjacent, so inward reads can overlap
+        # Arrays carry no stencil order, so only a nonnegative gap is assumed
         assert a.distance(b) == b.distance(a) == (expected,)
 
     @pytest.mark.parametrize('side', ['left', 'right'])
-    def test_overlap_subdimension_bounds(self, side):
+    def test_unseparated_subdimension_bounds(self, side):
         grid = Grid(shape=(32, 32))
         x, y = grid.dimensions
-        xl = SubDimension.left('xl', x, 8, overlap=side == 'left')
-        xr = SubDimension.right('xr', x, 8, overlap=side == 'right')
+        xl = SubDimension.left('xl', x, 8, separated=side != 'left')
+        xr = SubDimension.right('xr', x, 8, separated=side != 'right')
         yl = SubDimension.left('yl', y, 4)
         ym = SubDimension.middle('ym', y, 4, 0)
 
@@ -639,6 +639,26 @@ def test_overlap_subdimension_bounds(self, side):
         # Permitting overlap along x does not manufacture a dependence along y
         assert a.distance(c) == c.distance(a) == (S.ImaginaryUnit,)
 
+    @pytest.mark.parametrize('shift,expected', [
+        (-4, S.ImaginaryUnit), (-5, S.Infinity)
+    ])
+    @pytest.mark.parametrize('bundle', [False, True])
+    def test_subdimension_stencil_gap(self, shift, expected, bundle):
+        grid = Grid(shape=(32,))
+        x, = grid.dimensions
+        xl = SubDimension.left('xl', x, 8)
+        xr = SubDimension.right('xr', x, 8)
+
+        f = Function(name='f', grid=grid, space_order=4)
+        if bundle:
+            f = Bundle(name='fg', components=(f, f.func(name='g')))
+
+        a = TimedAccess(f[xl], 'W', 0, IterationSpace([Interval(xl)]))
+        b = TimedAccess(f[xr + shift], 'R', 1, IterationSpace([Interval(xr)]))
+
+        # The promise covers the compiled order, not arbitrary larger shifts
+        assert a.distance(b) == b.distance(a) == (expected,)
+
 
 class TestSpace:
 
diff --git a/tests/test_pickle.py b/tests/test_pickle.py
index fca861c1444..7295788cbeb 100644
--- a/tests/test_pickle.py
+++ b/tests/test_pickle.py
@@ -333,7 +333,7 @@ def test_array(self, pickle):
         assert new_pa.array.name == 'a'
 
     def test_sub_dimension(self, pickle):
-        di = SubDimension.middle('di', Dimension(name='d'), 1, 1, overlap=True)
+        di = SubDimension.middle('di', Dimension(name='d'), 1, 1, separated=False)
 
         pkl_di = pickle.dumps(di)
         new_di = pickle.loads(pkl_di)
@@ -343,21 +343,21 @@ def test_sub_dimension(self, pickle):
         assert di.parent.name == new_di.parent.name
         assert di._thickness == new_di._thickness
         assert di._interval == new_di._interval
-        assert new_di.overlap
-        assert all(t.overlap for t in new_di.thickness)
+        assert not new_di.separated
+        assert all(not t.separated for t in new_di.thickness)
 
     def test_multi_sub_dimension(self, pickle):
-        di = MultiSubDimension('di', Dimension(name='d'), None, overlap=True)
+        di = MultiSubDimension('di', Dimension(name='d'), None, separated=False)
 
         new_di = pickle.loads(pickle.dumps(di))
 
-        assert new_di.overlap
-        assert all(t.overlap for t in new_di.thickness)
+        assert not new_di.separated
+        assert all(not t.separated for t in new_di.thickness)
 
-        rebuilt = new_di._rebuild(overlap=False)
+        rebuilt = new_di._rebuild(separated=True)
 
-        assert not rebuilt.overlap
-        assert all(not t.overlap for t in rebuilt.thickness)
+        assert rebuilt.separated
+        assert all(t.separated for t in rebuilt.thickness)
 
     def test_conditional_dimension(self, pickle):
         d = Dimension(name='d')
diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py
index 8f02caf841f..7e4b22fb4e9 100644
--- a/tests/test_subdomains.py
+++ b/tests/test_subdomains.py
@@ -153,7 +153,7 @@ class sd0(SubDomain):
             def define(self, dimensions):
                 x, y = dimensions
                 return {x: ('middle', 2, 2),
-                        y: SubDimension.right('iy', y, 10, overlap=True)}
+                        y: SubDimension.right('iy', y, 10, separated=False)}
 
         class sd1(SubDomain):
             name = 'sd1'
@@ -1984,6 +1984,9 @@ def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
             with pytest.raises(InvalidArgument,
                                match='Cannot override SubDimension thickness'):
                 op.arguments(time_M=0, **kwargs)
+        elif not middle and margin < 0:
+            with pytest.raises(InvalidArgument, match='gap of at least 8'):
+                op.arguments(time_M=0, **kwargs)
         else:
             op.arguments(time_M=0, **kwargs)
 
@@ -1991,7 +1994,7 @@ def test_stencil_gap(self, middle, left_shift, right_shift, margin, override):
     @pytest.mark.parametrize('space_order', [4, 8, 12])
     @pytest.mark.parametrize('margin', [-1, 0, 1])
     def test_runtime_space_order(self, middle, space_order, margin):
-        """Neither compiled nor override space_order constrains region widths."""
+        """The compiled order constrains the gap, not the runtime override's order."""
         op, fields = self._make_operator(middle=middle)
 
         required = 8
@@ -2000,7 +2003,11 @@ def test_runtime_space_order(self, middle, space_order, margin):
         kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
                                        space_order=space_order) for f in fields}
 
-        op.arguments(time_M=0, **kwargs)
+        if not middle and margin < 0:
+            with pytest.raises(InvalidArgument, match='gap of at least 8'):
+                op.arguments(time_M=0, **kwargs)
+        else:
+            op.arguments(time_M=0, **kwargs)
 
     @pytest.mark.parametrize('side', ['left', 'right'])
     def test_empty_subdimension(self, side):
@@ -2022,7 +2029,7 @@ def test_before_autotuning(self):
 
         assert 'autotuning' not in op._state
 
-        grid = Grid(shape=(8, 16))
+        grid = Grid(shape=(8, 24))
 
         kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
                                        space_order=8) for f in fields}
@@ -2043,12 +2050,12 @@ def test_distributed_middle(self, left, mode):
 
     @pytest.mark.parallel(mode=[(2, 'basic')])
     def test_collective_rejection(self, mode):
-        grid = Grid(shape=(32, 16), topology=(2, 1))
+        grid = Grid(shape=(32, 24), topology=(2, 1))
 
         right = 9 if grid.distributor.myrank == 0 else 8
         op, _ = self._make_operator(right_shift=-4, grid=grid, thickness=(8, right))
 
-        # Only rank 0 declares left/right SubDimensions too thick for the global domain
+        # Only rank 0 leaves an insufficient global stencil gap
         with pytest.raises(InvalidArgument, match='combined thickness'):
             op.arguments(time_M=0)
 
@@ -2156,13 +2163,13 @@ def test_collective_thickness_rejection(self, mode):
     def test_left_right_partition(self):
         op, fields = self._make_operator()
 
-        for size in (15, 16):
+        for size in (23, 24):
             grid = Grid(shape=(8, size))
 
             kwargs = {f.name: TimeFunction(name=f'runtime_{f.name}', grid=grid,
                                            space_order=8) for f in fields}
 
-            if size < 16:
+            if size < 24:
                 with pytest.raises(InvalidArgument, match='combined thickness'):
                     op.arguments(time_M=0, **kwargs)
             else:
@@ -2208,8 +2215,8 @@ def test_left_right_pairing(self):
     def test_left_right_overlap(self, side):
         grid = Grid(shape=(10,))
         x, = grid.dimensions
-        xl = SubDimension.left('xl', x, 6, overlap=side == 'left')
-        xr = SubDimension.right('xr', x, 6, overlap=side == 'right')
+        xl = SubDimension.left('xl', x, 6, separated=side != 'left')
+        xr = SubDimension.right('xr', x, 6, separated=side != 'right')
 
         f = Function(name='f', grid=grid, space_order=0)
 
@@ -2222,6 +2229,48 @@ def test_left_right_overlap(self, side):
         assert np.all(f.data[:4] == 1)
         assert np.all(f.data[4:] == 2)
 
+    def test_stencil_gap_axes(self):
+        grid = Grid(shape=(20, 32))
+        x, y = grid.dimensions
+        xl = SubDimension.left('xl', x, 8)
+        xr = SubDimension.right('xr', x, 8)
+
+        f = Function(name='f', grid=grid, dimensions=(x,), shape=(20,), space_order=4)
+        g = Function(name='g', grid=grid, dimensions=(y,), shape=(32,), space_order=8)
+
+        eqs = [Eq(f[xl], 1), Eq(f[xr], 2), Eq(g, g + 1)]
+
+        op = Operator(eqs, name='stencil_gap_axes')
+
+        # The higher-order Function on y must not enlarge the required x gap
+        op.apply()
+
+        assert np.all(f.data[:8] == 1)
+        assert np.all(f.data[-8:] == 2)
+
+    @pytest.mark.parametrize('space_order', [4, 8])
+    def test_stencil_gap_functions_on_subdomains(self, space_order):
+        grid = Grid(shape=(8, 20))
+        left = ReducedDomain(None, ('left', 8), grid=grid)
+        right = ReducedDomain(None, ('right', 8), grid=grid)
+
+        f = Function(name='f', grid=left, space_order=space_order)
+        g = Function(name='g', grid=right, space_order=space_order)
+
+        eqs = [Eq(f, 1), Eq(g, 2)]
+
+        op = Operator(eqs, name='stencil_gap_functions_on_subdomains')
+
+        # The SubDimension axes must contribute their Functions' stencil order
+        if space_order > 4:
+            with pytest.raises(InvalidArgument, match='gap of at least 8'):
+                op.arguments()
+        else:
+            op.apply()
+
+            assert np.all(f.data == 1)
+            assert np.all(g.data == 2)
+
     @pytest.mark.parallel(mode=[(2, 'basic')])
     def test_distributed_left_right(self, mode):
         grid = Grid(shape=(8, 32), topology=(1, 2))