Repository: devitocodes/devito
Issue: MPI with SparseFunction on the mpi domain boundary

There seem to be an issue when the source is on or near the MPI halo in some cases.

MFE  to be made.

Changed files:
devito/types/sparse.py, tests/test_caching.py, tests/test_interpolation.py, devito/ir/equations/algorithms.py, devito/operations/interpolators.py

diff --git a/devito/ir/equations/algorithms.py b/devito/ir/equations/algorithms.py
index 0f8bcceff..90ba26282 100644
--- a/devito/ir/equations/algorithms.py
+++ b/devito/ir/equations/algorithms.py
@@ -1,5 +1,4 @@
 from collections.abc import Iterable
-
 from operator import attrgetter
 
 from devito.symbolics import (retrieve_functions, retrieve_indexed, split_affine,
diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py
index dc013cb8e..f4c079075 100644
--- a/devito/operations/interpolators.py
+++ b/devito/operations/interpolators.py
@@ -187,11 +187,16 @@ class LinearInterpolator(GenericInterpolator):
             # Track Indexed substitutions
             idx_subs.append(mapper)
 
-        # Temporaries for the indirection dimensions
+        # Temporaries for the position
         temps = [Eq(v, k, implicit_dims=self.sfunction.dimensions)
-                 for k, v in points.items()]
+                 for k, v in self.sfunction._position_map.items()]
+        # Temporaries for the indirection dimensions
+        temps.extend([Eq(v, k.subs(self.sfunction._position_map),
+                         implicit_dims=self.sfunction.dimensions)
+                      for k, v in points.items()])
         # Temporaries for the coefficients
-        temps.extend([Eq(p, c, implicit_dims=self.sfunction.dimensions)
+        temps.extend([Eq(p, c.subs(self.sfunction._position_map),
+                         implicit_dims=self.sfunction.dimensions)
                       for p, c in zip(self.sfunction._point_symbols,
                                       self.sfunction._coordinate_bases(field_offset))])
 
diff --git a/devito/types/sparse.py b/devito/types/sparse.py
index c4b1eb9e3..23548eb38 100644
--- a/devito/types/sparse.py
+++ b/devito/types/sparse.py
@@ -8,7 +8,8 @@ from cached_property import cached_property
 from devito.finite_differences import generate_fd_shortcuts
 from devito.mpi import MPI, SparseDistributor
 from devito.operations import LinearInterpolator, PrecomputedInterpolator
-from devito.symbolics import INT, cast_mapper, indexify, retrieve_function_carriers
+from devito.symbolics import (INT, FLOOR, cast_mapper, indexify,
+                              retrieve_function_carriers)
 from devito.tools import (ReducerMap, as_tuple, flatten, prod, filter_ordered,
                           memoized_meth)
 from devito.types.dense import DiscreteFunction, Function, SubFunction
@@ -535,6 +536,27 @@ class SparseFunction(AbstractSparseFunction):
         return tuple(Scalar(name='p%s' % d, dtype=self.dtype)
                      for d in self.grid.dimensions)
 
+    @cached_property
+    def _position_map(self):
+        """
+        Symbols map for the position of the sparse points relative to the grid
+        origin.
+
+        Notes
+        -----
+        The expression `(coord - origin)/spacing` could also be computed in the
+        mathematically equivalent expanded form `coord/spacing -
+        origin/spacing`. This particular form is problematic when a sparse
+        point is in close proximity of the grid origin, since due to a larger
+        machine precision error it may cause a +-1 error in the computation of
+        the position. We mitigate this problem by computing the positions
+        individually (hence the need for a position map).
+        """
+        symbols = [Scalar(name='pos%s' % d, dtype=self.dtype)
+                   for d in self.grid.dimensions]
+        return OrderedDict([(c - o, p) for p, c, o in
+                            zip(symbols, self._coordinate_symbols, self.grid.origin)])
+
     @cached_property
     def _point_increments(self):
         """Index increments in each dimension for each point symbol."""
@@ -550,19 +572,18 @@ class SparseFunction(AbstractSparseFunction):
     @cached_property
     def _coordinate_indices(self):
         """Symbol for each grid index according to the coordinates."""
-        indices = self.grid.dimensions
-        return tuple([INT(sympy.Function('floor')((c - o) / i.spacing))
-                      for c, o, i in zip(self._coordinate_symbols, self.grid.origin,
-                                         indices[:self.grid.dim])])
+        return tuple([INT(FLOOR((c - o) / i.spacing))
+                      for c, o, i in zip(self._coordinate_symbols,
+                                         self.grid.origin,
+                                         self.grid.dimensions[:self.grid.dim])])
 
     def _coordinate_bases(self, field_offset):
         """Symbol for the base coordinates of the reference grid point."""
-        indices = self.grid.dimensions
         return tuple([cast_mapper[self.dtype](c - o - idx * i.spacing)
                       for c, o, idx, i, of in zip(self._coordinate_symbols,
                                                   self.grid.origin,
                                                   self._coordinate_indices,
-                                                  indices[:self.grid.dim],
+                                                  self.grid.dimensions[:self.grid.dim],
                                                   field_offset)])
 
     @memoized_meth
@@ -629,9 +650,13 @@ class SparseFunction(AbstractSparseFunction):
                          if f.is_SparseFunction}
             out = indexify(expr).xreplace({f._sparse_dim: cd for f in functions})
 
-        # Temporaries for the indirection dimensions
+        # Temporaries for the position
         temps = [Eq(v, k, implicit_dims=self.dimensions)
-                 for k, v in points.items() if v in conditions]
+                 for k, v in self._position_map.items()]
+        # Temporaries for the indirection dimensions
+        temps.extend([Eq(v, k.subs(self._position_map),
+                         implicit_dims=self.dimensions)
+                      for k, v in points.items() if v in conditions])
 
         return out, temps
 
diff --git a/tests/test_caching.py b/tests/test_caching.py
index 446fc304a..ac5b7293d 100644
--- a/tests/test_caching.py
+++ b/tests/test_caching.py
@@ -523,8 +523,8 @@ class TestCaching(object):
         i = u.inject(expr=u, field=u)
 
         # created: ii_u_0*2 (Symbol and ConditionalDimension), ii_u_1*2, ii_u_2*2,
-        # ii_u_3*2, px, py, u_coords (as indexified),
-        ncreated = 2+2+2+2+1+1+1
+        # ii_u_3*2, px, py, posx, posy, u_coords (as indexified),
+        ncreated = 2+2+2+2+2+1+1+1
         # Note that injection is now lazy so no new symbols should be created
         assert len(_SymbolCache) == cur_cache_size
         i.evaluate
diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py
index aacc94469..cb966683d 100644
--- a/tests/test_interpolation.py
+++ b/tests/test_interpolation.py
@@ -472,3 +472,24 @@ def test_position(shape):
                                  o_x=ox_g, o_y=oy_g, o_z=oz_g)
 
     assert(np.allclose(rec.data, rec1.data, atol=1e-5))
+
+
+def test_edge_sparse():
+    """
+    Test that interpolation uses the correct point for the edge case
+    where the sparse point is at the origin with non rational grid spacing.
+    Due to round up error the interpolation would use the halo point instead of
+    the point (0, 0) without the factorizaion of the expressions.
+    """
+    grid = Grid(shape=(16, 16), extent=(225., 225.), origin=(25., 35.))
+    u = unit_box(shape=(16, 16), grid=grid)
+    u._data_with_outhalo[:u.space_order, :] = -1
+    u._data_with_outhalo[:, :u.space_order] = -1
+    sf1 = SparseFunction(name='s', grid=u.grid, npoint=1)
+    sf1.coordinates.data[0, :] = (25.0, 35.0)
+
+    expr = sf1.interpolate(u)
+    subs = {d.spacing: v for d, v in zip(u.grid.dimensions, u.grid.spacing)}
+    op = Operator(expr, subs=subs)
+    op()
+    assert sf1.data[0] == 0