From edeb9897207e263075c6cbe1434015b85fb82b18 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Fri, 17 Jul 2026 10:26:24 +0100
Subject: [PATCH 1/9] compiler: Generalize fuse-tasks grouping

---
 devito/core/operator.py             |  7 +++
 devito/passes/clusters/buffering.py | 58 +++++++++++++++++---
 tests/test_gpu_common.py            | 82 +++++++++++++++++++++++++++++
 3 files changed, 140 insertions(+), 7 deletions(-)

diff --git a/devito/core/operator.py b/devito/core/operator.py
index 2494dbd7cf..ba1ac94f5e 100644
--- a/devito/core/operator.py
+++ b/devito/core/operator.py
@@ -247,6 +247,13 @@ def _check_kwargs(cls, **kwargs):
         if oo['mpi'] and oo['mpi'] not in cls.MPI_MODES:
             raise InvalidOperator(f"Unsupported MPI mode `{oo['mpi']}`")
 
+        fuse_tasks = oo['fuse-tasks']
+        if not isinstance(fuse_tasks, bool) and \
+           (not is_integer(fuse_tasks) or fuse_tasks <= 0):
+            raise InvalidOperator(
+                "`fuse-tasks` must be a bool or a positive integer"
+            )
+
         if oo['cire-maxpar'] not in (False, 'basic', 'compact'):
             raise InvalidOperator("Illegal `cire-maxpar` value")
 
diff --git a/devito/passes/clusters/buffering.py b/devito/passes/clusters/buffering.py
index 3bf69ef1d3..8ee0c1a33f 100644
--- a/devito/passes/clusters/buffering.py
+++ b/devito/passes/clusters/buffering.py
@@ -18,6 +18,7 @@
     timed_pass
 )
 from devito.types import Array, CustomDimension, Eq, ModuloDimension
+from devito.warnings import warn
 
 __all__ = ['buffering']
 
@@ -39,7 +40,7 @@ def buffering(clusters, key, sregistry, options, **kwargs):
         The symbol registry, to create unique names for buffers and Dimensions.
     options : dict
         The optimization options.
-        Accepted: ['buf-async-degree'].
+        Accepted: ['buf-async-degree', 'buf-reuse', 'fuse-tasks'].
         * 'buf-async-degree': Specify the size of the buffer. By default, the
           buffer size is the minimal one, inferred from the memory accesses in
           the ``clusters`` themselves. An asynchronous degree equals to `k`
@@ -49,6 +50,9 @@ def buffering(clusters, key, sregistry, options, **kwargs):
           implemented by other passes).
         * 'buf-reuse': If True, the pass will try to reuse existing Buffers for
           different buffered Functions. By default, False.
+        * 'fuse-tasks': If True, fuse all asynchronous tasks. If a positive
+          integer, fuse tasks into balanced groups whose size is at most the
+          supplied value. By default, False.
     **kwargs
         Additional compilation options.
         Accepted: ['opt_init_onwrite', 'opt_buffer'].
@@ -252,15 +256,19 @@ def callback(self, clusters, prefix):
                 processed.append(Cluster(expr, ispace, guards, properties, syncs))
 
         # Lift {write,read}-only buffers into separate IterationSpaces
-        if not self.options['fuse-tasks']:
-            processed = self._optimize(processed, descriptors)
+        fuse_tasks = self.options['fuse-tasks']
+        if fuse_tasks is not True:
+            processed = self._optimize(processed, descriptors, fuse_tasks)
 
         if self.options['buf-reuse']:
             init, processed = self._reuse(init, processed, descriptors)
 
         return init + processed
 
-    def _optimize(self, clusters, descriptors):
+    def _optimize(self, clusters, descriptors, fuse_tasks=False):
+        if fuse_tasks:
+            stamps = self._make_task_groups(descriptors, int(fuse_tasks))
+
         for b, v in descriptors.items():
             if v.is_writeonly:
                 # `b` might be written by multiple, potentially mutually
@@ -268,20 +276,25 @@ def _optimize(self, clusters, descriptors):
                 # will have complementary guards, hence only one will be
                 # executed. In such a case, we can split the equations over
                 # separate IterationSpaces
-                key0 = lambda: Stamp()
+                if fuse_tasks:
+                    stamp = stamps[b]
+                    key0 = lambda: stamp  # noqa: B023
+                else:
+                    key0 = lambda: Stamp()
             elif v.is_readonly:
                 # `b` is read multiple times -- this could just be the case of
                 # coupled equations, so we more cautiously perform a
                 # "buffer-wise" splitting of the IterationSpaces (i.e., only
                 # relevant if there are at least two read-only buffers)
-                stamp = Stamp()
+                stamp = stamps[b] if fuse_tasks else Stamp()
                 key0 = lambda: stamp  # noqa: B023
             else:
                 continue
 
             processed = []
             for c in clusters:
-                if b not in c.functions:
+                function = v.f if fuse_tasks else b
+                if function not in c.functions:
                     processed.append(c)
                     continue
 
@@ -294,6 +307,37 @@ def _optimize(self, clusters, descriptors):
 
         return clusters
 
+    def _make_task_groups(self, descriptors, size):
+        stamps = {}
+
+        task_sets = (
+            [b for b, v in descriptors.items() if v.is_writeonly],
+            [b for b, v in descriptors.items() if v.is_readonly],
+        )
+        for tasks in task_sets:
+            if not tasks:
+                continue
+
+            ntasks = len(tasks)
+            ngroups = (ntasks + size - 1) // size
+            base, remainder = divmod(ntasks, ngroups)
+            sizes = (base + 1,) * remainder + (base,) * (ngroups - remainder)
+
+            if ntasks % size:
+                warn(
+                    f"`fuse-tasks={size}` does not make sense for {ntasks} tasks; "
+                    f"using balanced groups of sizes {sizes} instead"
+                )
+
+            start = 0
+            for n in sizes:
+                stamp = Stamp()
+                for b in tasks[start:start + n]:
+                    stamps[b] = stamp
+                start += n
+
+        return stamps
+
     def _reuse(self, init, clusters, descriptors):
         """
         Reuse existing Buffers for buffering candidates.
diff --git a/tests/test_gpu_common.py b/tests/test_gpu_common.py
index 079b5bdf1b..b83aa7136e 100644
--- a/tests/test_gpu_common.py
+++ b/tests/test_gpu_common.py
@@ -18,6 +18,7 @@
 )
 from devito.passes.iet.languages.openmp import OmpIteration
 from devito.types import DeviceID, DeviceRM, Lock, NPThreads, PThreadArray, Symbol
+from devito.warnings import DevitoWarning
 
 pytestmark = skipif(['nodevice'], whole_module=True)
 
@@ -678,6 +679,39 @@ def test_streaming_two_buffers(self, opt, ntmps):
         assert np.all(u.data[0] == 56)
         assert np.all(u.data[1] == 72)
 
+    def test_streaming_fuse_groups(self):
+        nt = 4
+        grid = Grid(shape=(4, 4))
+
+        u = TimeFunction(name='u', grid=grid)
+        fsaves = [TimeFunction(name=f'fsave{i}', grid=grid, save=nt)
+                  for i in range(9)]
+
+        for i, f in enumerate(fsaves):
+            for t in range(nt):
+                f.data[t, :] = i + t
+
+        eqn = Eq(u.forward, u + sum(fsaves))
+
+        op0 = Operator(eqn, opt=('noop', {'gpu-fit': tuple(fsaves)}))
+        op1 = Operator(
+            eqn,
+            opt=('buffering', 'streaming', 'fuse', 'orchestrate',
+                 {'fuse-tasks': 3})
+        )
+
+        symbols = FindSymbols().visit(op1)
+        threads = [i for i in symbols if isinstance(i, PThreadArray)]
+        assert len(threads) == 3
+        assert all(i.size == 1 for i in threads)
+
+        op0.apply(time_M=nt-2)
+
+        u1 = TimeFunction(name='u', grid=grid)
+        op1.apply(time_M=nt-2, u=u1)
+
+        assert np.all(u.data == u1.data)
+
     def test_streaming_fused(self):
         nt = 10
         grid = Grid(shape=(4, 4))
@@ -1002,6 +1036,54 @@ def test_composite_buffering_tasking_multi_output(self):
             assert np.all(usave.data == usave1.data)
             assert np.all(vsave.data == vsave1.data)
 
+    @pytest.mark.parametrize('fuse_tasks', [3, 4])
+    def test_composite_buffering_tasking_fuse_groups(self, fuse_tasks):
+        nt = 4
+        bundle0 = Bundle()
+        grid = Grid(shape=(4, 4, 4), subdomains=bundle0)
+
+        u = TimeFunction(name='u', grid=grid, time_order=2)
+        fsaves = [TimeFunction(name=f'fsave{i}', grid=grid, save=nt)
+                  for i in range(9)]
+
+        eqns = [Eq(u.forward, u + 1)]
+        eqns.extend(Eq(f, u + i, subdomain=bundle0)
+                    for i, f in enumerate(fsaves))
+
+        op0 = Operator(eqns, opt=('noop', {'gpu-fit': tuple(fsaves)}))
+
+        opt = ('buffering', 'tasking', 'topofuse', 'orchestrate',
+               {'fuse-tasks': fuse_tasks})
+        if fuse_tasks == 4:
+            with pytest.warns(DevitoWarning, match='fuse-tasks=4.*9 tasks'):
+                op1 = Operator(eqns, opt=opt)
+        else:
+            op1 = Operator(eqns, opt=opt)
+
+        symbols = FindSymbols().visit(op1)
+        threads = [i for i in symbols if isinstance(i, PThreadArray)]
+        assert len(threads) == 3
+        assert all(i.size == 1 for i in threads)
+
+        tasks = [v.root for k, v in op1._func_table.items()
+                 if k.startswith('copy_to_host')]
+        # The three pthreads reuse the same parameterized task callable
+        task, = tasks
+        writes = {i.write for i in FindNodes(Expression).visit(task)}
+        assert len([i for i in writes if i.is_TimeFunction]) == 3
+
+        op0.apply(time_M=nt-1)
+
+        u1 = TimeFunction(name='u', grid=grid, time_order=2)
+        fsaves1 = [TimeFunction(name=f'fsave{i}', grid=grid, save=nt)
+                   for i in range(9)]
+        kwargs = {f.name: f1 for f, f1 in zip(fsaves, fsaves1)}
+        op1.apply(time_M=nt-1, u=u1, **kwargs)
+
+        assert np.all(u.data == u1.data)
+        assert all(np.all(f.data == f1.data)
+                   for f, f1 in zip(fsaves, fsaves1))
+
     def test_composite_full_0(self):
         nt = 10
         grid = Grid(shape=(10, 10, 10))

From c4631a58353be902deb72dcd1443ed8f16e78c30 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Sat, 18 Jul 2026 15:50:47 +0100
Subject: [PATCH 2/9] compiler: Replace fuse-tasks with npthreads

---
 devito/core/cpu.py                  |  4 +--
 devito/core/gpu.py                  |  4 +--
 devito/core/operator.py             | 10 ++++---
 devito/passes/clusters/buffering.py | 35 ++++++++++++------------
 devito/passes/clusters/fusion.py    |  2 +-
 tests/test_gpu_common.py            | 41 ++++++++++++++---------------
 tests/test_operator.py              |  6 +++++
 7 files changed, 54 insertions(+), 48 deletions(-)

diff --git a/devito/core/cpu.py b/devito/core/cpu.py
index 79d1b409e7..3435a472e6 100644
--- a/devito/core/cpu.py
+++ b/devito/core/cpu.py
@@ -47,8 +47,8 @@ def _normalize_kwargs(cls, **kwargs):
         o['buf-async-degree'] = oo.pop('buf-async-degree', None)
         o['buf-reuse'] = oo.pop('buf-reuse', None)
 
-        # Fusion
-        o['fuse-tasks'] = oo.pop('fuse-tasks', False)
+        # Tasking
+        o['npthreads'] = oo.pop('npthreads', None)
 
         # Flops minimization
         o['cse-min-cost'] = oo.pop('cse-min-cost', cls.CSE_MIN_COST)
diff --git a/devito/core/gpu.py b/devito/core/gpu.py
index f907d4be87..0f0354f663 100644
--- a/devito/core/gpu.py
+++ b/devito/core/gpu.py
@@ -61,8 +61,8 @@ def _normalize_kwargs(cls, **kwargs):
         o['buf-async-degree'] = oo.pop('buf-async-degree', None)
         o['buf-reuse'] = oo.pop('buf-reuse', None)
 
-        # Fusion
-        o['fuse-tasks'] = oo.pop('fuse-tasks', False)
+        # Tasking
+        o['npthreads'] = oo.pop('npthreads', None)
 
         # Flops minimization
         o['cse-min-cost'] = oo.pop('cse-min-cost', cls.CSE_MIN_COST)
diff --git a/devito/core/operator.py b/devito/core/operator.py
index ba1ac94f5e..f8f34b0a24 100644
--- a/devito/core/operator.py
+++ b/devito/core/operator.py
@@ -247,11 +247,13 @@ def _check_kwargs(cls, **kwargs):
         if oo['mpi'] and oo['mpi'] not in cls.MPI_MODES:
             raise InvalidOperator(f"Unsupported MPI mode `{oo['mpi']}`")
 
-        fuse_tasks = oo['fuse-tasks']
-        if not isinstance(fuse_tasks, bool) and \
-           (not is_integer(fuse_tasks) or fuse_tasks <= 0):
+        npthreads = oo['npthreads']
+        if npthreads is not None and (
+            isinstance(npthreads, bool) or
+            not is_integer(npthreads) or npthreads <= 0
+        ):
             raise InvalidOperator(
-                "`fuse-tasks` must be a bool or a positive integer"
+                "`npthreads` must be a positive integer"
             )
 
         if oo['cire-maxpar'] not in (False, 'basic', 'compact'):
diff --git a/devito/passes/clusters/buffering.py b/devito/passes/clusters/buffering.py
index 8ee0c1a33f..0875226db1 100644
--- a/devito/passes/clusters/buffering.py
+++ b/devito/passes/clusters/buffering.py
@@ -40,7 +40,7 @@ def buffering(clusters, key, sregistry, options, **kwargs):
         The symbol registry, to create unique names for buffers and Dimensions.
     options : dict
         The optimization options.
-        Accepted: ['buf-async-degree', 'buf-reuse', 'fuse-tasks'].
+        Accepted: ['buf-async-degree', 'buf-reuse', 'npthreads'].
         * 'buf-async-degree': Specify the size of the buffer. By default, the
           buffer size is the minimal one, inferred from the memory accesses in
           the ``clusters`` themselves. An asynchronous degree equals to `k`
@@ -50,9 +50,9 @@ def buffering(clusters, key, sregistry, options, **kwargs):
           implemented by other passes).
         * 'buf-reuse': If True, the pass will try to reuse existing Buffers for
           different buffered Functions. By default, False.
-        * 'fuse-tasks': If True, fuse all asynchronous tasks. If a positive
-          integer, fuse tasks into balanced groups whose size is at most the
-          supplied value. By default, False.
+        * 'npthreads': Number of pthreads for asynchronous tasks. The tasks are
+          divided into this many balanced groups. By default, None, which uses
+          one pthread per task.
     **kwargs
         Additional compilation options.
         Accepted: ['opt_init_onwrite', 'opt_buffer'].
@@ -256,18 +256,16 @@ def callback(self, clusters, prefix):
                 processed.append(Cluster(expr, ispace, guards, properties, syncs))
 
         # Lift {write,read}-only buffers into separate IterationSpaces
-        fuse_tasks = self.options['fuse-tasks']
-        if fuse_tasks is not True:
-            processed = self._optimize(processed, descriptors, fuse_tasks)
+        processed = self._optimize(processed, descriptors, self.options['npthreads'])
 
         if self.options['buf-reuse']:
             init, processed = self._reuse(init, processed, descriptors)
 
         return init + processed
 
-    def _optimize(self, clusters, descriptors, fuse_tasks=False):
-        if fuse_tasks:
-            stamps = self._make_task_groups(descriptors, int(fuse_tasks))
+    def _optimize(self, clusters, descriptors, npthreads=None):
+        if npthreads:
+            stamps = self._make_task_groups(descriptors, npthreads)
 
         for b, v in descriptors.items():
             if v.is_writeonly:
@@ -276,7 +274,7 @@ def _optimize(self, clusters, descriptors, fuse_tasks=False):
                 # will have complementary guards, hence only one will be
                 # executed. In such a case, we can split the equations over
                 # separate IterationSpaces
-                if fuse_tasks:
+                if npthreads:
                     stamp = stamps[b]
                     key0 = lambda: stamp  # noqa: B023
                 else:
@@ -286,14 +284,14 @@ def _optimize(self, clusters, descriptors, fuse_tasks=False):
                 # coupled equations, so we more cautiously perform a
                 # "buffer-wise" splitting of the IterationSpaces (i.e., only
                 # relevant if there are at least two read-only buffers)
-                stamp = stamps[b] if fuse_tasks else Stamp()
+                stamp = stamps[b] if npthreads else Stamp()
                 key0 = lambda: stamp  # noqa: B023
             else:
                 continue
 
             processed = []
             for c in clusters:
-                function = v.f if fuse_tasks else b
+                function = v.f if npthreads else b
                 if function not in c.functions:
                     processed.append(c)
                     continue
@@ -307,7 +305,8 @@ def _optimize(self, clusters, descriptors, fuse_tasks=False):
 
         return clusters
 
-    def _make_task_groups(self, descriptors, size):
+    def _make_task_groups(self, descriptors, npthreads):
+        """Assign task buffers to at most ``npthreads`` balanced groups."""
         stamps = {}
 
         task_sets = (
@@ -319,14 +318,14 @@ def _make_task_groups(self, descriptors, size):
                 continue
 
             ntasks = len(tasks)
-            ngroups = (ntasks + size - 1) // size
+            ngroups = min(npthreads, ntasks)
             base, remainder = divmod(ntasks, ngroups)
             sizes = (base + 1,) * remainder + (base,) * (ngroups - remainder)
 
-            if ntasks % size:
+            if npthreads > ntasks:
                 warn(
-                    f"`fuse-tasks={size}` does not make sense for {ntasks} tasks; "
-                    f"using balanced groups of sizes {sizes} instead"
+                    f"`npthreads={npthreads}` exceeds the {ntasks} available "
+                    f"tasks; using `npthreads={ntasks}` instead"
                 )
 
             start = 0
diff --git a/devito/passes/clusters/fusion.py b/devito/passes/clusters/fusion.py
index fcc3d9affa..bbd4694614 100644
--- a/devito/passes/clusters/fusion.py
+++ b/devito/passes/clusters/fusion.py
@@ -152,7 +152,7 @@ def __init__(self, toposort, options=None):
         options = options or {}
 
         self.toposort = toposort
-        self.fuse_tasks = options.get('fuse-tasks', False)
+        self.fuse_tasks = options.get('npthreads') is not None
 
         super().__init__()
 
diff --git a/tests/test_gpu_common.py b/tests/test_gpu_common.py
index b83aa7136e..8b499deb2a 100644
--- a/tests/test_gpu_common.py
+++ b/tests/test_gpu_common.py
@@ -18,7 +18,6 @@
 )
 from devito.passes.iet.languages.openmp import OmpIteration
 from devito.types import DeviceID, DeviceRM, Lock, NPThreads, PThreadArray, Symbol
-from devito.warnings import DevitoWarning
 
 pytestmark = skipif(['nodevice'], whole_module=True)
 
@@ -497,7 +496,7 @@ def test_tasking_forcefuse(self):
                 Eq(v.forward, tmp1, subdomain=bundle0)]
 
         op = Operator(eqns, opt=('tasking', 'fuse', 'orchestrate',
-                                 {'fuse-tasks': True, 'linearize': False}))
+                                 {'npthreads': 1, 'linearize': False}))
 
         # Check generated code
         assert len(retrieve_iteration_tree(op)) == 3
@@ -647,7 +646,7 @@ def test_streaming_basic(self, opt, ntmps):
 
     @pytest.mark.parametrize('opt,ntmps', [
         (('buffering', 'streaming', 'orchestrate'), 14),
-        (('buffering', 'streaming', 'fuse', 'orchestrate', {'fuse-tasks': True}), 8),
+        (('buffering', 'streaming', 'fuse', 'orchestrate', {'npthreads': 1}), 8),
     ])
     def test_streaming_two_buffers(self, opt, ntmps):
         nt = 10
@@ -679,7 +678,8 @@ def test_streaming_two_buffers(self, opt, ntmps):
         assert np.all(u.data[0] == 56)
         assert np.all(u.data[1] == 72)
 
-    def test_streaming_fuse_groups(self):
+    @pytest.mark.parametrize('npthreads', [1, 3, 4])
+    def test_streaming_fuse_groups(self, npthreads):
         nt = 4
         grid = Grid(shape=(4, 4))
 
@@ -697,13 +697,13 @@ def test_streaming_fuse_groups(self):
         op1 = Operator(
             eqn,
             opt=('buffering', 'streaming', 'fuse', 'orchestrate',
-                 {'fuse-tasks': 3})
+                 {'npthreads': npthreads})
         )
 
         symbols = FindSymbols().visit(op1)
         threads = [i for i in symbols if isinstance(i, PThreadArray)]
-        assert len(threads) == 3
         assert all(i.size == 1 for i in threads)
+        assert op1.npthreads == npthreads
 
         op0.apply(time_M=nt-2)
 
@@ -995,7 +995,7 @@ def test_composite_buffering_tasking_multi_output(self):
                                   {'buf-async-degree': async_degree}))
         op2 = Operator(eqns, opt=('buffering', 'tasking', 'topofuse', 'orchestrate',
                                   {'buf-async-degree': async_degree,
-                                   'fuse-tasks': True}))
+                                   'npthreads': 1}))
 
         # Check generated code -- thanks to buffering only expect 1 lock!
         assert len(retrieve_iteration_tree(op0)) == 2
@@ -1036,8 +1036,8 @@ def test_composite_buffering_tasking_multi_output(self):
             assert np.all(usave.data == usave1.data)
             assert np.all(vsave.data == vsave1.data)
 
-    @pytest.mark.parametrize('fuse_tasks', [3, 4])
-    def test_composite_buffering_tasking_fuse_groups(self, fuse_tasks):
+    @pytest.mark.parametrize('npthreads', [1, 3, 4])
+    def test_composite_buffering_tasking_fuse_groups(self, npthreads):
         nt = 4
         bundle0 = Bundle()
         grid = Grid(shape=(4, 4, 4), subdomains=bundle0)
@@ -1053,24 +1053,23 @@ def test_composite_buffering_tasking_fuse_groups(self, fuse_tasks):
         op0 = Operator(eqns, opt=('noop', {'gpu-fit': tuple(fsaves)}))
 
         opt = ('buffering', 'tasking', 'topofuse', 'orchestrate',
-               {'fuse-tasks': fuse_tasks})
-        if fuse_tasks == 4:
-            with pytest.warns(DevitoWarning, match='fuse-tasks=4.*9 tasks'):
-                op1 = Operator(eqns, opt=opt)
-        else:
-            op1 = Operator(eqns, opt=opt)
+               {'npthreads': npthreads})
+        op1 = Operator(eqns, opt=opt)
 
         symbols = FindSymbols().visit(op1)
         threads = [i for i in symbols if isinstance(i, PThreadArray)]
-        assert len(threads) == 3
         assert all(i.size == 1 for i in threads)
+        assert op1.npthreads == npthreads
 
         tasks = [v.root for k, v in op1._func_table.items()
                  if k.startswith('copy_to_host')]
-        # The three pthreads reuse the same parameterized task callable
-        task, = tasks
-        writes = {i.write for i in FindNodes(Expression).visit(task)}
-        assert len([i for i in writes if i.is_TimeFunction]) == 3
+        # Pthreads handling the same number of tasks reuse a Callable
+        widths = []
+        for task in tasks:
+            writes = {i.write for i in FindNodes(Expression).visit(task)}
+            widths.append(len([i for i in writes if i.is_TimeFunction]))
+        expected = {1: [9], 3: [3], 4: [2, 3]}
+        assert sorted(widths) == expected[npthreads]
 
         op0.apply(time_M=nt-1)
 
@@ -1409,7 +1408,7 @@ def test_streaming_complete(self):
         op1 = Operator(eqns, opt=('buffering', 'streaming', 'orchestrate'))
         op2 = Operator(eqns, opt=('buffering', 'streaming', 'fuse', 'orchestrate'))
         op3 = Operator(eqns, opt=('buffering', 'streaming', 'fuse', 'orchestrate',
-                                  {'fuse-tasks': True}))
+                                  {'npthreads': 1}))
 
         # Check generated code
         diff = int(configuration['language'] == 'openmp')
diff --git a/tests/test_operator.py b/tests/test_operator.py
index 806b091e24..b3f7ea295e 100644
--- a/tests/test_operator.py
+++ b/tests/test_operator.py
@@ -135,6 +135,12 @@ def test_opt_options(self):
         except InvalidOperator:
             assert True
 
+        # `npthreads` is a direct count, not a boolean switch
+        for npthreads in (False, True, 0, -1, 1.5):
+            with pytest.raises(InvalidOperator, match='positive integer'):
+                Operator(Eq(u, u + 1),
+                         opt=('advanced', {'npthreads': npthreads}))
+
     def test_compiler_uniqueness(self):
         grid = Grid(shape=(3, 3, 3))
 

From 8ad3d2aceb718c74f6bd45ae882525c5d11a6430 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Mon, 20 Jul 2026 09:02:38 +0100
Subject: [PATCH 3/9] compiler: Polish buffering

---
 devito/passes/clusters/buffering.py | 49 ++++++++++++++++++-----------
 1 file changed, 30 insertions(+), 19 deletions(-)

diff --git a/devito/passes/clusters/buffering.py b/devito/passes/clusters/buffering.py
index 0875226db1..5f72b8d746 100644
--- a/devito/passes/clusters/buffering.py
+++ b/devito/passes/clusters/buffering.py
@@ -256,16 +256,17 @@ def callback(self, clusters, prefix):
                 processed.append(Cluster(expr, ispace, guards, properties, syncs))
 
         # Lift {write,read}-only buffers into separate IterationSpaces
-        processed = self._optimize(processed, descriptors, self.options['npthreads'])
+        processed = self._optimize(processed, descriptors)
 
-        if self.options['buf-reuse']:
-            init, processed = self._reuse(init, processed, descriptors)
+        # Reuse existing Buffers for buffering candidates, if requested
+        init, processed = self._reuse(init, processed, descriptors)
 
         return init + processed
 
-    def _optimize(self, clusters, descriptors, npthreads=None):
-        if npthreads:
-            stamps = self._make_task_groups(descriptors, npthreads)
+    def _optimize(self, clusters, descriptors):
+        npthreads = self.options['npthreads']
+
+        stamps = self._make_task_groups(descriptors)
 
         for b, v in descriptors.items():
             if v.is_writeonly:
@@ -274,25 +275,22 @@ def _optimize(self, clusters, descriptors, npthreads=None):
                 # will have complementary guards, hence only one will be
                 # executed. In such a case, we can split the equations over
                 # separate IterationSpaces
-                if npthreads:
-                    stamp = stamps[b]
-                    key0 = lambda: stamp  # noqa: B023
-                else:
-                    key0 = lambda: Stamp()
+                key0 = lambda: stamps[b] if npthreads else Stamp()  # noqa: B023
             elif v.is_readonly:
                 # `b` is read multiple times -- this could just be the case of
                 # coupled equations, so we more cautiously perform a
                 # "buffer-wise" splitting of the IterationSpaces (i.e., only
                 # relevant if there are at least two read-only buffers)
-                stamp = stamps[b] if npthreads else Stamp()
-                key0 = lambda: stamp  # noqa: B023
+                stamp_fixed = Stamp()
+                key0 = lambda: stamps[b] if npthreads else stamp_fixed  # noqa: B023
             else:
                 continue
 
             processed = []
             for c in clusters:
-                function = v.f if npthreads else b
-                if function not in c.functions:
+                f = v.f if npthreads else b
+
+                if f not in c.functions:
                     processed.append(c)
                     continue
 
@@ -305,18 +303,27 @@ def _optimize(self, clusters, descriptors, npthreads=None):
 
         return clusters
 
-    def _make_task_groups(self, descriptors, npthreads):
-        """Assign task buffers to at most ``npthreads`` balanced groups."""
+    def _make_task_groups(self, descriptors):
+        """
+        Assign task buffers to at most `npthreads` balanced groups.
+        """
+        npthreads = self.options['npthreads']
+
         stamps = {}
+        if not npthreads:
+            return stamps
 
         task_sets = (
             [b for b, v in descriptors.items() if v.is_writeonly],
             [b for b, v in descriptors.items() if v.is_readonly],
         )
+
         for tasks in task_sets:
             if not tasks:
                 continue
 
+            # Calculate the number of groups and their sizes, so that the tasks
+            # are divided into balanced groups
             ntasks = len(tasks)
             ngroups = min(npthreads, ntasks)
             base, remainder = divmod(ntasks, ngroups)
@@ -328,11 +335,12 @@ def _make_task_groups(self, descriptors, npthreads):
                     f"tasks; using `npthreads={ntasks}` instead"
                 )
 
+            # Create a unique Stamp for each group and assign it to the tasks
+            # in that group
             start = 0
             for n in sizes:
                 stamp = Stamp()
-                for b in tasks[start:start + n]:
-                    stamps[b] = stamp
+                stamps.update({b: stamp for b in tasks[start:start + n]})
                 start += n
 
         return stamps
@@ -343,6 +351,9 @@ def _reuse(self, init, clusters, descriptors):
         """
         buf_reuse = self.options['buf-reuse']
 
+        if not buf_reuse:
+            return init, clusters
+
         if callable(buf_reuse):
             cbk = lambda v: [i for i in v if buf_reuse(descriptors[i])]
         else:

From 26b03bcc8fc04cfb388fbceb04b238ab256be742 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Mon, 20 Jul 2026 13:46:25 +0100
Subject: [PATCH 4/9] compiler: Restructure SyncSpot

---
 devito/ir/iet/algorithms.py          |  2 +-
 devito/ir/iet/nodes.py               | 16 +++++++++-------
 devito/ir/support/syncs.py           | 13 +++++++++----
 devito/passes/clusters/asynchrony.py | 19 +++++++++++++++++--
 devito/passes/iet/orchestration.py   |  4 ++--
 5 files changed, 38 insertions(+), 16 deletions(-)

diff --git a/devito/ir/iet/algorithms.py b/devito/ir/iet/algorithms.py
index e805c34f42..847b206581 100644
--- a/devito/ir/iet/algorithms.py
+++ b/devito/ir/iet/algorithms.py
@@ -58,7 +58,7 @@ def iet_build(stree):
                 body = HaloSpot(None, i.halo_scheme)
 
         elif i.is_Sync:
-            body = SyncSpot(i.sync_ops, body=queues.pop(i, None))
+            body = SyncSpot((i.sync_ops,), body=queues.pop(i, None))
 
         queues.setdefault(i.parent, []).append(body)
 
diff --git a/devito/ir/iet/nodes.py b/devito/ir/iet/nodes.py
index bc4b5dc48c..9b36972478 100644
--- a/devito/ir/iet/nodes.py
+++ b/devito/ir/iet/nodes.py
@@ -1511,18 +1511,20 @@ class BusyWait(While):
 class SyncSpot(List):
 
     """
-    A node representing one or more synchronization operations, e.g., WaitLock,
-    withLock, etc.
+    A node coupling synchronization operations with an IET body.
+
+    `ops` contains either one group applying to the whole `body`, or one group
+    per body item, with `ops[i]` applying to `body[i]`.
     """
 
     is_SyncSpot = True
 
-    def __init__(self, sync_ops, body=None):
+    def __init__(self, ops, body=None):
         super().__init__(body=body)
-        self.sync_ops = sync_ops
+        self.ops = tuple(tuple(i) for i in ops)
 
     def __repr__(self):
-        return f"<SyncSpot ({','.join(str(i) for i in self.sync_ops)})>"
+        return f"<SyncSpot ({','.join(str(i) for i in flatten(self.ops))})>"
 
     @property
     def is_async_op(self):
@@ -1531,11 +1533,11 @@ def is_async_op(self):
         If False, the SyncSpot may for example represent a wait on a lock.
         """
         return any(isinstance(s, (WithLock, PrefetchUpdate))
-                   for s in self.sync_ops)
+                   for s in flatten(self.ops))
 
     @property
     def functions(self):
-        ret = [(s.lock, s.function, s.target) for s in self.sync_ops]
+        ret = [(s.lock, s.function, s.target) for s in flatten(self.ops)]
         ret = tuple(filter_ordered(f for f in flatten(ret) if f is not None))
         return ret
 
diff --git a/devito/ir/support/syncs.py b/devito/ir/support/syncs.py
index d7f4f86bfe..3a39a417dd 100644
--- a/devito/ir/support/syncs.py
+++ b/devito/ir/support/syncs.py
@@ -27,10 +27,12 @@
 class SyncOp(Pickable):
 
     __rargs__ = ('handle', 'target')
-    __rkwargs__ = ('tindex', 'function', 'findex', 'dim', 'size', 'origin')
+    __rkwargs__ = (
+        'tindex', 'function', 'findex', 'dim', 'size', 'origin', 'gid'
+    )
 
     def __init__(self, handle, target, tindex=None, function=None, findex=None,
-                 dim=None, size=1, origin=None):
+                 dim=None, size=1, origin=None, gid=None):
         self.handle = handle
         self.target = target
 
@@ -40,6 +42,7 @@ def __init__(self, handle, target, tindex=None, function=None, findex=None,
         self.dim = dim
         self.size = size
         self.origin = origin
+        self.gid = gid
 
     def __eq__(self, other):
         return (type(self) is type(other) and
@@ -50,11 +53,13 @@ def __eq__(self, other):
                 self.findex == other.findex and
                 self.dim is other.dim and
                 self.size == other.size and
-                self.origin == other.origin)
+                self.origin == other.origin and
+                self.gid == other.gid)
 
     def __hash__(self):
         return hash((self.__class__, self.handle, self.target, self.tindex,
-                     self.function, self.findex, self.dim, self.size, self.origin))
+                     self.function, self.findex, self.dim, self.size, self.origin,
+                     self.gid))
 
     def __repr__(self):
         return f"{self.__class__.__name__}<{self.handle.name}>"
diff --git a/devito/passes/clusters/asynchrony.py b/devito/passes/clusters/asynchrony.py
index e32190ddef..7a92484680 100644
--- a/devito/passes/clusters/asynchrony.py
+++ b/devito/passes/clusters/asynchrony.py
@@ -45,6 +45,16 @@ def memcpy_key(c):
     return task_key, memcpy_key
 
 
+def _task_gid(ispace, d):
+    """Return the group ID carried by the non-trigger Intervals."""
+    gids = {i.stamp for i in ispace if not i.dim._defines & d._defines}
+    if len(gids) != 1:
+        return None
+
+    gid, = gids
+    return gid
+
+
 @timed_pass(name='tasking')
 def tasking(clusters, key0, sregistry):
     """
@@ -147,6 +157,8 @@ def _schedule_waitlocks(self, c0, d, clusters, locks, syncs):
         return protected
 
     def _schedule_withlocks(self, c0, d, protected, locks, syncs):
+        gid = _task_gid(c0.ispace, d)
+
         for target in protected:
             lock = locks[target]
 
@@ -169,7 +181,7 @@ def _schedule_withlocks(self, c0, d, protected, locks, syncs):
             for i in indices:
                 syncs[c0][d].update([
                     ReleaseLock(lock[i], target),
-                    WithLock(lock[i], target, i, function, findex, d)
+                    WithLock(lock[i], target, i, function, findex, d, gid=gid)
                 ])
 
 
@@ -274,9 +286,12 @@ def _actions_from_update_memcpy(c, d, clusters, actions, sregistry):
     guard1 = GuardBoundNext(function.indices[d], direction)
     guards = c.guards.impose(d, guard0 & guard1)
 
+    gid = _task_gid(ispace, d)
+
     syncs = {d: [
         ReleaseLock(handle, target),
-        PrefetchUpdate(handle, target, tindex, function, findex, d, 1, e.rhs)
+        PrefetchUpdate(handle, target, tindex, function, findex, d,
+                       origin=e.rhs, gid=gid)
     ]}
     syncs = {**c.syncs, **syncs}
 
diff --git a/devito/passes/iet/orchestration.py b/devito/passes/iet/orchestration.py
index 3bea70fb0f..a697abf11e 100644
--- a/devito/passes/iet/orchestration.py
+++ b/devito/passes/iet/orchestration.py
@@ -15,7 +15,7 @@
 from devito.passes.iet.engine import iet_pass
 from devito.passes.iet.langbase import LangBB
 from devito.symbolics import CondEq, CondNe
-from devito.tools import DAG, as_mapper, as_tuple
+from devito.tools import DAG, as_mapper, as_tuple, flatten
 from devito.types import HostLayer
 
 __init__ = ['Orchestrator']
@@ -151,7 +151,7 @@ def process(self, iet):
                 break
 
             n0 = ordered(sync_spots).pop(0)
-            mapper = as_mapper(n0.sync_ops, lambda i: type(i))
+            mapper = as_mapper(flatten(n0.ops), lambda i: type(i))
 
             subs = {}
             for t in sorted(mapper, key=key):

From 9b38f496506c614e70158b50293ff50ed5fbc046 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Mon, 20 Jul 2026 17:14:17 +0100
Subject: [PATCH 5/9] compiler: Refactor fuse-tasks lowering

---
 devito/ir/iet/algorithms.py        |   2 +-
 devito/ir/iet/nodes.py             |  37 +++++--
 devito/passes/iet/orchestration.py | 171 ++++++++++++++++++++++++-----
 3 files changed, 171 insertions(+), 39 deletions(-)

diff --git a/devito/ir/iet/algorithms.py b/devito/ir/iet/algorithms.py
index 847b206581..e805c34f42 100644
--- a/devito/ir/iet/algorithms.py
+++ b/devito/ir/iet/algorithms.py
@@ -58,7 +58,7 @@ def iet_build(stree):
                 body = HaloSpot(None, i.halo_scheme)
 
         elif i.is_Sync:
-            body = SyncSpot((i.sync_ops,), body=queues.pop(i, None))
+            body = SyncSpot(i.sync_ops, body=queues.pop(i, None))
 
         queues.setdefault(i.parent, []).append(body)
 
diff --git a/devito/ir/iet/nodes.py b/devito/ir/iet/nodes.py
index 9b36972478..ff3ea9004a 100644
--- a/devito/ir/iet/nodes.py
+++ b/devito/ir/iet/nodes.py
@@ -61,6 +61,7 @@
     'Section',
     'Switch',
     'SyncSpot',
+    'SyncSpotRegion',
     'TimedList',
     'Transfer',
     'Using',
@@ -92,6 +93,7 @@ class Node(Signer):
     is_HaloSpot = False
     is_ExpressionBundle = False
     is_SyncSpot = False
+    is_SyncSpotRegion = False
 
     _traversable = []
     """
@@ -1512,19 +1514,16 @@ class SyncSpot(List):
 
     """
     A node coupling synchronization operations with an IET body.
-
-    `ops` contains either one group applying to the whole `body`, or one group
-    per body item, with `ops[i]` applying to `body[i]`.
     """
 
     is_SyncSpot = True
 
-    def __init__(self, ops, body=None):
+    def __init__(self, sync_ops, body=None):
         super().__init__(body=body)
-        self.ops = tuple(tuple(i) for i in ops)
+        self.sync_ops = sync_ops
 
     def __repr__(self):
-        return f"<SyncSpot ({','.join(str(i) for i in flatten(self.ops))})>"
+        return f"<SyncSpot ({','.join(str(i) for i in self.sync_ops)})>"
 
     @property
     def is_async_op(self):
@@ -1533,15 +1532,37 @@ def is_async_op(self):
         If False, the SyncSpot may for example represent a wait on a lock.
         """
         return any(isinstance(s, (WithLock, PrefetchUpdate))
-                   for s in flatten(self.ops))
+                   for s in self.sync_ops)
 
     @property
     def functions(self):
-        ret = [(s.lock, s.function, s.target) for s in flatten(self.ops)]
+        ret = [(s.lock, s.function, s.target) for s in self.sync_ops]
         ret = tuple(filter_ordered(f for f in flatten(ret) if f is not None))
         return ret
 
 
+class SyncSpotRegion(List):
+
+    """A sequence of SyncSpots treated as one unit during orchestration."""
+
+    is_SyncSpotRegion = True
+
+    def __init__(self, body):
+        body = as_tuple(body)
+        assert body and all(isinstance(i, SyncSpot) for i in body)
+        super().__init__(body=body)
+
+    @property
+    def sync_spots(self):
+        return self.body
+
+    @cached_property
+    def optype(self):
+        """The common type of the synchronization operations in this region."""
+        optype, = {type(j) for i in self.sync_spots for j in i.sync_ops}
+        return optype
+
+
 class CBlankLine(List):
 
     def __init__(self, **kwargs):
diff --git a/devito/passes/iet/orchestration.py b/devito/passes/iet/orchestration.py
index a697abf11e..e360c86162 100644
--- a/devito/passes/iet/orchestration.py
+++ b/devito/passes/iet/orchestration.py
@@ -1,4 +1,4 @@
-from collections import OrderedDict
+from collections import OrderedDict, defaultdict, namedtuple
 from contextlib import suppress
 from functools import singledispatch
 
@@ -6,25 +6,32 @@
 
 from devito.exceptions import CompilationError
 from devito.ir.iet import (
-    AsyncCall, AsyncCallable, BlankLine, BusyWait, Call, Callable, DummyExpr, FindNodes,
-    List, SyncSpot, Transformer, derive_parameters, make_callable
+    AsyncCall, AsyncCallable, BlankLine, BusyWait, Call, Callable, Conditional,
+    DummyExpr, FindNodes, List, SyncSpot, SyncSpotRegion, Transformer,
+    derive_parameters, make_callable
 )
+from devito.ir.iet.visitors import LazyVisitor
 from devito.ir.support import (
     InitArray, PrefetchUpdate, ReleaseLock, SnapIn, SnapOut, SyncArray, WaitLock, WithLock
 )
 from devito.passes.iet.engine import iet_pass
 from devito.passes.iet.langbase import LangBB
 from devito.symbolics import CondEq, CondNe
-from devito.tools import DAG, as_mapper, as_tuple, flatten
+from devito.tools import DAG, as_mapper
 from devito.types import HostLayer
 
 __init__ = ['Orchestrator']
 
 
+_Task = namedtuple(
+    '_Task', 'spot condition anchor iteration gid task_type sync_ops releases'
+)
+
+
 class Orchestrator:
 
     """
-    Lower the SyncSpot in IET for efficient host-device asynchronous computation.
+    Lower synchronization nodes for efficient host-device asynchronous computation.
     """
 
     langbb = LangBB
@@ -32,8 +39,47 @@ class Orchestrator:
     The language used to implement host-device data movements.
     """
 
-    def __init__(self, sregistry=None, **kwargs):
+    def __init__(self, sregistry=None, options=None, **kwargs):
         self.sregistry = sregistry
+        self.npthreads = (options or {}).get('npthreads')
+
+    def _fuse_tasks(self, iet):
+        """Group compatible task SyncSpots into SyncSpotRegions."""
+        mapper = as_mapper(
+            _FindTasks().visit(iet),
+            lambda task: (task.iteration, task.gid, task.task_type,
+                          isinstance(task.anchor, SyncSpot))
+        )
+
+        insertions = defaultdict(list)
+        subs1 = {}
+
+        for tasks in mapper.values():
+            spots = [SyncSpot(task.sync_ops,
+                              body=Conditional(task.condition.condition,
+                                               task.spot.body))
+                     for task in tasks]
+            insertions[tasks[-1].anchor].append(SyncSpotRegion(spots))
+
+            for task in tasks:
+                subs1[task.spot] = SyncSpot(task.releases)
+
+        if not subs1:
+            return iet
+
+        # These substitutions cannot be merged because a Transformer does not
+        # revisit a replacement, while task spots may be nested below an anchor
+        subs0 = {}
+        for anchor, regions in insertions.items():
+            if isinstance(anchor, SyncSpot):
+                subs0[anchor] = anchor._rebuild(body=anchor.body + tuple(regions))
+            else:
+                subs0[anchor] = List(body=(anchor, *regions))
+
+        iet = Transformer(subs0).visit(iet)
+        iet = Transformer(subs1).visit(iet)
+
+        return iet
 
     def _make_waitlock(self, iet, sync_ops, *args):
         waitloop = List(
@@ -53,6 +99,9 @@ def _make_releaselock(self, iet, sync_ops, *args):
         parameters = derive_parameters(pre, ordering='canonical')
         efunc = Callable(name, pre, 'void', parameters, 'static')
 
+        if isinstance(iet, SyncSpot) and not iet.body:
+            iet = List()
+
         iet = List(body=[Call(name, efunc.parameters)] + [iet])
 
         return iet, [efunc]
@@ -60,17 +109,15 @@ def _make_releaselock(self, iet, sync_ops, *args):
     def _make_withlock(self, iet, sync_ops, layer):
         body, prefix = withlock(layer, iet, sync_ops, self.langbb, self.sregistry)
 
-        # Turn `iet` into an AsyncCallable so that subsequent passes know
-        # that we're happy for this Callable to be executed asynchronously
+        return self._make_async_callable(body, prefix)
+
+    def _make_async_callable(self, body, prefix):
         name = self.sregistry.make_name(prefix=prefix)
         body = List(body=body)
         parameters = derive_parameters(body, ordering='canonical')
         efunc = AsyncCallable(name, body, parameters=parameters)
 
-        # The corresponding AsyncCall
-        iet = AsyncCall(name, efunc.parameters)
-
-        return iet, [efunc]
+        return AsyncCall(name, efunc.parameters), [efunc]
 
     def _make_callable(self, name, iet, *args):
         name = self.sregistry.make_name(prefix=name)
@@ -108,23 +155,36 @@ def _make_syncarray(self, iet, sync_ops, layer):
     def _make_prefetchupdate(self, iet, sync_ops, layer):
         body, prefix = prefetchupdate(layer, iet, sync_ops, self.langbb, self.sregistry)
 
-        # Turn `iet` into an AsyncCallable so that subsequent passes know
-        # that we're happy for this Callable to be executed asynchronously
-        name = self.sregistry.make_name(prefix=prefix)
-        body = List(body=body)
-        parameters = derive_parameters(body, ordering='canonical')
-        efunc = AsyncCallable(name, body, parameters=parameters)
+        return self._make_async_callable(body, prefix)
 
-        # The corresponding AsyncCall
-        iet = AsyncCall(name, efunc.parameters)
+    def _make_region(self, iet):
+        """Lower a SyncSpotRegion into one asynchronous callable."""
+        sync_ops = tuple(j for i in iet.sync_spots for j in i.sync_ops)
+        callback = {
+            PrefetchUpdate: prefetchupdate,
+            WithLock: withlock
+        }[iet.optype]
 
-        return iet, [efunc]
+        layers = {infer_layer(i.function) for i in sync_ops}
+        if len(layers) != 1:
+            raise CompilationError("Unsupported streaming case")
+        layer = layers.pop()
+
+        body = []
+        for spot in iet.sync_spots:
+            condition, = spot.body
+            task_body, prefix = callback(
+                layer, List(body=condition.then_body), spot.sync_ops, self.langbb,
+                self.sregistry
+            )
+            body.append(condition._rebuild(then_body=task_body))
+
+        return self._make_async_callable(body, prefix)
 
     @iet_pass
     def process(self, iet):
-        sync_spots = FindNodes(SyncSpot).visit(iet)
-        if not sync_spots:
-            return iet, {}
+        if self.npthreads is not None:
+            iet = self._fuse_tasks(iet)
 
         # The SyncOps are to be processed in a given order
         callbacks = OrderedDict([
@@ -139,19 +199,32 @@ def process(self, iet):
         ])
         key = lambda s: list(callbacks).index(s)
 
+        # A region consumes its immediate SyncSpots as one unit, so lower all
+        # regions before looking for ordinary SyncSpots
+        efuncs = []
+        while True:
+            sync_regions = FindNodes(SyncSpotRegion).visit(iet)
+            if not sync_regions:
+                break
+
+            n0 = ordered(sync_regions).pop(0)
+            n1, v = self._make_region(n0)
+
+            iet = Transformer({n0: n1}).visit(iet)
+            efuncs.extend(v)
+
         # The SyncSpots may be nested, so we compute a topological ordering
         # so that they are processed in a bottom-up fashion. This is necessary
         # because e.g. an inner SyncSpot may generate new objects (e.g., a new
         # Queue), which in turn must be visible to the outer SyncSpot to
         # generate the correct parameters list
-        efuncs = []
         while True:
             sync_spots = FindNodes(SyncSpot).visit(iet)
             if not sync_spots:
                 break
 
             n0 = ordered(sync_spots).pop(0)
-            mapper = as_mapper(flatten(n0.ops), lambda i: type(i))
+            mapper = as_mapper(n0.sync_ops, type)
 
             subs = {}
             for t in sorted(mapper, key=key):
@@ -172,15 +245,53 @@ def process(self, iet):
         return iet, {'efuncs': efuncs}
 
 
-def ordered(sync_spots):
-    dag = DAG(nodes=sync_spots)
-    for n0 in sync_spots:
-        for n1 in as_tuple(FindNodes(SyncSpot).visit(n0.body)):
+def ordered(sync_nodes):
+    dag = DAG(nodes=sync_nodes)
+    for n0 in sync_nodes:
+        for n1 in FindNodes(type(n0)).visit(n0.body):
             dag.add_edge(n1, n0)
 
     return dag.topological_sort()
 
 
+class _FindTasks(LazyVisitor):
+
+    """Find guarded asynchronous SyncSpots and their structural context."""
+
+    def visit_Iteration(self, o, **kwargs):
+        kwargs['iteration'] = o
+        yield from self._visit(o.children, **kwargs)
+
+    def visit_Conditional(self, o, **kwargs):
+        kwargs['condition'] = o
+        yield from self._visit(o.children, **kwargs)
+
+    def visit_SyncSpot(self, o, iteration=None, condition=None, anchor=None):
+        if iteration is not None and condition is not None and not condition.else_body:
+            syncs = as_mapper(o.sync_ops, type)
+            sync_types = set(syncs)
+            for task_type in (PrefetchUpdate, WithLock):
+                if sync_types != {task_type, ReleaseLock}:
+                    continue
+
+                sync_ops = syncs[task_type]
+                gids = {i.gid for i in sync_ops}
+
+                if len(gids) == 1 and None not in gids:
+                    gid, = gids
+                    yield _Task(o, condition, anchor or condition, iteration,
+                                gid, task_type, sync_ops, syncs[ReleaseLock])
+                break
+
+        if any(isinstance(i, SnapOut) for i in o.sync_ops):
+            # Keep composite task calls inside the `SnapOut` scope
+            anchor = o
+
+        yield from self._visit(
+            o.children, iteration=iteration, condition=condition, anchor=anchor
+        )
+
+
 # Task handlers
 
 layer_host = HostLayer('host')

From 7fe79ff87ffb3c20bbb6f79d82a985d03c0c03a4 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Mon, 20 Jul 2026 17:54:45 +0100
Subject: [PATCH 6/9] compiler: Simplify npthreads lowering

---
 devito/ir/iet/nodes.py             |  3 +-
 devito/passes/iet/orchestration.py | 74 +++++++++++++++---------------
 2 files changed, 38 insertions(+), 39 deletions(-)

diff --git a/devito/ir/iet/nodes.py b/devito/ir/iet/nodes.py
index ff3ea9004a..fac1e42e7b 100644
--- a/devito/ir/iet/nodes.py
+++ b/devito/ir/iet/nodes.py
@@ -1559,7 +1559,8 @@ def sync_spots(self):
     @cached_property
     def optype(self):
         """The common type of the synchronization operations in this region."""
-        optype, = {type(j) for i in self.sync_spots for j in i.sync_ops}
+        optype, = {type(op) for spot in self.sync_spots
+                   for op in spot.sync_ops}
         return optype
 
 
diff --git a/devito/passes/iet/orchestration.py b/devito/passes/iet/orchestration.py
index e360c86162..a86f768b57 100644
--- a/devito/passes/iet/orchestration.py
+++ b/devito/passes/iet/orchestration.py
@@ -6,9 +6,9 @@
 
 from devito.exceptions import CompilationError
 from devito.ir.iet import (
-    AsyncCall, AsyncCallable, BlankLine, BusyWait, Call, Callable, Conditional,
-    DummyExpr, FindNodes, List, SyncSpot, SyncSpotRegion, Transformer,
-    derive_parameters, make_callable
+    AsyncCall, AsyncCallable, BlankLine, BusyWait, Call, Callable, Conditional, DummyExpr,
+    FindNodes, List, SyncSpot, SyncSpotRegion, Transformer, derive_parameters,
+    make_callable
 )
 from devito.ir.iet.visitors import LazyVisitor
 from devito.ir.support import (
@@ -20,11 +20,11 @@
 from devito.tools import DAG, as_mapper
 from devito.types import HostLayer
 
-__init__ = ['Orchestrator']
+__all__ = ['Orchestrator']
 
 
-_Task = namedtuple(
-    '_Task', 'spot condition anchor iteration gid task_type sync_ops releases'
+TaskMeta = namedtuple(
+    'TaskMeta', 'spot condition anchor iteration gid optype sync_ops releases'
 )
 
 
@@ -45,16 +45,16 @@ def __init__(self, sregistry=None, options=None, **kwargs):
 
     def _fuse_tasks(self, iet):
         """Group compatible task SyncSpots into SyncSpotRegions."""
-        mapper = as_mapper(
-            _FindTasks().visit(iet),
-            lambda task: (task.iteration, task.gid, task.task_type,
+        groups = as_mapper(
+            CollectTasks().visit(iet),
+            lambda task: (task.iteration, task.gid, task.optype,
                           isinstance(task.anchor, SyncSpot))
         )
 
         insertions = defaultdict(list)
         subs1 = {}
 
-        for tasks in mapper.values():
+        for tasks in groups.values():
             spots = [SyncSpot(task.sync_ops,
                               body=Conditional(task.condition.condition,
                                                task.spot.body))
@@ -159,7 +159,8 @@ def _make_prefetchupdate(self, iet, sync_ops, layer):
 
     def _make_region(self, iet):
         """Lower a SyncSpotRegion into one asynchronous callable."""
-        sync_ops = tuple(j for i in iet.sync_spots for j in i.sync_ops)
+        sync_ops = tuple(op for spot in iet.sync_spots
+                         for op in spot.sync_ops)
         callback = {
             PrefetchUpdate: prefetchupdate,
             WithLock: withlock
@@ -186,6 +187,17 @@ def process(self, iet):
         if self.npthreads is not None:
             iet = self._fuse_tasks(iet)
 
+        # Lower regions first so their member SyncSpots are not processed
+        # independently by the generic SyncSpot lowering below
+        efuncs = []
+        subs = {}
+        for region in FindNodes(SyncSpotRegion).visit(iet):
+            call, efuncs1 = self._make_region(region)
+            subs[region] = call
+            efuncs.extend(efuncs1)
+
+        iet = Transformer(subs).visit(iet)
+
         # The SyncOps are to be processed in a given order
         callbacks = OrderedDict([
             (WaitLock, self._make_waitlock),
@@ -197,21 +209,7 @@ def process(self, iet):
             (PrefetchUpdate, self._make_prefetchupdate),
             (ReleaseLock, self._make_releaselock),
         ])
-        key = lambda s: list(callbacks).index(s)
-
-        # A region consumes its immediate SyncSpots as one unit, so lower all
-        # regions before looking for ordinary SyncSpots
-        efuncs = []
-        while True:
-            sync_regions = FindNodes(SyncSpotRegion).visit(iet)
-            if not sync_regions:
-                break
-
-            n0 = ordered(sync_regions).pop(0)
-            n1, v = self._make_region(n0)
-
-            iet = Transformer({n0: n1}).visit(iet)
-            efuncs.extend(v)
+        key = tuple(callbacks).index
 
         # The SyncSpots may be nested, so we compute a topological ordering
         # so that they are processed in a bottom-up fashion. This is necessary
@@ -245,18 +243,18 @@ def process(self, iet):
         return iet, {'efuncs': efuncs}
 
 
-def ordered(sync_nodes):
-    dag = DAG(nodes=sync_nodes)
-    for n0 in sync_nodes:
-        for n1 in FindNodes(type(n0)).visit(n0.body):
+def ordered(sync_spots):
+    dag = DAG(nodes=sync_spots)
+    for n0 in sync_spots:
+        for n1 in FindNodes(SyncSpot).visit(n0.body):
             dag.add_edge(n1, n0)
 
     return dag.topological_sort()
 
 
-class _FindTasks(LazyVisitor):
+class CollectTasks(LazyVisitor):
 
-    """Find guarded asynchronous SyncSpots and their structural context."""
+    """Collect guarded asynchronous SyncSpots and their structural metadata."""
 
     def visit_Iteration(self, o, **kwargs):
         kwargs['iteration'] = o
@@ -269,18 +267,18 @@ def visit_Conditional(self, o, **kwargs):
     def visit_SyncSpot(self, o, iteration=None, condition=None, anchor=None):
         if iteration is not None and condition is not None and not condition.else_body:
             syncs = as_mapper(o.sync_ops, type)
-            sync_types = set(syncs)
-            for task_type in (PrefetchUpdate, WithLock):
-                if sync_types != {task_type, ReleaseLock}:
+            optypes = set(syncs)
+            for optype in (PrefetchUpdate, WithLock):
+                if optypes != {optype, ReleaseLock}:
                     continue
 
-                sync_ops = syncs[task_type]
+                sync_ops = syncs[optype]
                 gids = {i.gid for i in sync_ops}
 
                 if len(gids) == 1 and None not in gids:
                     gid, = gids
-                    yield _Task(o, condition, anchor or condition, iteration,
-                                gid, task_type, sync_ops, syncs[ReleaseLock])
+                    yield TaskMeta(o, condition, anchor or condition, iteration,
+                                   gid, optype, sync_ops, syncs[ReleaseLock])
                 break
 
         if any(isinstance(i, SnapOut) for i in o.sync_ops):

From a5a6e1f68a1aa78c5f1d1fde09fc8ca6d20613e2 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Tue, 21 Jul 2026 09:22:41 +0100
Subject: [PATCH 7/9] compiler: Simplify task group collection

---
 devito/passes/clusters/asynchrony.py | 15 ++++----
 devito/passes/iet/orchestration.py   | 56 ++++++++++++++++------------
 2 files changed, 41 insertions(+), 30 deletions(-)

diff --git a/devito/passes/clusters/asynchrony.py b/devito/passes/clusters/asynchrony.py
index 7a92484680..e37e334927 100644
--- a/devito/passes/clusters/asynchrony.py
+++ b/devito/passes/clusters/asynchrony.py
@@ -45,14 +45,15 @@ def memcpy_key(c):
     return task_key, memcpy_key
 
 
-def _task_gid(ispace, d):
-    """Return the group ID carried by the non-trigger Intervals."""
+def make_gid(ispace, d):
+    """
+    Make a group ID from the stamps carried by the non-trigger Intervals.
+    """
     gids = {i.stamp for i in ispace if not i.dim._defines & d._defines}
     if len(gids) != 1:
         return None
-
-    gid, = gids
-    return gid
+    else:
+        return gids.pop()
 
 
 @timed_pass(name='tasking')
@@ -157,7 +158,7 @@ def _schedule_waitlocks(self, c0, d, clusters, locks, syncs):
         return protected
 
     def _schedule_withlocks(self, c0, d, protected, locks, syncs):
-        gid = _task_gid(c0.ispace, d)
+        gid = make_gid(c0.ispace, d)
 
         for target in protected:
             lock = locks[target]
@@ -286,7 +287,7 @@ def _actions_from_update_memcpy(c, d, clusters, actions, sregistry):
     guard1 = GuardBoundNext(function.indices[d], direction)
     guards = c.guards.impose(d, guard0 & guard1)
 
-    gid = _task_gid(ispace, d)
+    gid = make_gid(ispace, d)
 
     syncs = {d: [
         ReleaseLock(handle, target),
diff --git a/devito/passes/iet/orchestration.py b/devito/passes/iet/orchestration.py
index a86f768b57..2d421fd41c 100644
--- a/devito/passes/iet/orchestration.py
+++ b/devito/passes/iet/orchestration.py
@@ -23,9 +23,7 @@
 __all__ = ['Orchestrator']
 
 
-TaskMeta = namedtuple(
-    'TaskMeta', 'spot condition anchor iteration gid optype sync_ops releases'
-)
+Task = namedtuple('Task', 'spot guard anchor sync_ops releases')
 
 
 class Orchestrator:
@@ -44,25 +42,24 @@ def __init__(self, sregistry=None, options=None, **kwargs):
         self.npthreads = (options or {}).get('npthreads')
 
     def _fuse_tasks(self, iet):
-        """Group compatible task SyncSpots into SyncSpotRegions."""
-        groups = as_mapper(
-            CollectTasks().visit(iet),
-            lambda task: (task.iteration, task.gid, task.optype,
-                          isinstance(task.anchor, SyncSpot))
-        )
+        """
+        Group compatible task SyncSpots into SyncSpotRegions.
+        """
+        if self.npthreads is None:
+            return iet
+
+        groups = CollectTasks().visit(iet)
 
         insertions = defaultdict(list)
         subs1 = {}
 
         for tasks in groups.values():
             spots = [SyncSpot(task.sync_ops,
-                              body=Conditional(task.condition.condition,
-                                               task.spot.body))
+                              body=Conditional(task.guard, task.spot.body))
                      for task in tasks]
-            insertions[tasks[-1].anchor].append(SyncSpotRegion(spots))
 
-            for task in tasks:
-                subs1[task.spot] = SyncSpot(task.releases)
+            insertions[tasks[-1].anchor].append(SyncSpotRegion(spots))
+            subs1.update({task.spot: SyncSpot(task.releases) for task in tasks})
 
         if not subs1:
             return iet
@@ -184,8 +181,9 @@ def _make_region(self, iet):
 
     @iet_pass
     def process(self, iet):
-        if self.npthreads is not None:
-            iet = self._fuse_tasks(iet)
+        # Group compatible task SyncSpots into SyncSpotRegions, if requested
+        # by the user
+        iet = self._fuse_tasks(iet)
 
         # Lower regions first so their member SyncSpots are not processed
         # independently by the generic SyncSpot lowering below
@@ -254,7 +252,13 @@ def ordered(sync_spots):
 
 class CollectTasks(LazyVisitor):
 
-    """Collect guarded asynchronous SyncSpots and their structural metadata."""
+    """Collect and group guarded asynchronous SyncSpots."""
+
+    def _post_visit(self, ret):
+        groups = defaultdict(list)
+        for key, task in ret:
+            groups[key].append(task)
+        return groups
 
     def visit_Iteration(self, o, **kwargs):
         kwargs['iteration'] = o
@@ -265,20 +269,26 @@ def visit_Conditional(self, o, **kwargs):
         yield from self._visit(o.children, **kwargs)
 
     def visit_SyncSpot(self, o, iteration=None, condition=None, anchor=None):
-        if iteration is not None and condition is not None and not condition.else_body:
+        if iteration is not None and condition is not None:
             syncs = as_mapper(o.sync_ops, type)
+
             optypes = set(syncs)
             for optype in (PrefetchUpdate, WithLock):
                 if optypes != {optype, ReleaseLock}:
                     continue
 
+                # Task SyncSpots inherit a guard without an `else` branch from
+                # the originating Cluster
+                assert not condition.else_body
+
                 sync_ops = syncs[optype]
-                gids = {i.gid for i in sync_ops}
 
-                if len(gids) == 1 and None not in gids:
-                    gid, = gids
-                    yield TaskMeta(o, condition, anchor or condition, iteration,
-                                   gid, optype, sync_ops, syncs[ReleaseLock])
+                gid, = {i.gid for i in sync_ops}
+                if gid is not None:
+                    task = Task(o, condition.condition, anchor or condition,
+                                sync_ops, syncs[ReleaseLock])
+                    yield (iteration, gid, optype, anchor is not None), task
+
                 break
 
         if any(isinstance(i, SnapOut) for i in o.sync_ops):

From cae7e6a79898643923be8336693c7d0ab9f2d413 Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Tue, 21 Jul 2026 10:24:42 +0100
Subject: [PATCH 8/9] compiler: Support unguarded task groups

---
 devito/passes/iet/orchestration.py | 36 +++++++++++++++++++++---------
 1 file changed, 25 insertions(+), 11 deletions(-)

diff --git a/devito/passes/iet/orchestration.py b/devito/passes/iet/orchestration.py
index 2d421fd41c..7b856157bf 100644
--- a/devito/passes/iet/orchestration.py
+++ b/devito/passes/iet/orchestration.py
@@ -6,8 +6,8 @@
 
 from devito.exceptions import CompilationError
 from devito.ir.iet import (
-    AsyncCall, AsyncCallable, BlankLine, BusyWait, Call, Callable, Conditional, DummyExpr,
-    FindNodes, List, SyncSpot, SyncSpotRegion, Transformer, derive_parameters,
+    AsyncCall, AsyncCallable, BlankLine, Block, BusyWait, Call, Callable, Conditional,
+    DummyExpr, FindNodes, List, SyncSpot, SyncSpotRegion, Transformer, derive_parameters,
     make_callable
 )
 from devito.ir.iet.visitors import LazyVisitor
@@ -54,12 +54,21 @@ def _fuse_tasks(self, iet):
         subs1 = {}
 
         for tasks in groups.values():
+            subs1.update({task.spot: SyncSpot(task.releases) for task in tasks})
+
+            # Unguarded tasks form separate groups and can share one SyncSpot
+            if tasks[0].guard is None:
+                sync_ops = tuple(op for task in tasks for op in task.sync_ops)
+                sync_ops += tuple(tasks[-1].releases)
+                # Blocks preserve the local scope of each task body
+                body = [Block(body=task.spot.body) for task in tasks]
+                subs1[tasks[-1].spot] = SyncSpot(sync_ops, body=body)
+                continue
+
             spots = [SyncSpot(task.sync_ops,
                               body=Conditional(task.guard, task.spot.body))
                      for task in tasks]
-
             insertions[tasks[-1].anchor].append(SyncSpotRegion(spots))
-            subs1.update({task.spot: SyncSpot(task.releases) for task in tasks})
 
         if not subs1:
             return iet
@@ -252,7 +261,7 @@ def ordered(sync_spots):
 
 class CollectTasks(LazyVisitor):
 
-    """Collect and group guarded asynchronous SyncSpots."""
+    """Collect and group compatible asynchronous SyncSpots."""
 
     def _post_visit(self, ret):
         groups = defaultdict(list)
@@ -269,7 +278,7 @@ def visit_Conditional(self, o, **kwargs):
         yield from self._visit(o.children, **kwargs)
 
     def visit_SyncSpot(self, o, iteration=None, condition=None, anchor=None):
-        if iteration is not None and condition is not None:
+        if iteration is not None:
             syncs = as_mapper(o.sync_ops, type)
 
             optypes = set(syncs)
@@ -277,17 +286,22 @@ def visit_SyncSpot(self, o, iteration=None, condition=None, anchor=None):
                 if optypes != {optype, ReleaseLock}:
                     continue
 
-                # Task SyncSpots inherit a guard without an `else` branch from
-                # the originating Cluster
-                assert not condition.else_body
+                guard = None
+                if condition is not None:
+                    # Task SyncSpots inherit a guard without an `else` branch
+                    # from the originating Cluster
+                    assert not condition.else_body
+                    guard = condition.condition
 
                 sync_ops = syncs[optype]
 
                 gid, = {i.gid for i in sync_ops}
                 if gid is not None:
-                    task = Task(o, condition.condition, anchor or condition,
+                    task = Task(o, guard, anchor or condition,
                                 sync_ops, syncs[ReleaseLock])
-                    yield (iteration, gid, optype, anchor is not None), task
+                    key = (iteration, gid, optype, anchor is not None,
+                           condition is not None)
+                    yield key, task
 
                 break
 

From 3afe7c7172f4620b75bdd1c45935b63fe135a16a Mon Sep 17 00:00:00 2001
From: Fabio Luporini <fabio@devitocodes.com>
Date: Tue, 21 Jul 2026 14:00:54 +0100
Subject: [PATCH 9/9] compiler: Unify task group lowering

---
 devito/ir/iet/nodes.py             |  25 ------
 devito/ir/support/syncs.py         |  29 +++++++
 devito/passes/iet/orchestration.py | 129 ++++++++++++++++++-----------
 tests/test_gpu_common.py           |   4 +-
 4 files changed, 111 insertions(+), 76 deletions(-)

diff --git a/devito/ir/iet/nodes.py b/devito/ir/iet/nodes.py
index fac1e42e7b..4e63438dad 100644
--- a/devito/ir/iet/nodes.py
+++ b/devito/ir/iet/nodes.py
@@ -61,7 +61,6 @@
     'Section',
     'Switch',
     'SyncSpot',
-    'SyncSpotRegion',
     'TimedList',
     'Transfer',
     'Using',
@@ -93,7 +92,6 @@ class Node(Signer):
     is_HaloSpot = False
     is_ExpressionBundle = False
     is_SyncSpot = False
-    is_SyncSpotRegion = False
 
     _traversable = []
     """
@@ -1541,29 +1539,6 @@ def functions(self):
         return ret
 
 
-class SyncSpotRegion(List):
-
-    """A sequence of SyncSpots treated as one unit during orchestration."""
-
-    is_SyncSpotRegion = True
-
-    def __init__(self, body):
-        body = as_tuple(body)
-        assert body and all(isinstance(i, SyncSpot) for i in body)
-        super().__init__(body=body)
-
-    @property
-    def sync_spots(self):
-        return self.body
-
-    @cached_property
-    def optype(self):
-        """The common type of the synchronization operations in this region."""
-        optype, = {type(op) for spot in self.sync_spots
-                   for op in spot.sync_ops}
-        return optype
-
-
 class CBlankLine(List):
 
     def __init__(self, **kwargs):
diff --git a/devito/ir/support/syncs.py b/devito/ir/support/syncs.py
index 3a39a417dd..57280fa28f 100644
--- a/devito/ir/support/syncs.py
+++ b/devito/ir/support/syncs.py
@@ -26,6 +26,35 @@
 
 class SyncOp(Pickable):
 
+    """
+    Metadata for a synchronization operation attached to a `Cluster` or `SyncSpot`.
+
+    Parameters
+    ----------
+    handle : object
+        The symbolic object identifying or controlling the operation, such as
+        an entry in a `Lock`. May be None when no handle is required.
+    target : AbstractFunction
+        The `Function` whose access is synchronized. For buffered data movements,
+        this is the compiler-generated buffer.
+    tindex : Expr, optional
+        The index into `target` involved in the operation.
+    function : AbstractFunction, optional
+        The original `Function` represented by a compiler-generated `target`. It
+        is the source of a `SyncCopyIn` and the destination of a `SyncCopyOut`.
+    findex : Expr, optional
+        The index into `function` corresponding to `tindex`.
+    dim : Dimension, optional
+        The `Dimension` along which `tindex` and `findex` are defined.
+    size : int, optional
+        The extent associated with the operation along `dim`. Defaults to 1.
+    origin : Indexed, optional
+        The original `Indexed` access from which the operation was derived.
+    gid : Stamp, optional
+        The `Stamp` identifying the asynchronous task group to which the operation
+        belongs.
+    """
+
     __rargs__ = ('handle', 'target')
     __rkwargs__ = (
         'tindex', 'function', 'findex', 'dim', 'size', 'origin', 'gid'
diff --git a/devito/passes/iet/orchestration.py b/devito/passes/iet/orchestration.py
index 7b856157bf..52c0821992 100644
--- a/devito/passes/iet/orchestration.py
+++ b/devito/passes/iet/orchestration.py
@@ -1,14 +1,13 @@
 from collections import OrderedDict, defaultdict, namedtuple
 from contextlib import suppress
-from functools import singledispatch
+from functools import cached_property, singledispatch
 
 from sympy import Or
 
 from devito.exceptions import CompilationError
 from devito.ir.iet import (
     AsyncCall, AsyncCallable, BlankLine, Block, BusyWait, Call, Callable, Conditional,
-    DummyExpr, FindNodes, List, SyncSpot, SyncSpotRegion, Transformer, derive_parameters,
-    make_callable
+    DummyExpr, FindNodes, List, SyncSpot, Transformer, derive_parameters, make_callable
 )
 from devito.ir.iet.visitors import LazyVisitor
 from devito.ir.support import (
@@ -23,7 +22,32 @@
 __all__ = ['Orchestrator']
 
 
-Task = namedtuple('Task', 'spot guard anchor sync_ops releases')
+Task = namedtuple('Task', 'spot guard sync_ops releases')
+
+
+class SyncSpotRegion(List):
+
+    """A sequence of SyncSpots treated as one unit during orchestration."""
+
+    def __init__(self, body):
+        super().__init__(body=body)
+        assert self.body and all(isinstance(i, SyncSpot) for i in self.body)
+
+    @property
+    def sync_spots(self):
+        return self.body
+
+    @cached_property
+    def optype(self):
+        """
+        The type of the synchronization operation in this region.
+        """
+        optypes = {type(op) for spot in self.sync_spots for op in spot.sync_ops}
+        assert len(optypes) == 1, (
+            "Expected a SyncSpotRegion to contain exactly one type of "
+            "synchronization operation"
+        )
+        return optypes.pop()
 
 
 class Orchestrator:
@@ -48,39 +72,29 @@ def _fuse_tasks(self, iet):
         if self.npthreads is None:
             return iet
 
-        groups = CollectTasks().visit(iet)
-
         insertions = defaultdict(list)
         subs1 = {}
 
-        for tasks in groups.values():
+        for tasks, anchor in CollectTasks().visit(iet):
+            spots = []
+            for task in tasks:
+                if task.guard is None:
+                    # Preserve a local scope when there is no Conditional
+                    scope = Block(body=task.spot.body)
+                else:
+                    scope = Conditional(task.guard, task.spot.body)
+                spots.append(SyncSpot(task.sync_ops, body=scope))
+
+            insertions[anchor].append(SyncSpotRegion(spots))
             subs1.update({task.spot: SyncSpot(task.releases) for task in tasks})
 
-            # Unguarded tasks form separate groups and can share one SyncSpot
-            if tasks[0].guard is None:
-                sync_ops = tuple(op for task in tasks for op in task.sync_ops)
-                sync_ops += tuple(tasks[-1].releases)
-                # Blocks preserve the local scope of each task body
-                body = [Block(body=task.spot.body) for task in tasks]
-                subs1[tasks[-1].spot] = SyncSpot(sync_ops, body=body)
-                continue
-
-            spots = [SyncSpot(task.sync_ops,
-                              body=Conditional(task.guard, task.spot.body))
-                     for task in tasks]
-            insertions[tasks[-1].anchor].append(SyncSpotRegion(spots))
-
         if not subs1:
             return iet
 
         # These substitutions cannot be merged because a Transformer does not
         # revisit a replacement, while task spots may be nested below an anchor
-        subs0 = {}
-        for anchor, regions in insertions.items():
-            if isinstance(anchor, SyncSpot):
-                subs0[anchor] = anchor._rebuild(body=anchor.body + tuple(regions))
-            else:
-                subs0[anchor] = List(body=(anchor, *regions))
+        subs0 = {anchor: List(body=(anchor, *regions))
+                 for anchor, regions in insertions.items()}
 
         iet = Transformer(subs0).visit(iet)
         iet = Transformer(subs1).visit(iet)
@@ -172,19 +186,22 @@ def _make_region(self, iet):
             WithLock: withlock
         }[iet.optype]
 
-        layers = {infer_layer(i.function) for i in sync_ops}
-        if len(layers) != 1:
-            raise CompilationError("Unsupported streaming case")
-        layer = layers.pop()
+        layer = infer_sync_layer(sync_ops)
 
         body = []
         for spot in iet.sync_spots:
-            condition, = spot.body
+            scope, = spot.body
+            if isinstance(scope, Conditional):
+                task_body = scope.then_body
+            else:
+                assert isinstance(scope, Block)
+                task_body = scope.body
+
             task_body, prefix = callback(
-                layer, List(body=condition.then_body), spot.sync_ops, self.langbb,
+                layer, List(body=task_body), spot.sync_ops, self.langbb,
                 self.sregistry
             )
-            body.append(condition._rebuild(then_body=task_body))
+            body.append(scope._rebuild(task_body))
 
         return self._make_async_callable(body, prefix)
 
@@ -235,10 +252,7 @@ def process(self, iet):
             for t in sorted(mapper, key=key):
                 sync_ops = mapper[t]
 
-                layers = {infer_layer(s.function) for s in sync_ops}
-                if len(layers) != 1:
-                    raise CompilationError("Unsupported streaming case")
-                layer = layers.pop()
+                layer = infer_sync_layer(sync_ops)
 
                 n1, v = callbacks[t](subs.get(n0, n0), sync_ops, layer)
 
@@ -265,9 +279,12 @@ class CollectTasks(LazyVisitor):
 
     def _post_visit(self, ret):
         groups = defaultdict(list)
-        for key, task in ret:
+        anchors = {}
+        for key, task, anchor in ret:
             groups[key].append(task)
-        return groups
+            # Activate the group after its last task in program order
+            anchors[key] = anchor
+        return tuple((tasks, anchors[key]) for key, tasks in groups.items())
 
     def visit_Iteration(self, o, **kwargs):
         kwargs['iteration'] = o
@@ -277,7 +294,8 @@ def visit_Conditional(self, o, **kwargs):
         kwargs['condition'] = o
         yield from self._visit(o.children, **kwargs)
 
-    def visit_SyncSpot(self, o, iteration=None, condition=None, anchor=None):
+    def visit_SyncSpot(self, o, iteration=None, condition=None,
+                       in_snapshot=False):
         if iteration is not None:
             syncs = as_mapper(o.sync_ops, type)
 
@@ -297,20 +315,19 @@ def visit_SyncSpot(self, o, iteration=None, condition=None, anchor=None):
 
                 gid, = {i.gid for i in sync_ops}
                 if gid is not None:
-                    task = Task(o, guard, anchor or condition,
-                                sync_ops, syncs[ReleaseLock])
-                    key = (iteration, gid, optype, anchor is not None,
-                           condition is not None)
-                    yield key, task
+                    task = Task(o, guard, sync_ops, syncs[ReleaseLock])
+                    key = (iteration, gid, optype, in_snapshot)
+                    yield key, task, condition or o
 
                 break
 
         if any(isinstance(i, SnapOut) for i in o.sync_ops):
-            # Keep composite task calls inside the `SnapOut` scope
-            anchor = o
+            # Do not mix composite tasks with other compatible groups
+            in_snapshot = True
 
         yield from self._visit(
-            o.children, iteration=iteration, condition=condition, anchor=anchor
+            o.children, iteration=iteration, condition=condition,
+            in_snapshot=in_snapshot
         )
 
 
@@ -327,6 +344,20 @@ def infer_layer(f):
     return layer_host
 
 
+def infer_sync_layer(sync_ops):
+    """
+    Infer the unique storage layer used by a sequence of SyncOps.
+    """
+    layers = {infer_layer(i.function) for i in sync_ops}
+    if len(layers) != 1:
+        found = ', '.join(sorted(str(i) for i in layers)) or 'none'
+        raise CompilationError(
+            "Expected synchronization operations to use exactly one storage "
+            f"layer, but found: {found}"
+        )
+    return layers.pop()
+
+
 @singledispatch
 def withlock(layer, iet, sync_ops, lang, sregistry):
     raise NotImplementedError
diff --git a/tests/test_gpu_common.py b/tests/test_gpu_common.py
index 8b499deb2a..e3297f221d 100644
--- a/tests/test_gpu_common.py
+++ b/tests/test_gpu_common.py
@@ -1076,12 +1076,12 @@ def test_composite_buffering_tasking_fuse_groups(self, npthreads):
         u1 = TimeFunction(name='u', grid=grid, time_order=2)
         fsaves1 = [TimeFunction(name=f'fsave{i}', grid=grid, save=nt)
                    for i in range(9)]
-        kwargs = {f.name: f1 for f, f1 in zip(fsaves, fsaves1)}
+        kwargs = {f.name: f1 for f, f1 in zip(fsaves, fsaves1, strict=True)}
         op1.apply(time_M=nt-1, u=u1, **kwargs)
 
         assert np.all(u.data == u1.data)
         assert all(np.all(f.data == f1.data)
-                   for f, f1 in zip(fsaves, fsaves1))
+                   for f, f1 in zip(fsaves, fsaves1, strict=True))
 
     def test_composite_full_0(self):
         nt = 10