compiler: Fix DDA with degenerating ModuloDimensions

Author: FabioLuporini

devito/ir/support/basic.py

@@ -382,32 +382,36 @@ def distance(self, other):
                     if disjoint_test(self[n], other[n], sai, sit):
                         return Vector(S.ImaginaryUnit)
 
+            # Compute the distance along the current IterationInterval
             if self.function._mem_shared:
                 # Special case: the distance between two regular, thread-shared
-                # objects fallbacks to zero, as any other value would be nonsensical
+                # objects fallbacks to zero, as any other value would be
+                # nonsensical
+                ret.append(S.Zero)
+            elif degenerating_dimensions(sai, oai):
+                # Special case: `sai` and `oai` may be different symbolic objects
+                # but they can be proved to systematically generate the same value
                 ret.append(S.Zero)
-
             elif sai and oai and sai._defines & sit.dim._defines:
-                # E.g., `self=R<f,[t + 1, x]>`, `self.itintervals=(time, x)`, `ai=t`
+                # E.g., `self=R<f,[t + 1, x]>`, `self.itintervals=(time, x)`,
+                # and `ai=t`
                 if sit.direction is Backward:
                     ret.append(other[n] - self[n])
                 else:
                     ret.append(self[n] - other[n])
-
             elif not sai and not oai:
                 # E.g., `self=R<a,[3]>` and `other=W<a,[4]>`
                 if self[n] - other[n] == 0:
                     ret.append(S.Zero)
                 else:
                     break
-
             elif sai in self.ispace and oai in other.ispace:
                 # E.g., `self=R<f,[x, y]>`, `sai=time`,
                 #       `self.itintervals=(time, x, y)`, `n=0`
                 continue
-
             else:
-                # E.g., `self=R<u,[t+1, ii_src_0+1, ii_src_1+2]>`, `fi=p_src`, `n=1`
+                # E.g., `self=R<u,[t+1, ii_src_0+1, ii_src_1+2]>`, `fi=p_src`,
+                # and `n=1`
                 return vinf(ret)
 
         n = len(ret)
@@ -1408,3 +1412,19 @@ def disjoint_test(e0, e1, d, it):
     i1 = sympy.Interval(min(p10, p11), max(p10, p11))
 
     return not bool(i0.intersect(i1))
+
+
+def degenerating_dimensions(d0, d1):
+    """
+    True if `d0` and `d1` are Dimensions that are possibly symbolically
+    different, but they can be proved to systematically degenerate to the
+    same value, False otherwise.
+    """
+    # Case 1: ModuloDimensions of size 1
+    try:
+        if d0.is_Modulo and d1.is_Modulo and d0.modulo == d1.modulo == 1:
+            return True
+    except AttributeError:
+        pass
+
+    return False

devito/passes/iet/mpi.py

@@ -386,7 +386,9 @@ def _rule0(dep, hs, loc_indices):
 
 def _rule1(dep, hs, loc_indices):
     # E.g., `dep=W<f,[t1, x+1]> -> R<f,[t1, xl+1]>` and `loc_indices={t: t0}` => True
-    return any(dep.distance_mapper[d] == 0 and dep.source[d] is not v
+    return any(dep.distance_mapper[d] == 0 and
+               dep.source[d] is not v and
+               dep.sink[d] is not v
                for d, v in loc_indices.items())
 
 

tests/test_mpi.py

@@ -1832,6 +1832,31 @@ def test_haloupdate_buffer1(self, mode):
         # ModuloDimensions
         assert len([i for i in FindSymbols('dimensions').visit(op) if i.is_Modulo]) == 0
 
+    @pytest.mark.parallel(mode=1)
+    def test_haloupdate_buffer1_v2(self, mode):
+        grid = Grid((65, 65, 65), topology=('*', 1, '*'))
+
+        v1 = TimeFunction(name='v1', grid=grid, space_order=2, time_order=1,
+                          save=Buffer(1))
+        v2 = TimeFunction(name='v2', grid=grid, space_order=2, time_order=1,
+                          save=Buffer(1))
+
+        rec = SparseTimeFunction(name='rec', grid=grid, nt=500, npoint=65)
+
+        eqns = [Eq(v1.forward, v2.dx2 + v1),
+                Eq(v2.forward, v1.forward.dx2 + v2),
+                rec.interpolate(v2)]
+
+        op = Operator(eqns)
+        op.cfunction
+
+        # Ensure there's a halo exchange over v2 before the rec interpolation
+        section1 = op.body.body[-1].body[1].nodes[1]
+        assert section1.is_Section
+        calls = FindNodes(HaloUpdateCall).visit(section1)
+        assert len(calls) == 1
+        assert calls[0].arguments[0] is v2
+
 
 class TestOperatorAdvanced: