diff --git a/.github/workflows/pytest-core-mpi.yml b/.github/workflows/pytest-core-mpi.yml index 05100d4d2d7..b0fb73abc26 100644 --- a/.github/workflows/pytest-core-mpi.yml +++ b/.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 diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index dddb2274e79..38a3571b17f 100644 --- a/devito/operations/interpolators.py +++ b/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 diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 98e0e3a2dd3..921c115a22a 100644 --- a/devito/types/sparse.py +++ b/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') diff --git a/docker/entrypoint.sh b/docker/entrypoint.sh index 2a3fab2bbd2..077a84cd05e 100644 --- a/docker/entrypoint.sh +++ b/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 diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py index 8a341f8a5ce..aced327cdca 100644 --- a/examples/seismic/acoustic/acoustic_example.py +++ b/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)]) diff --git a/examples/seismic/utils.py b/examples/seismic/utils.py index ef517d7ce75..42ca0ced721 100644 --- a/examples/seismic/utils.py +++ b/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)