Merge pull request #2557 from devitocodes/coordinates

dsl: ensure SparseFunction coordinates and point symbols are always real

Author: mloubout

devito/operations/interpolators.py

@@ -472,7 +472,8 @@ def _weights(self, subdomain=None):
     @cached_property
     def _point_symbols(self):
         """Symbol for coordinate value in each Dimension of the point."""
-        return DimensionTuple(*(Symbol(name='p%s' % d, dtype=self.sfunction.dtype)
+        dtype = self.sfunction.coordinates.dtype
+        return DimensionTuple(*(Symbol(name=f'p{d}', dtype=dtype)
                                 for d in self.grid.dimensions),
                               getters=self.grid.dimensions)
 

devito/types/sparse.py

@@ -202,6 +202,10 @@ def __subfunc_setup__(self, suffix, keys, dtype=None, inkwargs=False, **kwargs):
             else:
                 dtype = dtype or self.dtype
 
+        # Complex coordinates are not valid, so fall back to corresponding
+        # real floating point type if dtype is complex.
+        dtype = dtype(0).real.__class__
+
         sf = SparseSubFunction(
             name=name, dtype=dtype, dimensions=dimensions,
             shape=shape, space_order=0, initializer=key, alias=self.alias,

tests/test_interpolation.py

@@ -841,6 +841,18 @@ def test_sinc_accuracy(r, tol):
     assert err_lin > 0.01
 
 
+@pytest.mark.parametrize('dtype, expected', [(np.complex64, np.float32),
+                                             (np.complex128, np.float64)])
+def test_point_symbol_types(dtype, expected):
+    """Test that positions are always real"""
+    grid = Grid(shape=(11,))
+    s = SparseFunction(name='src', npoint=1,
+                       grid=grid, dtype=dtype)
+    point_symbol = s.interpolator._point_symbols[0]
+
+    assert point_symbol.dtype is expected
+
+
 class SD0(SubDomain):
     name = 'sd0'
 

tests/test_sparse.py

@@ -497,6 +497,20 @@ def test_mpi_no_data(self, mode):
         ftest.data[:] = expected
         assert np.all(m.data[0, :, :] == ftest.data[:])
 
+    @pytest.mark.parametrize('dtype, expected', [(np.complex64, np.float32),
+                                                 (np.complex128, np.float64),
+                                                 (np.float16, np.float16)])
+    def test_coordinate_type(self, dtype, expected):
+        """
+        Test that coordinates are always real and SparseFunction dtype is
+        otherwise preserved.
+        """
+        grid = Grid(shape=(11,))
+        s = SparseFunction(name='src', npoint=1,
+                           grid=grid, dtype=dtype)
+
+        assert s.coordinates.dtype is expected
+
 
 if __name__ == "__main__":
     TestMatrixSparseTimeFunction().test_mpi_no_data()