operator: Fix SubDomain set bug. See issue #1453.

devito/operator/operator.py

@@ -254,7 +254,6 @@ def _add_implicit(cls, expressions):
         but instead are requisites of some specified functionality.
         """
         processed = []
-        seen = set()
         for e in expressions:
             if e.subdomain:
                 try:
@@ -263,11 +262,9 @@ def _add_implicit(cls, expressions):
                     sub_dims.append(e.subdomain.implicit_dimension)
                     dims = [d for d in dims if d not in frozenset(sub_dims)]
                     dims.append(e.subdomain.implicit_dimension)
-                    if e.subdomain not in seen:
-                        grid = list(retrieve_functions(e, mode='unique'))[0].grid
-                        processed.extend([i.func(*i.args, implicit_dims=dims) for i in
-                                          e.subdomain._create_implicit_exprs(grid)])
-                        seen.add(e.subdomain)
+                    grid = list(retrieve_functions(e, mode='unique'))[0].grid
+                    processed.extend([i.func(*i.args, implicit_dims=dims) for i in
+                                      e.subdomain._create_implicit_exprs(grid)])
                     dims.extend(e.subdomain.dimensions)
                     new_e = Eq(e.lhs, e.rhs, subdomain=e.subdomain, implicit_dims=dims)
                     processed.append(new_e)

tests/test_subdomains.py

@@ -356,3 +356,65 @@ class Inner(SubDomainSet):
         fex.data[:] = np.transpose(expected)
 
         assert((np.array(result) == np.array(fex.data[:])).all())
+
+    def test_multi_sets_eq(self):
+        """
+        Check functionality for when multiple subdomain sets are present, each
+        with multiple equations.
+        """
+
+        Nx = 10
+        Ny = Nx
+        n_domains = 2
+
+        n = Dimension(name='n')
+        m = Dimension(name='m')
+
+        class MySubdomains1(SubDomainSet):
+            name = 'mydomains1'
+            implicit_dimension = n
+
+        class MySubdomains2(SubDomainSet):
+            name = 'mydomains2'
+            implicit_dimension = m
+
+        bounds_xm = np.array([1, Nx/2+1], dtype=np.int32)
+        bounds_xM = np.array([Nx/2+1, 1], dtype=np.int32)
+        bounds_ym = int(1)
+        bounds_yM = int(Ny/2+1)
+        bounds1 = (bounds_xm, bounds_xM, bounds_ym, bounds_yM)
+
+        bounds_xm = np.array([1, Nx/2+1], dtype=np.int32)
+        bounds_xM = np.array([Nx/2+1, 1], dtype=np.int32)
+        bounds_ym = int(Ny/2+1)
+        bounds_yM = int(1)
+        bounds2 = (bounds_xm, bounds_xM, bounds_ym, bounds_yM)
+
+        my_sd1 = MySubdomains1(N=n_domains, bounds=bounds1)
+        my_sd2 = MySubdomains2(N=n_domains, bounds=bounds2)
+
+        grid = Grid(extent=(Nx, Ny), shape=(Nx, Ny), subdomains=(my_sd1, my_sd2))
+
+        f = Function(name='f', grid=grid, dtype=np.int32)
+        g = Function(name='g', grid=grid, dtype=np.int32)
+
+        eq1 = Eq(f, f+2, subdomain=grid.subdomains['mydomains1'])
+        eq2 = Eq(g, g+2, subdomain=grid.subdomains['mydomains2'])
+        eq3 = Eq(f, f-1, subdomain=grid.subdomains['mydomains1'])
+        eq4 = Eq(g, g+1, subdomain=grid.subdomains['mydomains2'])
+
+        op = Operator([eq1, eq2, eq3, eq4])
+        op.apply()
+
+        expected = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 1, 1, 1, 0, 0, 3, 3, 3, 0],
+                             [0, 1, 1, 1, 0, 0, 3, 3, 3, 0],
+                             [0, 1, 1, 1, 0, 0, 3, 3, 3, 0],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                             [0, 1, 1, 1, 0, 0, 3, 3, 3, 0],
+                             [0, 1, 1, 1, 0, 0, 3, 3, 3, 0],
+                             [0, 1, 1, 1, 0, 0, 3, 3, 3, 0],
+                             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.int32)
+
+        assert((np.array(f.data[:]+g.data[:]) == expected).all())