From c3abf22435dc70b7660c6edfdb1fa8dcbf77c2be Mon Sep 17 00:00:00 2001
From: Andrew <AndrewCheng827827@gmail.com>
Date: Fri, 22 Sep 2023 11:42:55 +0100
Subject: [PATCH] compiler: Fix issue 2194

---
 devito/ir/support/basic.py | 13 +++++++++-
 tests/test_operator.py     | 52 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 64 insertions(+), 1 deletion(-)

diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py
index 98ba9da51a..0ca1da8796 100644
--- a/devito/ir/support/basic.py
+++ b/devito/ir/support/basic.py
@@ -426,7 +426,18 @@ def distance(self, other):
             else:
                 v = i - j
                 if v.is_Number and v.is_finite:
-                    return Vector(S.ImaginaryUnit)
+                    if i.is_Number and j.is_Number:
+                        return Vector(S.ImaginaryUnit)
+                    else:
+                        # For example:
+                        # self=W<u,[0,y]> and other=R<u,[0,y+1]>
+                        ret.append(v)
+
+                # Writing (reading) over an entire dimension, reading (writing)
+                # from one point. For example:
+                # self=R<u,[1,2]> and other=W<u,[1, y+1]>
+                elif (not i.is_Number or not j.is_Number):
+                    ret.append(S.Infinity)
 
         return Vector(*ret)
 
diff --git a/tests/test_operator.py b/tests/test_operator.py
index 3064565e3c..b2188d007c 100644
--- a/tests/test_operator.py
+++ b/tests/test_operator.py
@@ -1936,6 +1936,58 @@ def test_topofuse_w_numeric_dim(self):
 
         assert_structure(op, ['r,i', 'r'], 'r,i')
 
+    @pytest.mark.parametrize('eqns, expected, exp_trees, exp_iters', [
+        (['Eq(u[0, x], 1)',
+            'Eq(u[1, x], u[0, x + h_x] + u[0, x - h_x] - 2*u[0, x])'],
+            np.array([[1., 1., 1.], [-1., 0., -1.]]),
+            ['x', 'x'], 'x,x')
+    ])
+    def test_2194(self, eqns, expected, exp_trees, exp_iters):
+        grid = Grid(shape=(3, ))
+        u = TimeFunction(name='u', grid=grid)
+        x = grid.dimensions[0]
+        h_x = x.spacing  # noqa: F841
+
+        for i, e in enumerate(list(eqns)):
+            eqns[i] = eval(e)
+
+        op = Operator(eqns)
+        assert_structure(op, exp_trees, exp_iters)
+
+        op.apply()
+        assert(np.all(u.data[:] == expected[:]))
+
+    @pytest.mark.parametrize('eqns, expected, exp_trees, exp_iters', [
+        (['Eq(u[0, y], 1)', 'Eq(u[1, y], u[0, y + 1])'],
+            np.array([[1., 1.], [1., 0.]]),
+            ['y', 'y'], 'y,y'),
+        (['Eq(u[0, y], 1)', 'Eq(u[1, y], u[0, 2])'],
+            np.array([[1., 1.], [0., 0.]]),
+            ['y', 'y'], 'y,y'),
+        (['Eq(u[0, y], 1)', 'Eq(u[1, y], u[0, 1])'],
+            np.array([[1., 1.], [1., 1.]]),
+            ['y', 'y'], 'y,y'),
+        (['Eq(u[0, y], 1)', 'Eq(u[1, y], u[0, y + 1])'],
+            np.array([[1., 1.], [1., 0.]]),
+            ['y', 'y'], 'y,y'),
+        (['Eq(u[0, 1], 1)', 'Eq(u[x, y], u[0, y])'],
+            np.array([[0., 1.], [0., 1.]]),
+            ['xy'], 'x,y')
+    ])
+    def test_2194_v2(self, eqns, expected, exp_trees, exp_iters):
+        grid = Grid(shape=(2, 2))
+        u = Function(name='u', grid=grid)
+        x, y = grid.dimensions
+
+        for i, e in enumerate(list(eqns)):
+            eqns[i] = eval(e)
+
+        op = Operator(eqns)
+        assert_structure(op, exp_trees, exp_iters)
+
+        op.apply()
+        assert(np.all(u.data[:] == expected[:]))
+
 
 class TestInternals(object):