From d69cc43b7c380dba637769067aed4cab1f74d62d Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Mon, 15 Jul 2024 17:04:00 +0100
Subject: [PATCH] compiler: Patch bug in SubDomainSet lowering

---
 devito/passes/clusters/implicit.py |  2 ++
 tests/test_subdomains.py           | 30 ++++++++++++++++++++++++++++++
 2 files changed, 32 insertions(+)

diff --git a/devito/passes/clusters/implicit.py b/devito/passes/clusters/implicit.py
index 590b0a0377..715a97fc27 100644
--- a/devito/passes/clusters/implicit.py
+++ b/devito/passes/clusters/implicit.py
@@ -92,6 +92,8 @@ def callback(self, clusters, prefix):
                 d = None
             if d is None:
                 processed.append(c)
+                # If no MultiSubDomain present in this cluster, then tip should be reset
+                tip = None
                 continue
 
             # Get all MultiSubDimensions in the cluster and get the dynamic thickness
diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py
index 5ee2076fd9..6ff13bc918 100644
--- a/tests/test_subdomains.py
+++ b/tests/test_subdomains.py
@@ -307,6 +307,36 @@ class MySubdomains(SubDomainSet):
         reads = set().union(*[e.reads for e in exprs])
         assert len(reads) == 4  # f, g, h, mydomains
 
+    def test_multi_eq_split(self):
+        """
+        Test cases where two loops over the same SubDomainSet will be
+        separated by another loop.
+        """
+        # Note: a bug was found where this would cause SubDomainSet
+        # bounds expressions not to be generated in the second loop over
+        # the SubDomainSet
+        class MSD(SubDomainSet):
+            name = 'msd'
+
+        msd = MSD(N=1, bounds=(1, 1, 1, 1))
+
+        grid = Grid(shape=(11, 11), subdomains=(msd,))
+
+        f = Function(name='f', grid=grid)
+        g = Function(name='g', grid=grid)
+
+        eq0 = Eq(f, 1, subdomain=msd)
+        eq1 = Eq(f, g)  # Dependency needed to fix equation order
+        eq2 = Eq(g, 1, subdomain=msd)
+
+        op = Operator([eq0, eq1, eq2])
+
+        # Ensure the loop structure is correct
+        # Note the two 'n0' correspond to the thickness definitions
+        assert_structure(op,
+                         ['n0', 'n0xy', 'xy', 'n0', 'n0xy'],
+                         'n0xyxyn0xy')
+
     def test_multi_sets(self):
         """
         Check functionality for when multiple subdomain sets are present.