api: fix sparse position index for negative location

.github/workflows/pytest-core-mpi.yml

@@ -51,7 +51,8 @@ jobs:
     - name: Test examples
       run: |
         python3 scripts/clear_devito_cache.py
-        python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml examples/seismic
+        DEVITO_MPI=1 mpirun -n 2 python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml examples/seismic/acoustic
+        DEVITO_MPI=1 mpirun -n 2 python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml examples/seismic/tti
 
     - name: Upload coverage to Codecov
       uses: codecov/codecov-action@v4
@@ -87,4 +88,5 @@ jobs:
 
       - name: Test examples with MPI
         run: |
-          docker run --rm --env DEVITO_MPI=1 --name examplerun mpiexec -n 2 pytest examples/seismic
+          docker run --rm -e DEVITO_MPI=1 -e OMP_NUM_THREADS=1 --name examplerun devito_img mpiexec -n 2 pytest examples/seismic/acoustic
+          docker run --rm -e DEVITO_MPI=1 -e OMP_NUM_THREADS=1 --name examplerun devito_img mpiexec -n 2 pytest examples/seismic/tti

devito/operations/interpolators.py

@@ -23,7 +23,7 @@ def wrapper(interp, *args, **kwargs):
         funcs = set().union(*[retrieve_functions(a) for a in args])
         so = min({f.space_order for f in funcs if not f.is_SparseFunction} or {r})
         if so < r:
-            raise ValueError("Space order %d smaller than interpolation r %d" % (so, r))
+            raise ValueError("Space order %d too small for interpolation r %d" % (so, r))
         return func(interp, *args, **kwargs)
     return wrapper
 
@@ -216,7 +216,7 @@ def _positions(self, implicit_dims):
         return [Eq(v, INT(floor(k)), implicit_dims=implicit_dims)
                 for k, v in self.sfunction._position_map.items()]
 
-    def _interp_idx(self, variables, implicit_dims=None):
+    def _interp_idx(self, variables, implicit_dims=None, pos_only=None):
         """
         Generate interpolation indices for the DiscreteFunctions in ``variables``.
         """
@@ -235,6 +235,14 @@ def _interp_idx(self, variables, implicit_dims=None):
                     for ((k, c), p) in zip(mapper.items(), pos)})
                     for v in variables}
 
+        # Position only replacement, not radius dependent.
+        # E.g src.inject(vp(x)*src) needs to use vp[posx] at all points
+        # not vp[posx + rx]
+        if pos_only is not None:
+            idx_subs.update({v: v.subs({k: p
+                             for (k, p) in zip(mapper, pos)})
+                             for v in pos_only})
+
         return idx_subs, temps
 
     @check_radius
@@ -359,10 +367,9 @@ def _inject(self, field, expr, implicit_dims=None):
         # summing temp that wouldn't allow collapsing
         implicit_dims = implicit_dims + tuple(r.parent for r in self._rdim)
 
-        variables = variables + list(fields)
-
         # List of indirection indices for all adjacent grid points
-        idx_subs, temps = self._interp_idx(variables, implicit_dims=implicit_dims)
+        idx_subs, temps = self._interp_idx(list(fields), implicit_dims=implicit_dims,
+                                           pos_only=variables)
 
         # Substitute coordinate base symbols into the interpolation coefficients
         eqns = [Inc(_field.xreplace(idx_subs),
@@ -424,7 +431,7 @@ def _positions(self, implicit_dims):
         if self.sfunction.gridpoints is None:
             return super()._positions(implicit_dims)
         # No position temp as we have directly the gridpoints
-        return [Eq(p, k, implicit_dims=implicit_dims)
+        return [Eq(p, floor(k), implicit_dims=implicit_dims)
                 for (k, p) in self.sfunction._position_map.items()]
 
     @property
@@ -484,13 +491,7 @@ def _arg_defaults(self, coords=None, sfunc=None):
         if coords is None or sfunc is None:
             raise ValueError("No coordinates or sparse function provided")
         # Coords to indices
-        if sfunc.grid._distributor.nprocs == 1:
-            origin = sfunc.grid.origin
-        else:
-            # Already shifted to zero through scatter
-            origin = tuple([0]*sfunc.grid.dim)
-
-        coords = (coords - np.array(origin)) / np.array(sfunc.grid.spacing)
+        coords = coords / np.array(sfunc.grid.spacing)
         coords = coords - np.floor(coords)
 
         # Precompute sinc

devito/types/sparse.py

@@ -811,15 +811,17 @@ class SparseFunction(AbstractSparseFunction):
     __rkwargs__ = AbstractSparseFunction.__rkwargs__ + ('coordinates', 'interpolation')
 
     def __init_finalize__(self, *args, **kwargs):
-        super().__init_finalize__(*args, **kwargs)
-
-        interp = kwargs.get('interpolation', 'linear')
+        # Interpolation method
+        interp = kwargs.pop('interpolation', 'linear')
         self.interpolation = interp
         self.interpolator = _interpolators[interp](self)
-        self._radius = kwargs.get('r', _default_radius[interp])
+        self._radius = kwargs.pop('r', _default_radius[interp])
         if interp == 'sinc' and self._radius < 2:
             raise ValueError("The 'sinc' interpolator requires a radius of at least 2")
 
+        # Set space ordert to `r` for safety
+        super().__init_finalize__(*args, **kwargs)
+
         # Set up sparse point coordinates
         coordinates = kwargs.get('coordinates', kwargs.get('coordinates_data'))
         self._coordinates = self.__subfunc_setup__(coordinates, 'coords')

docker/entrypoint.sh

@@ -10,7 +10,7 @@ fi
 
 if [[ "$DEVITO_ARCH" = "icx" || "$DEVITO_ARCH" = "icc" ]]; then
     echo "Initializing oneapi environement"
-    source /opt/intel/oneapi/setvars.sh
+    source /opt/intel/oneapi/setvars.sh intel64
 fi
 
 if [[ -z "${DEPLOY_ENV}" ]]; then

examples/seismic/acoustic/acoustic_example.py

@@ -60,7 +60,7 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
     solver.jacobian(dm, autotune=autotune)
     info("Applying Gradient")
     solver.jacobian_adjoint(rec, u, autotune=autotune, checkpointing=checkpointing)
-    return summary.gflopss, summary.oi, summary.timings, [rec, u.data]
+    return summary.gflopss, summary.oi, summary.timings, [rec, u.data._local]
 
 
 @pytest.mark.parametrize('shape', [(101,), (51, 51), (16, 16, 16)])

examples/seismic/utils.py

@@ -13,9 +13,11 @@
 def setup_geometry(model, tn, f0=0.010, interpolation='linear', **kwargs):
     # Source and receiver geometries
     src_coordinates = np.empty((1, model.dim))
-    src_coordinates[0, :] = np.array(model.domain_size) * .5
     if model.dim > 1:
+        src_coordinates[0, :] = np.array(model.domain_size) * .5
         src_coordinates[0, -1] = model.origin[-1] + model.spacing[-1]
+    else:
+        src_coordinates[0, 0] = 2 * model.spacing[0]
 
     rec_coordinates = setup_rec_coords(model)