From cf4b65ee84721c41540752a03707b31c44d0a297 Mon Sep 17 00:00:00 2001 From: mloubout Date: Thu, 28 Mar 2024 14:23:58 -0400 Subject: [PATCH 1/8] api: support sinc interpolation --- devito/operations/interpolators.py | 65 ++++++++++++++++++- devito/types/sparse.py | 24 +++++-- examples/seismic/acoustic/acoustic_example.py | 2 +- examples/seismic/source.py | 2 + tests/test_interpolation.py | 6 +- 5 files changed, 91 insertions(+), 8 deletions(-) diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index 98e0105e8a..ced9fda21a 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -5,11 +5,15 @@ from cached_property import cached_property from devito.finite_differences.differentiable import Mul +<<<<<<< HEAD from devito.finite_differences.elementary import floor +======= +from devito.finite_differences.elementary import floor, sqrt +>>>>>>> 5815a494d (api: support sinc interpolation) from devito.symbolics import retrieve_function_carriers, retrieve_functions, INT from devito.tools import as_tuple, flatten, filter_ordered from devito.types import (ConditionalDimension, Eq, Inc, Evaluable, Symbol, - CustomDimension) + CustomDimension, SubFunction) from devito.types.utils import DimensionTuple __all__ = ['LinearInterpolator', 'PrecomputedInterpolator'] @@ -133,6 +137,9 @@ def inject(self, *args, **kwargs): def interpolate(self, *args, **kwargs): pass + def _arg_defaults(self, **args): + return {} + class WeightedInterpolator(GenericInterpolator): @@ -421,3 +428,59 @@ def _weights(self): for (ri, rd) in enumerate(self._rdim)] return Mul(*[self.interpolation_coeffs.subs(mapper) for mapper in mappers]) + + +class SincInterpolator(PrecomputedInterpolator): + """ + Hicks windowed sinc interpolation scheme. + + Arbitrary source and receiver positioning in finite‐difference schemes + using Kaiser windowed sinc functions + + https://library.seg.org/doi/10.1190/1.1451454 + + """ + + # Table 1 + _b_table = {1: 0.0, 2: 1.84, 3: 3.04, + 4: 4.14, 5: 5.26, 6: 6.40, + 7: 7.51, 8: 8.56, 9: 9.56, 10: 10.64} + + @property + def interpolation_coeffs(self): + coeffs = {} + shape = (self.sfunction.npoint, 2 * self.r) + for r in self._rdim: + dimensions = (self.sfunction._sparse_dim, r.parent) + sf = SubFunction(name="wsinc%s" % r.name, dtype=self.sfunction.dtype, + shape=shape, dimensions=dimensions, + space_order=0, alias=self.sfunction.alias, + distributor=self.sfunction._distributor, + parent=self.sfunction) + coeffs[r] = sf + return coeffs + + @property + def _weights(self): + return Mul(*[w._subs(rd, rd-rd.parent.symbolic_min) + for (rd, w) in self.interpolation_coeffs.items()]) + + def _arg_defaults(self, coords=None, sfunc=None): + args = {} + b = self._b_table[self.r] + b0 = np.i0(b) + if coords is None or sfunc is None: + raise ValueError("No coordinates or sparse function provided") + # Coords to indices + coords = (coords - np.array(sfunc.grid.origin)) / np.array(sfunc.grid.spacing) + coords = np.floor(coords) - coords + 1 - self.r + # Precompute sinc + for j, r in enumerate(self._rdim): + data = np.zeros((coords.shape[0], 2*self.r), dtype=sfunc.dtype) + for ri in range(2*self.r): + rpos = ri + coords[:, j] + num = np.i0(b*sqrt(1 - (rpos/self.r)**2)) + data[:, ri] = num / b0 * np.sinc(rpos) + args[self.interpolation_coeffs[r].name] = data + + return args diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 6829698582..602c137c7f 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -12,7 +12,8 @@ from devito.finite_differences import generate_fd_shortcuts from devito.mpi import MPI, SparseDistributor -from devito.operations import LinearInterpolator, PrecomputedInterpolator +from devito.operations import (LinearInterpolator, PrecomputedInterpolator, + SincInterpolator) from devito.symbolics import indexify, retrieve_function_carriers from devito.tools import (ReducerMap, as_tuple, flatten, prod, filter_ordered, is_integer, dtype_to_mpidtype) @@ -29,6 +30,10 @@ 'PrecomputedSparseTimeFunction', 'MatrixSparseTimeFunction'] +_interpolators = {'linear': LinearInterpolator, 'sinc': SincInterpolator} +_default_radius = {'linear': 1, 'sinc': 2} + + class AbstractSparseFunction(DiscreteFunction): """ @@ -799,16 +804,18 @@ class SparseFunction(AbstractSparseFunction): is_SparseFunction = True - _radius = 1 """The radius of the stencil operators provided by the SparseFunction.""" _sub_functions = ('coordinates',) - __rkwargs__ = AbstractSparseFunction.__rkwargs__ + ('coordinates',) + __rkwargs__ = AbstractSparseFunction.__rkwargs__ + ('coordinates', 'interpolator') def __init_finalize__(self, *args, **kwargs): super().__init_finalize__(*args, **kwargs) - self.interpolator = LinearInterpolator(self) + + interp = kwargs.get('interpolator', 'linear') + self.interpolator = _interpolators[interp](self) + self._radius = kwargs.get('r', _default_radius[interp]) # Set up sparse point coordinates coordinates = kwargs.get('coordinates', kwargs.get('coordinates_data')) @@ -827,6 +834,15 @@ def _decomposition(self): mapper = {self._sparse_dim: self._distributor.decomposition[self._sparse_dim]} return tuple(mapper.get(d) for d in self.dimensions) + def _arg_defaults(self, alias=None): + defaults = super()._arg_defaults(alias=alias) + + key = alias or self + coords = + defaults.update(key.interpolator._arg_defaults(coords=coords, + sfunc=key)) + return defaults + class SparseTimeFunction(AbstractSparseTimeFunction, SparseFunction): """ diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py index 7e922cbb6e..f2376105e0 100644 --- a/examples/seismic/acoustic/acoustic_example.py +++ b/examples/seismic/acoustic/acoustic_example.py @@ -4,7 +4,7 @@ from devito.logger import info from devito import Constant, Function, smooth, norm from examples.seismic.acoustic import AcousticWaveSolver -from examples.seismic import demo_model, setup_geometry, seismic_args +from examples.seismic import demo_model, setup_geometry, seismic_args, Receiver def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), diff --git a/examples/seismic/source.py b/examples/seismic/source.py index a736c5bb4f..d0cdef2613 100644 --- a/examples/seismic/source.py +++ b/examples/seismic/source.py @@ -117,6 +117,8 @@ def __args_setup__(cls, *args, **kwargs): raise TypeError("Need either `npoint` or `coordinates`") kwargs['npoint'] = coordinates.shape[0] + kwargs.setdefault('r', 1) + kwargs.setdefault('interpolator', 'linear') return args, kwargs def __init_finalize__(self, *args, **kwargs): diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 506a3b88de..ae013572c8 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -752,12 +752,14 @@ def test_inject_function(): assert u.data[1, i, j] == 0 -def test_interpolation_radius(): +@pytest.mark.parametrize('r, interp', [(2, 'linear'), (4, 'cubic')]) +def test_interpolation_radius(rm interp): nt = 11 grid = Grid(shape=(5, 5)) u = TimeFunction(name="u", grid=grid, space_order=0) - src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1) + src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1, + r=r, interpolator=interp) try: src.interpolate(u) assert False From 83bbef339ec0cb343bfd4cadabc49e32c38bab45 Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 1 Apr 2024 12:56:00 -0400 Subject: [PATCH 2/8] api: numpify sinc precompute and add test --- devito/operations/interpolators.py | 33 +++++--- devito/types/sparse.py | 11 ++- examples/seismic/acoustic/acoustic_example.py | 21 +++-- examples/seismic/acoustic/operators.py | 24 ++---- examples/seismic/source.py | 2 - examples/seismic/utils.py | 31 ++++++-- tests/test_interpolation.py | 79 ++++++++++++++++++- 7 files changed, 150 insertions(+), 51 deletions(-) diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index ced9fda21a..6871b69da3 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -2,21 +2,18 @@ from functools import wraps import sympy +import numpy as np from cached_property import cached_property from devito.finite_differences.differentiable import Mul -<<<<<<< HEAD from devito.finite_differences.elementary import floor -======= -from devito.finite_differences.elementary import floor, sqrt ->>>>>>> 5815a494d (api: support sinc interpolation) from devito.symbolics import retrieve_function_carriers, retrieve_functions, INT from devito.tools import as_tuple, flatten, filter_ordered from devito.types import (ConditionalDimension, Eq, Inc, Evaluable, Symbol, CustomDimension, SubFunction) from devito.types.utils import DimensionTuple -__all__ = ['LinearInterpolator', 'PrecomputedInterpolator'] +__all__ = ['LinearInterpolator', 'PrecomputedInterpolator', 'SincInterpolator'] def check_radius(func): @@ -129,6 +126,8 @@ class GenericInterpolator(ABC): Abstract base class defining the interface for an interpolator. """ + _name = "generic" + @abstractmethod def inject(self, *args, **kwargs): pass @@ -137,6 +136,10 @@ def inject(self, *args, **kwargs): def interpolate(self, *args, **kwargs): pass + @property + def name(self): + return self._name + def _arg_defaults(self, **args): return {} @@ -149,6 +152,8 @@ class WeightedInterpolator(GenericInterpolator): and multiplied at a given point: `w[x, y] = wx[x] * wy[y]` """ + _name = 'weighted' + def __init__(self, sfunction): self.sfunction = sfunction @@ -377,6 +382,9 @@ class LinearInterpolator(WeightedInterpolator): ---------- sfunction: The SparseFunction that this Interpolator operates on. """ + + _name = 'linear' + @property def _weights(self): c = [(1 - p) * (1 - r) + p * r @@ -410,6 +418,8 @@ class PrecomputedInterpolator(WeightedInterpolator): sfunction: The SparseFunction that this Interpolator operates on. """ + _name = 'precomp' + def _positions(self, implicit_dims): if self.sfunction.gridpoints is None: return super()._positions(implicit_dims) @@ -441,12 +451,14 @@ class SincInterpolator(PrecomputedInterpolator): """ + _name = 'sinc' + # Table 1 - _b_table = {1: 0.0, 2: 1.84, 3: 3.04, + _b_table = {2: 2.94, 3: 4.53, 4: 4.14, 5: 5.26, 6: 6.40, 7: 7.51, 8: 8.56, 9: 9.56, 10: 10.64} - @property + @cached_property def interpolation_coeffs(self): coeffs = {} shape = (self.sfunction.npoint, 2 * self.r) @@ -473,13 +485,14 @@ def _arg_defaults(self, coords=None, sfunc=None): raise ValueError("No coordinates or sparse function provided") # Coords to indices coords = (coords - np.array(sfunc.grid.origin)) / np.array(sfunc.grid.spacing) - coords = np.floor(coords) - coords + 1 - self.r + coords = coords - np.floor(coords) + # Precompute sinc for j, r in enumerate(self._rdim): data = np.zeros((coords.shape[0], 2*self.r), dtype=sfunc.dtype) for ri in range(2*self.r): - rpos = ri + coords[:, j] - num = np.i0(b*sqrt(1 - (rpos/self.r)**2)) + rpos = ri - self.r + 1 - coords[:, j] + num = np.i0(b*np.sqrt(1 - (rpos/self.r)**2)) data[:, ri] = num / b0 * np.sinc(rpos) args[self.interpolation_coeffs[r].name] = data diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 602c137c7f..98e0e3a2dd 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -31,7 +31,7 @@ _interpolators = {'linear': LinearInterpolator, 'sinc': SincInterpolator} -_default_radius = {'linear': 1, 'sinc': 2} +_default_radius = {'linear': 1, 'sinc': 4} class AbstractSparseFunction(DiscreteFunction): @@ -808,14 +808,17 @@ class SparseFunction(AbstractSparseFunction): _sub_functions = ('coordinates',) - __rkwargs__ = AbstractSparseFunction.__rkwargs__ + ('coordinates', 'interpolator') + __rkwargs__ = AbstractSparseFunction.__rkwargs__ + ('coordinates', 'interpolation') def __init_finalize__(self, *args, **kwargs): super().__init_finalize__(*args, **kwargs) - interp = kwargs.get('interpolator', 'linear') + interp = kwargs.get('interpolation', 'linear') + self.interpolation = interp self.interpolator = _interpolators[interp](self) self._radius = kwargs.get('r', _default_radius[interp]) + if interp == 'sinc' and self._radius < 2: + raise ValueError("The 'sinc' interpolator requires a radius of at least 2") # Set up sparse point coordinates coordinates = kwargs.get('coordinates', kwargs.get('coordinates_data')) @@ -838,7 +841,7 @@ def _arg_defaults(self, alias=None): defaults = super()._arg_defaults(alias=alias) key = alias or self - coords = + coords = defaults.get(key.coordinates.name, key.coordinates.data) defaults.update(key.interpolator._arg_defaults(coords=coords, sfunc=key)) return defaults diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py index f2376105e0..8a341f8a5c 100644 --- a/examples/seismic/acoustic/acoustic_example.py +++ b/examples/seismic/acoustic/acoustic_example.py @@ -4,7 +4,7 @@ from devito.logger import info from devito import Constant, Function, smooth, norm from examples.seismic.acoustic import AcousticWaveSolver -from examples.seismic import demo_model, setup_geometry, seismic_args, Receiver +from examples.seismic import demo_model, setup_geometry, seismic_args def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), @@ -15,7 +15,7 @@ def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), fs=fs, **kwargs) # Source and receiver geometries - geometry = setup_geometry(model, tn) + geometry = setup_geometry(model, tn, **kwargs) # Create solver object to provide relevant operators solver = AcousticWaveSolver(model, geometry, kernel=kernel, @@ -65,16 +65,21 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0, @pytest.mark.parametrize('shape', [(101,), (51, 51), (16, 16, 16)]) @pytest.mark.parametrize('k', ['OT2', 'OT4']) -def test_isoacoustic_stability(shape, k): +@pytest.mark.parametrize('interp', ['linear', 'sinc']) +def test_isoacoustic_stability(shape, k, interp): spacing = tuple([20]*len(shape)) - _, _, _, [rec, _] = run(shape=shape, spacing=spacing, tn=20000.0, nbl=0, kernel=k) + _, _, _, [rec, _] = run(shape=shape, spacing=spacing, tn=20000.0, nbl=0, kernel=k, + interpolation=interp) assert np.isfinite(norm(rec)) -@pytest.mark.parametrize('fs, normrec, dtype', [(True, 369.955, np.float32), - (False, 459.1678, np.float64)]) -def test_isoacoustic(fs, normrec, dtype): - _, _, _, [rec, _] = run(fs=fs, dtype=dtype) +@pytest.mark.parametrize('fs, normrec, dtype, interp', [ + (True, 369.955, np.float32, 'linear'), + (False, 459.1678, np.float64, 'linear'), + (True, 402.216, np.float32, 'sinc'), + (False, 509.0681, np.float64, 'sinc')]) +def test_isoacoustic(fs, normrec, dtype, interp): + _, _, _, [rec, _] = run(fs=fs, dtype=dtype, interpolation=interp) assert np.isclose(norm(rec), normrec, rtol=1e-3, atol=0) diff --git a/examples/seismic/acoustic/operators.py b/examples/seismic/acoustic/operators.py index 3623915a28..64e7f63284 100644 --- a/examples/seismic/acoustic/operators.py +++ b/examples/seismic/acoustic/operators.py @@ -1,6 +1,5 @@ from devito import Eq, Operator, Function, TimeFunction, Inc, solve, sign from devito.symbolics import retrieve_functions, INT -from examples.seismic import PointSource, Receiver def freesurface(model, eq): @@ -119,11 +118,8 @@ def ForwardOperator(model, geometry, space_order=4, u = TimeFunction(name='u', grid=model.grid, save=geometry.nt if save else None, time_order=2, space_order=space_order) - src = PointSource(name='src', grid=geometry.grid, time_range=geometry.time_axis, - npoint=geometry.nsrc) - - rec = Receiver(name='rec', grid=geometry.grid, time_range=geometry.time_axis, - npoint=geometry.nrec) + src = geometry.src + rec = geometry.rec s = model.grid.stepping_dim.spacing eqn = iso_stencil(u, model, kernel) @@ -160,10 +156,8 @@ def AdjointOperator(model, geometry, space_order=4, v = TimeFunction(name='v', grid=model.grid, save=None, time_order=2, space_order=space_order) - srca = PointSource(name='srca', grid=model.grid, time_range=geometry.time_axis, - npoint=geometry.nsrc) - rec = Receiver(name='rec', grid=model.grid, time_range=geometry.time_axis, - npoint=geometry.nrec) + srca = geometry.new_src(name='srca', src_type=None) + rec = geometry.rec s = model.grid.stepping_dim.spacing eqn = iso_stencil(v, model, kernel, forward=False) @@ -206,8 +200,7 @@ def GradientOperator(model, geometry, space_order=4, save=True, else None, time_order=2, space_order=space_order) v = TimeFunction(name='v', grid=model.grid, save=None, time_order=2, space_order=space_order) - rec = Receiver(name='rec', grid=model.grid, time_range=geometry.time_axis, - npoint=geometry.nrec) + rec = geometry.rec s = model.grid.stepping_dim.spacing eqn = iso_stencil(v, model, kernel, forward=False) @@ -244,11 +237,8 @@ def BornOperator(model, geometry, space_order=4, m = model.m # Create source and receiver symbols - src = Receiver(name='src', grid=model.grid, time_range=geometry.time_axis, - npoint=geometry.nsrc) - - rec = Receiver(name='rec', grid=model.grid, time_range=geometry.time_axis, - npoint=geometry.nrec) + src = geometry.src + rec = geometry.rec # Create wavefields and a dm field u = TimeFunction(name="u", grid=model.grid, save=None, diff --git a/examples/seismic/source.py b/examples/seismic/source.py index d0cdef2613..a736c5bb4f 100644 --- a/examples/seismic/source.py +++ b/examples/seismic/source.py @@ -117,8 +117,6 @@ def __args_setup__(cls, *args, **kwargs): raise TypeError("Need either `npoint` or `coordinates`") kwargs['npoint'] = coordinates.shape[0] - kwargs.setdefault('r', 1) - kwargs.setdefault('interpolator', 'linear') return args, kwargs def __init_finalize__(self, *args, **kwargs): diff --git a/examples/seismic/utils.py b/examples/seismic/utils.py index e6c0648a83..ef517d7ce7 100644 --- a/examples/seismic/utils.py +++ b/examples/seismic/utils.py @@ -3,13 +3,14 @@ from devito import error, configuration, warning from devito.tools import Pickable +from devito.types.sparse import _default_radius from .source import * __all__ = ['AcquisitionGeometry', 'setup_geometry', 'seismic_args'] -def setup_geometry(model, tn, f0=0.010): +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 @@ -18,8 +19,10 @@ def setup_geometry(model, tn, f0=0.010): rec_coordinates = setup_rec_coords(model) + r = kwargs.get('r', _default_radius[interpolation]) geometry = AcquisitionGeometry(model, rec_coordinates, src_coordinates, - t0=0.0, tn=tn, src_type='Ricker', f0=f0) + t0=0.0, tn=tn, src_type='Ricker', f0=f0, + interpolation=interpolation, r=r) return geometry @@ -58,7 +61,7 @@ class AcquisitionGeometry(Pickable): """ __rargs__ = ('grid', 'rec_positions', 'src_positions', 't0', 'tn') - __rkwargs__ = ('f0', 'src_type') + __rkwargs__ = ('f0', 'src_type', 'interpolation', 'r') def __init__(self, model, rec_positions, src_positions, t0, tn, **kwargs): """ @@ -85,6 +88,8 @@ def __init__(self, model, rec_positions, src_positions, t0, tn, **kwargs): self._dt = model.critical_dt self._t0 = t0 self._tn = tn + self._interpolation = kwargs.get('interpolation', 'linear') + self._r = kwargs.get('r', _default_radius[self.interpolation]) def resample(self, dt): self._dt = dt @@ -140,6 +145,14 @@ def nsrc(self): def dtype(self): return self.grid.dtype + @property + def r(self): + return self._r + + @property + def interpolation(self): + return self._interpolation + @property def rec(self): return self.new_rec() @@ -147,7 +160,8 @@ def rec(self): def new_rec(self, name='rec'): return Receiver(name=name, grid=self.grid, time_range=self.time_axis, npoint=self.nrec, - coordinates=self.rec_positions) + coordinates=self.rec_positions, + interpolation=self.interpolation, r=self._r) @property def adj_src(self): @@ -157,7 +171,8 @@ def adj_src(self): adj_src = sources[self.src_type](name='rec', grid=self.grid, f0=self.f0, time_range=self.time_axis, npoint=self.nrec, coordinates=self.rec_positions, - t0=self._t0w, a=self._a) + t0=self._t0w, a=self._a, + interpolation=self.interpolation, r=self._r) # Revert time axis to have a proper shot record and not compute on zeros for i in range(self.nrec): adj_src.data[:, i] = adj_src.wavelet[::-1] @@ -172,12 +187,14 @@ def new_src(self, name='src', src_type='self'): warning("No source type defined, returning uninitiallized (zero) source") return PointSource(name=name, grid=self.grid, time_range=self.time_axis, npoint=self.nsrc, - coordinates=self.src_positions) + coordinates=self.src_positions, + interpolation=self.interpolation, r=self._r) else: return sources[self.src_type](name=name, grid=self.grid, f0=self.f0, time_range=self.time_axis, npoint=self.nsrc, coordinates=self.src_positions, - t0=self._t0w, a=self._a) + t0=self._t0w, a=self._a, + interpolation=self.interpolation, r=self._r) sources = {'Wavelet': WaveletSource, 'Ricker': RickerSource, 'Gabor': GaborSource} diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index ae013572c8..664b882f4f 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -8,9 +8,10 @@ Function, TimeFunction, DefaultDimension, Eq, switchconfig, PrecomputedSparseFunction, PrecomputedSparseTimeFunction, MatrixSparseTimeFunction) +from devito.operations.interpolators import LinearInterpolator, SincInterpolator from examples.seismic import (demo_model, TimeAxis, RickerSource, Receiver, AcquisitionGeometry) -from examples.seismic.acoustic import AcousticWaveSolver +from examples.seismic.acoustic import AcousticWaveSolver, acoustic_setup import scipy.sparse @@ -753,15 +754,87 @@ def test_inject_function(): @pytest.mark.parametrize('r, interp', [(2, 'linear'), (4, 'cubic')]) -def test_interpolation_radius(rm interp): +def test_interpolation_radius(r, interp): nt = 11 grid = Grid(shape=(5, 5)) u = TimeFunction(name="u", grid=grid, space_order=0) src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1, - r=r, interpolator=interp) + r=r, interpolation=interp) try: src.interpolate(u) assert False except ValueError: assert True + + +def test_interp_default(): + nt = 3 + grid = Grid(shape=(5, 5)) + + src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1) + assert isinstance(src.interpolator, LinearInterpolator) + assert src.r == 1 + + src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1, interpolation='sinc') + assert isinstance(src.interpolator, SincInterpolator) + assert src.r == 2 + + src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1, + interpolation='sinc', r=6) + assert isinstance(src.interpolator, SincInterpolator) + assert src.r == 6 + + +@pytest.mark.parametrize('r, tol', [(2, 0.051), (3, 0.003), (4, 0.008), + (5, 0.002), (6, 0.0005), (7, 8e-5), + (8, 6e-5), (9, 5e-5), (10, 4.2e-5)]) +def test_sinc_accuracy(r, tol): + so = max(2, r) + solver_lin = acoustic_setup(preset='constant-isotropic', shape=(101, 101), + spacing=(10, 10), interpolation='linear', space_order=so) + solver_sinc = acoustic_setup(preset='constant-isotropic', shape=(101, 101), + spacing=(10, 10), interpolation='sinc', r=r, + space_order=so) + + # On node source + s_node = [500, 500] + src_n = solver_lin.geometry.src + src_n.coordinates.data[:] = s_node + + # Half node src + s_mid = [505, 505] + src_h = solver_lin.geometry.src + src_h.coordinates.data[:] = s_mid + + # On node rec + r_node = [750, 750] + rec_n = solver_lin.geometry.new_src(name='rec', src_type=None) + rec_n.coordinates.data[:] = r_node + + # Half node rec for linear + r_mid = [755, 755] + rec_hl = solver_lin.geometry.new_src(name='recl', src_type=None) + rec_hl.coordinates.data[:] = r_mid + + # Half node rec for sinc + r_mid = [755, 755] + rec_hs = solver_lin.geometry.new_src(name='recs', src_type=None) + rec_hs.coordinates.data[:] = r_mid + + # Reference solution, on node + _, un, _ = solver_lin.forward(src=src_n, rec=rec_n) + # Linear interp on half node + _, ul, _ = solver_lin.forward(src=src_h, rec=rec_hl) + # Sinc interp on half node + _, us, _ = solver_sinc.forward(src=src_h, rec=rec_hs) + + # Check sinc is more accuracte + nref = np.linalg.norm(rec_n.data) + err_lin = np.linalg.norm(rec_n.data - rec_hl.data)/nref + err_sinc = np.linalg.norm(rec_n.data - rec_hs.data)/nref + + print(f"Error linear: {err_lin}, Error sinc: {err_sinc}") + assert np.isclose(err_sinc, 0, rtol=0, atol=tol) + assert err_sinc < err_lin + assert err_lin > 0.01 From 4aadbf86383c38c182d69ee4a19be83e8bd57168 Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 1 Apr 2024 15:23:47 -0400 Subject: [PATCH 3/8] mpi: revamp conftest --- conftest.py | 56 ++++++++++++++++++++++++++----------- tests/test_autotuner.py | 1 - tests/test_builtins.py | 5 ---- tests/test_dse.py | 1 - tests/test_interpolation.py | 4 +-- tests/test_mpi.py | 4 +-- tests/test_operator.py | 1 - tests/test_pickle.py | 12 ++++---- 8 files changed, 47 insertions(+), 37 deletions(-) diff --git a/conftest.py b/conftest.py index cd74558788..3a4d6d4ef4 100644 --- a/conftest.py +++ b/conftest.py @@ -34,16 +34,12 @@ def skipif(items, whole_module=False): accepted.update({'device', 'device-C', 'device-openmp', 'device-openacc', 'device-aomp', 'cpu64-icc', 'cpu64-icx', 'cpu64-nvc', 'cpu64-arm', 'cpu64-icpx', 'chkpnt'}) - accepted.update({'nompi', 'nodevice'}) + accepted.update({'nodevice'}) unknown = sorted(set(items) - accepted) if unknown: raise ValueError("Illegal skipif argument(s) `%s`" % unknown) skipit = False for i in items: - # Skip if no MPI - if i == 'nompi' and MPI is None: - skipit = "mpi4py/MPI not installed" - break # Skip if won't run on GPUs if i == 'device' and isinstance(configuration['platform'], Device): skipit = "device `%s` unsupported" % configuration['platform'].name @@ -171,6 +167,9 @@ def parallel(item): os.environ['DEVITO_MPI'] = scheme try: check_call(call) + return True + except: + return False finally: os.environ['DEVITO_MPI'] = '0' @@ -189,18 +188,21 @@ def pytest_runtest_setup(item): partest = int(partest) except ValueError: pass - if item.get_closest_marker("parallel") and not partest: - # Blow away function arg in "master" process, to ensure - # this test isn't run on only one process - dummy_test = lambda *args, **kwargs: True - # For pytest <7 - if item.cls is not None: - attr = item.originalname or item.name - setattr(item.cls, attr, dummy_test) + if item.get_closest_marker("parallel"): + if MPI is None: + pytest.skip("mpi4py/MPI not installed") else: - item.obj = dummy_test - # For pytest >= 7 - setattr(item, '_obj', dummy_test) + # Blow away function arg in "master" process, to ensure + # this test isn't run on only one process + dummy_test = lambda *args, **kwargs: True + # For pytest <7 + if item.cls is not None: + attr = item.originalname or item.name + setattr(item.cls, attr, dummy_test) + else: + item.obj = dummy_test + # For pytest >= 7 + setattr(item, '_obj', dummy_test) def pytest_runtest_call(item): @@ -211,7 +213,27 @@ def pytest_runtest_call(item): pass if item.get_closest_marker("parallel") and not partest: # Spawn parallel processes to run test - parallel(item) + passed = parallel(item) + if not passed: + pytest.fail(f"{item} failed in parallel execution") + else: + pytest.skip(f"{item}t passed in parallel execution") + + +@pytest.hookimpl(tryfirst=True, hookwrapper=True) +def pytest_runtest_makereport(item, call): + outcome = yield + result = outcome.get_result() + + partest = os.environ.get('DEVITO_MPI', 0) + try: + partest = int(partest) + except ValueError: + pass + + if item.get_closest_marker("parallel") and not partest: + if call.when == 'call' and result.outcome == 'skipped': + result.outcome = 'passed' # A list of optimization options/pipelines to be used in testing diff --git a/tests/test_autotuner.py b/tests/test_autotuner.py index 496afa3666..72233d3fa0 100644 --- a/tests/test_autotuner.py +++ b/tests/test_autotuner.py @@ -180,7 +180,6 @@ def test_discarding_runs(): assert op._state['autotuning'][1]['tuned']['nthreads'] == 1 -@skipif('nompi') @pytest.mark.parallel(mode=[(2, 'diag'), (2, 'full')]) def test_at_w_mpi(): """Make sure autotuning works in presence of MPI. MPI ranks work diff --git a/tests/test_builtins.py b/tests/test_builtins.py index 4ffe02b552..d086415376 100644 --- a/tests/test_builtins.py +++ b/tests/test_builtins.py @@ -3,7 +3,6 @@ from scipy.ndimage import gaussian_filter from scipy.misc import ascent -from conftest import skipif from devito import ConditionalDimension, Grid, Function, TimeFunction, switchconfig from devito.builtins import (assign, norm, gaussian_smooth, initialize_function, inner, mmin, mmax, sum, sumall) @@ -92,7 +91,6 @@ def test_assign_subsampled_timefunction(self): assert np.all(f.data == 1) - @skipif('nompi') @pytest.mark.parallel(mode=4) def test_assign_parallel(self): a = np.arange(64).reshape((8, 8)) @@ -175,7 +173,6 @@ def test_gs_2d_float(self, sigma): assert np.amax(np.abs(sp_smoothed - np.array(dv_smoothed))) <= 1e-5 - @skipif('nompi') @pytest.mark.parallel(mode=[(4, 'full')]) def test_gs_parallel(self): a = np.arange(64).reshape((8, 8)) @@ -238,7 +235,6 @@ def test_nbl_zero(self): assert np.all(a[:] - np.array(f.data[:]) == 0) - @skipif('nompi') @pytest.mark.parallel(mode=4) def test_if_parallel(self): a = np.arange(36).reshape((6, 6)) @@ -294,7 +290,6 @@ def test_if_halo(self, ndim, nbl): assert np.take(f._data_with_outhalo, 0, axis=-1)[7] == 1 assert np.take(f._data_with_outhalo, -1, axis=-1)[7] == 3 - @skipif('nompi') @pytest.mark.parametrize('nbl', [0, 2]) @pytest.mark.parallel(mode=4) def test_if_halo_mpi(self, nbl): diff --git a/tests/test_dse.py b/tests/test_dse.py index acb3a1721b..9435d58c54 100644 --- a/tests/test_dse.py +++ b/tests/test_dse.py @@ -2799,7 +2799,6 @@ def test_fullopt(self): vexpanded = 2 if configuration['language'] == 'openmp' else 0 assert len(FindNodes(VExpanded).visit(pbs['x0_blk0'])) == vexpanded - @skipif(['nompi']) @switchconfig(profiling='advanced') @pytest.mark.parallel(mode=[(1, 'full')]) def test_fullopt_w_mpi(self): diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 664b882f4f..34bd39bbe0 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -753,7 +753,7 @@ def test_inject_function(): assert u.data[1, i, j] == 0 -@pytest.mark.parametrize('r, interp', [(2, 'linear'), (4, 'cubic')]) +@pytest.mark.parametrize('r, interp', [(2, 'linear'), (4, 'sinc')]) def test_interpolation_radius(r, interp): nt = 11 @@ -778,7 +778,7 @@ def test_interp_default(): src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1, interpolation='sinc') assert isinstance(src.interpolator, SincInterpolator) - assert src.r == 2 + assert src.r == 4 src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1, interpolation='sinc', r=6) diff --git a/tests/test_mpi.py b/tests/test_mpi.py index 12733de0a2..d6cb431a90 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -2,7 +2,7 @@ import pytest from cached_property import cached_property -from conftest import skipif, _R, assert_blocking, assert_structure +from conftest import _R, assert_blocking, assert_structure from devito import (Grid, Constant, Function, TimeFunction, SparseFunction, SparseTimeFunction, Dimension, ConditionalDimension, SubDimension, SubDomain, Eq, Ne, Inc, NODE, Operator, norm, inner, configuration, @@ -19,8 +19,6 @@ from examples.seismic.acoustic import acoustic_setup -pytestmark = skipif(['nompi'], whole_module=True) - class TestDistributor(object): diff --git a/tests/test_operator.py b/tests/test_operator.py index 73fb633a88..5698f8f913 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -1207,7 +1207,6 @@ def test_incomplete_override(self): except: assert False - @skipif('nompi') @pytest.mark.parallel(mode=1) def test_new_distributor(self): """ diff --git a/tests/test_pickle.py b/tests/test_pickle.py index 6e2e851ae7..1089ca5bf2 100644 --- a/tests/test_pickle.py +++ b/tests/test_pickle.py @@ -5,7 +5,6 @@ import numpy as np from sympy import Symbol -from conftest import skipif from devito import (Constant, Eq, Function, TimeFunction, SparseFunction, Grid, Dimension, SubDimension, ConditionalDimension, IncrDimension, TimeDimension, SteppingDimension, Operator, MPI, Min, solve, @@ -79,10 +78,12 @@ def test_function(self, pickle): assert f.dtype == new_f.dtype assert f.shape == new_f.shape - def test_sparse_function(self, pickle): + @pytest.mark.parametrize('interp', ['linear', 'sinc']) + def test_sparse_function(self, pickle, interp): grid = Grid(shape=(3,)) sf = SparseFunction(name='sf', grid=grid, npoint=3, space_order=2, - coordinates=[(0.,), (1.,), (2.,)]) + coordinates=[(0.,), (1.,), (2.,)], + interpolation=interp) sf.data[0] = 1. pkl_sf = pickle.dumps(sf) @@ -91,6 +92,7 @@ def test_sparse_function(self, pickle): # .data is initialized, so it should have been pickled too assert np.all(sf.data[0] == 1.) assert np.all(new_sf.data[0] == 1.) + assert new_sf.interpolation == interp # coordinates should also have been pickled assert np.all(sf.coordinates.data == new_sf.coordinates.data) @@ -633,7 +635,6 @@ def test_elemental(self, pickle): assert str(op) == str(new_op) - @skipif(['nompi']) @pytest.mark.parallel(mode=[1]) def test_mpi_objects(self, pickle): grid = Grid(shape=(4, 4, 4)) @@ -682,7 +683,6 @@ def test_threadid(self, pickle): assert tid.symbolic_min.name == new_tid.symbolic_min.name assert tid.symbolic_max.name == new_tid.symbolic_max.name - @skipif(['nompi']) @pytest.mark.parallel(mode=[2]) def test_mpi_grid(self, pickle): grid = Grid(shape=(4, 4, 4)) @@ -703,7 +703,6 @@ def test_mpi_grid(self, pickle): assert new_grid.distributor.comm.size == 1 MPI.COMM_WORLD.Barrier() - @skipif(['nompi']) @pytest.mark.parallel(mode=[(1, 'full')]) def test_mpi_fullmode_objects(self, pickle): grid = Grid(shape=(4, 4, 4)) @@ -743,7 +742,6 @@ def test_mpi_fullmode_objects(self, pickle): assert v[0] is d.symbolic_min assert v[1] == Min(d.symbolic_max, d.symbolic_min) - @skipif(['nompi']) @pytest.mark.parallel(mode=[(1, 'basic'), (1, 'full')]) def test_mpi_operator(self, pickle): grid = Grid(shape=(4,)) From c4ac44ee942253cab095374c502ee1ea7d90006e Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 1 Apr 2024 22:40:05 -0400 Subject: [PATCH 4/8] mpi: fix sinc setup with mpi and add back missing mpi example tests --- .github/workflows/pytest-core-mpi.yml | 11 ++++++++++- devito/operations/interpolators.py | 8 +++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pytest-core-mpi.yml b/.github/workflows/pytest-core-mpi.yml index 923971f89b..05100d4d2d 100644 --- a/.github/workflows/pytest-core-mpi.yml +++ b/.github/workflows/pytest-core-mpi.yml @@ -48,6 +48,11 @@ jobs: python3 scripts/clear_devito_cache.py python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml -m parallel tests/ + - name: Test examples + run: | + python3 scripts/clear_devito_cache.py + python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml examples/seismic + - name: Upload coverage to Codecov uses: codecov/codecov-action@v4 with: @@ -78,4 +83,8 @@ jobs: - name: Test with pytest run: | - docker run --rm -e CODECOV_TOKEN=${{ secrets.CODECOV_TOKEN }} -e OMP_NUM_THREADS=1 --name testrun devito_img pytest tests/test_mpi.py \ No newline at end of file + docker run --rm -e CODECOV_TOKEN=${{ secrets.CODECOV_TOKEN }} -e OMP_NUM_THREADS=1 --name testrun devito_img pytest tests/test_mpi.py + + - name: Test examples with MPI + run: | + docker run --rm --env DEVITO_MPI=1 --name examplerun mpiexec -n 2 pytest examples/seismic diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index 6871b69da3..dddb2274e7 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -484,7 +484,13 @@ 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 - coords = (coords - np.array(sfunc.grid.origin)) / np.array(sfunc.grid.spacing) + 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.floor(coords) # Precompute sinc From 4db1120225d0ebc703f15ced8f18e75a1660718a Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 2 Apr 2024 10:44:50 -0400 Subject: [PATCH 5/8] api: fix sparse position index for negative location --- .github/workflows/pytest-core-mpi.yml | 6 +++-- devito/operations/interpolators.py | 27 ++++++++++--------- devito/types/sparse.py | 10 ++++--- docker/entrypoint.sh | 2 +- examples/seismic/acoustic/acoustic_example.py | 2 +- examples/seismic/utils.py | 4 ++- 6 files changed, 29 insertions(+), 22 deletions(-) diff --git a/.github/workflows/pytest-core-mpi.yml b/.github/workflows/pytest-core-mpi.yml index 05100d4d2d..b0fb73abc2 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 dddb2274e7..38a3571b17 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 98e0e3a2dd..921c115a22 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 2a3fab2bbd..077a84cd05 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 8a341f8a5c..aced327cdc 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 ef517d7ce7..42ca0ced72 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) From 99448f1d6d970ca2a7aff2699bc93a5a5baa23dc Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 2 Apr 2024 23:15:25 -0400 Subject: [PATCH 6/8] examples: improved fs only check z derivs --- .github/workflows/pytest-core-mpi.yml | 5 ----- examples/seismic/acoustic/operators.py | 27 ++++++++++++++++++++------ 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/.github/workflows/pytest-core-mpi.yml b/.github/workflows/pytest-core-mpi.yml index b0fb73abc2..99d5f5fda2 100644 --- a/.github/workflows/pytest-core-mpi.yml +++ b/.github/workflows/pytest-core-mpi.yml @@ -85,8 +85,3 @@ jobs: - name: Test with pytest run: | docker run --rm -e CODECOV_TOKEN=${{ secrets.CODECOV_TOKEN }} -e OMP_NUM_THREADS=1 --name testrun devito_img pytest tests/test_mpi.py - - - name: Test examples with MPI - run: | - 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/examples/seismic/acoustic/operators.py b/examples/seismic/acoustic/operators.py index 64e7f63284..d06d572c3f 100644 --- a/examples/seismic/acoustic/operators.py +++ b/examples/seismic/acoustic/operators.py @@ -1,5 +1,5 @@ from devito import Eq, Operator, Function, TimeFunction, Inc, solve, sign -from devito.symbolics import retrieve_functions, INT +from devito.symbolics import retrieve_functions, INT, retrieve_derivatives def freesurface(model, eq): @@ -14,13 +14,21 @@ def freesurface(model, eq): eq : Eq Time-stepping stencil (time update) to mirror at the freesurface. """ - lhs, rhs = eq.evaluate.args + lhs, rhs = eq.args # Get vertical dimension and corresponding subdimension - zfs = model.grid.subdomains['fsdomain'].dimensions[-1] + fsdomain = model.grid.subdomains['fsdomain'] + zfs = fsdomain.dimensions[-1] z = zfs.parent - # Functions present in the stencil - funcs = retrieve_functions(rhs) + # Retrieve vertical derivatives + Dz = {d for d in retrieve_derivatives(rhs) if z in d.dims} + # Remove inner duplicate + Dz = Dz - {d for D in Dz for d in retrieve_derivatives(D.expr) if z in d.dims} + Dz = {d: d._eval_at(lhs).evaluate for d in Dz} + + # Finally get functions for evaluated derivatives + funcs = {f for f in retrieve_functions(Dz.values())} + mapper = {} # Antisymmetric mirror at negative indices # TODO: Make a proper "mirror_indices" tool function @@ -29,7 +37,14 @@ def freesurface(model, eq): if (zind - z).as_coeff_Mul()[0] < 0: s = sign(zind.subs({z: zfs, z.spacing: 1})) mapper.update({f: s * f.subs({zind: INT(abs(zind))})}) - return Eq(lhs, rhs.subs(mapper), subdomain=model.grid.subdomains['fsdomain']) + + # Mapper for vertical derivatives + dzmapper = {d: v.subs(mapper) for d, v in Dz.items()} + + fs_eq = [eq.func(lhs, rhs.subs(dzmapper), subdomain=fsdomain)] + fs_eq.append(eq.func(lhs._subs(z, 0), 0, subdomain=fsdomain)) + + return fs_eq def laplacian(field, model, kernel): From a83fac76aeb2e0158ad0e1495d90bfd330cf1fe6 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 3 Apr 2024 08:48:59 -0400 Subject: [PATCH 7/8] misc: review cleanup --- .github/workflows/pytest-core-mpi.yml | 2 +- devito/types/sparse.py | 21 ++++++++++++++------- examples/seismic/acoustic/operators.py | 10 +++++----- 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/.github/workflows/pytest-core-mpi.yml b/.github/workflows/pytest-core-mpi.yml index 99d5f5fda2..a3963fe866 100644 --- a/.github/workflows/pytest-core-mpi.yml +++ b/.github/workflows/pytest-core-mpi.yml @@ -48,7 +48,7 @@ jobs: python3 scripts/clear_devito_cache.py python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml -m parallel tests/ - - name: Test examples + - name: Test examples with MPI run: | python3 scripts/clear_devito_cache.py DEVITO_MPI=1 mpirun -n 2 python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml examples/seismic/acoustic diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 921c115a22..a5f711f7cf 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -812,14 +812,9 @@ class SparseFunction(AbstractSparseFunction): def __init_finalize__(self, *args, **kwargs): # Interpolation method - interp = kwargs.pop('interpolation', 'linear') - self.interpolation = interp - self.interpolator = _interpolators[interp](self) - 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") + self.__interp_setup__(**kwargs) - # Set space ordert to `r` for safety + # Initialization super().__init_finalize__(*args, **kwargs) # Set up sparse point coordinates @@ -827,6 +822,18 @@ def __init_finalize__(self, *args, **kwargs): self._coordinates = self.__subfunc_setup__(coordinates, 'coords') self._dist_origin = {self._coordinates: self.grid.origin_offset} + def __interp_setup__(self, interp='linear', r=None, **kwargs): + self.interpolation = interp + self.interpolator = _interpolators[interp](self) + self._radius = r or _default_radius[interp] + if interp == 'sinc': + if self._radius < 2: + raise ValueError("'sinc' interpolator requires a radius of at least 2") + elif self._radius > 10: + raise ValueError("'sinc' interpolator requires a radius of at most 10") + elif interp == 'linear' and self._radius != 1: + self._radius = 1 + @cached_property def _coordinate_symbols(self): """Symbol representing the coordinate values in each Dimension.""" diff --git a/examples/seismic/acoustic/operators.py b/examples/seismic/acoustic/operators.py index d06d572c3f..c8a2fdec5c 100644 --- a/examples/seismic/acoustic/operators.py +++ b/examples/seismic/acoustic/operators.py @@ -21,13 +21,13 @@ def freesurface(model, eq): z = zfs.parent # Retrieve vertical derivatives - Dz = {d for d in retrieve_derivatives(rhs) if z in d.dims} + dzs = {d for d in retrieve_derivatives(rhs) if z in d.dims} # Remove inner duplicate - Dz = Dz - {d for D in Dz for d in retrieve_derivatives(D.expr) if z in d.dims} - Dz = {d: d._eval_at(lhs).evaluate for d in Dz} + dzs = dzs - {d for D in dzs for d in retrieve_derivatives(D.expr) if z in d.dims} + dzs = {d: d._eval_at(lhs).evaluate for d in dzs} # Finally get functions for evaluated derivatives - funcs = {f for f in retrieve_functions(Dz.values())} + funcs = {f for f in retrieve_functions(dzs.values())} mapper = {} # Antisymmetric mirror at negative indices @@ -39,7 +39,7 @@ def freesurface(model, eq): mapper.update({f: s * f.subs({zind: INT(abs(zind))})}) # Mapper for vertical derivatives - dzmapper = {d: v.subs(mapper) for d, v in Dz.items()} + dzmapper = {d: v.subs(mapper) for d, v in dzs.items()} fs_eq = [eq.func(lhs, rhs.subs(dzmapper), subdomain=fsdomain)] fs_eq.append(eq.func(lhs._subs(z, 0), 0, subdomain=fsdomain)) From 353af4f7d3a38e855f3c88d91af8317d9eda0e10 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 3 Apr 2024 09:15:43 -0400 Subject: [PATCH 8/8] api: switch to recommended scipy i0 when available --- .github/workflows/pytest-core-mpi.yml | 8 ++++++ devito/operations/interpolators.py | 28 +++++++++++++------ devito/types/sparse.py | 12 ++++---- examples/seismic/acoustic/acoustic_example.py | 2 +- examples/seismic/utils.py | 2 ++ 5 files changed, 37 insertions(+), 15 deletions(-) diff --git a/.github/workflows/pytest-core-mpi.yml b/.github/workflows/pytest-core-mpi.yml index a3963fe866..62ae44f511 100644 --- a/.github/workflows/pytest-core-mpi.yml +++ b/.github/workflows/pytest-core-mpi.yml @@ -70,9 +70,12 @@ jobs: - name: gcc arch: gcc os: ubuntu-latest + mpiflag: "" - name: icx arch: icx os: ubuntu-latest + # Need safe math for icx due to inaccuracy with mpi+sinc interpolation + mpiflag: "-e DEVITO_SAFE_MATH=1" steps: - name: Checkout devito @@ -85,3 +88,8 @@ jobs: - name: Test with pytest run: | docker run --rm -e CODECOV_TOKEN=${{ secrets.CODECOV_TOKEN }} -e OMP_NUM_THREADS=1 --name testrun devito_img pytest tests/test_mpi.py + + - name: Test examples with MPI + run: | + docker run --rm ${{ matrix.mpiflag }} -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 \ No newline at end of file diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index 38a3571b17..3abf7669f5 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -5,8 +5,14 @@ import numpy as np from cached_property import cached_property +try: + from scipy.special import i0 +except ImportError: + from numpy import i0 + from devito.finite_differences.differentiable import Mul from devito.finite_differences.elementary import floor +from devito.logger import warning from devito.symbolics import retrieve_function_carriers, retrieve_functions, INT from devito.tools import as_tuple, flatten, filter_ordered from devito.types import (ConditionalDimension, Eq, Inc, Evaluable, Symbol, @@ -216,7 +222,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, pos_only=None): + def _interp_idx(self, variables, implicit_dims=None, pos_only=()): """ Generate interpolation indices for the DiscreteFunctions in ``variables``. """ @@ -238,10 +244,8 @@ def _interp_idx(self, variables, implicit_dims=None, pos_only=None): # 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}) + idx_subs.update({v: v.subs({k: p for (k, p) in zip(mapper, pos)}) + for v in pos_only}) return idx_subs, temps @@ -368,7 +372,7 @@ def _inject(self, field, expr, implicit_dims=None): implicit_dims = implicit_dims + tuple(r.parent for r in self._rdim) # List of indirection indices for all adjacent grid points - idx_subs, temps = self._interp_idx(list(fields), implicit_dims=implicit_dims, + idx_subs, temps = self._interp_idx(fields, implicit_dims=implicit_dims, pos_only=variables) # Substitute coordinate base symbols into the interpolation coefficients @@ -465,6 +469,14 @@ class SincInterpolator(PrecomputedInterpolator): 4: 4.14, 5: 5.26, 6: 6.40, 7: 7.51, 8: 8.56, 9: 9.56, 10: 10.64} + def __init__(self, sfunction): + if i0 is np.i0: + warning(""" +Using `numpy.i0`. We (and numpy) recommend to install scipy to improve the performance +of the SincInterpolator that uses i0 (Bessel function). +""") + super().__init__(sfunction) + @cached_property def interpolation_coeffs(self): coeffs = {} @@ -487,7 +499,7 @@ def _weights(self): def _arg_defaults(self, coords=None, sfunc=None): args = {} b = self._b_table[self.r] - b0 = np.i0(b) + b0 = i0(b) if coords is None or sfunc is None: raise ValueError("No coordinates or sparse function provided") # Coords to indices @@ -499,7 +511,7 @@ def _arg_defaults(self, coords=None, sfunc=None): data = np.zeros((coords.shape[0], 2*self.r), dtype=sfunc.dtype) for ri in range(2*self.r): rpos = ri - self.r + 1 - coords[:, j] - num = np.i0(b*np.sqrt(1 - (rpos/self.r)**2)) + num = i0(b*np.sqrt(1 - (rpos/self.r)**2)) data[:, ri] = num / b0 * np.sinc(rpos) args[self.interpolation_coeffs[r].name] = data diff --git a/devito/types/sparse.py b/devito/types/sparse.py index a5f711f7cf..85910ff955 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -822,16 +822,16 @@ def __init_finalize__(self, *args, **kwargs): self._coordinates = self.__subfunc_setup__(coordinates, 'coords') self._dist_origin = {self._coordinates: self.grid.origin_offset} - def __interp_setup__(self, interp='linear', r=None, **kwargs): - self.interpolation = interp - self.interpolator = _interpolators[interp](self) - self._radius = r or _default_radius[interp] - if interp == 'sinc': + def __interp_setup__(self, interpolation='linear', r=None, **kwargs): + self.interpolation = interpolation + self.interpolator = _interpolators[interpolation](self) + self._radius = r or _default_radius[interpolation] + if interpolation == 'sinc': if self._radius < 2: raise ValueError("'sinc' interpolator requires a radius of at least 2") elif self._radius > 10: raise ValueError("'sinc' interpolator requires a radius of at most 10") - elif interp == 'linear' and self._radius != 1: + elif interpolation == 'linear' and self._radius != 1: self._radius = 1 @cached_property diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py index aced327cdc..bbcc4a7c41 100644 --- a/examples/seismic/acoustic/acoustic_example.py +++ b/examples/seismic/acoustic/acoustic_example.py @@ -103,4 +103,4 @@ def test_isoacoustic(fs, normrec, dtype, interp): run(shape=shape, spacing=spacing, nbl=args.nbl, tn=tn, fs=args.fs, space_order=args.space_order, preset=preset, kernel=args.kernel, autotune=args.autotune, opt=args.opt, full_run=args.full, - checkpointing=args.checkpointing, dtype=args.dtype) + checkpointing=args.checkpointing, dtype=args.dtype, interpolation=args.interp) diff --git a/examples/seismic/utils.py b/examples/seismic/utils.py index 42ca0ced72..835961fa02 100644 --- a/examples/seismic/utils.py +++ b/examples/seismic/utils.py @@ -253,4 +253,6 @@ def __call__(self, parser, args, values, option_string=None): type=float, help="Simulation time in millisecond") parser.add_argument("-dtype", action=_dtype_store, dest="dtype", default=np.float32, choices=['float32', 'float64']) + parser.add_argument("-interp", dest="interp", default="linear", + choices=['linear', 'sinc']) return parser