diff --git a/.github/workflows/pytest-core-mpi.yml b/.github/workflows/pytest-core-mpi.yml index 923971f89b..62ae44f511 100644 --- a/.github/workflows/pytest-core-mpi.yml +++ b/.github/workflows/pytest-core-mpi.yml @@ -48,6 +48,12 @@ jobs: python3 scripts/clear_devito_cache.py python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml -m parallel tests/ + - 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 + 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 with: @@ -64,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 @@ -78,4 +87,9 @@ 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 ${{ 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/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/devito/operations/interpolators.py b/devito/operations/interpolators.py index 98e0105e8a..3abf7669f5 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -2,17 +2,24 @@ from functools import wraps import sympy +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, - CustomDimension) + CustomDimension, SubFunction) from devito.types.utils import DimensionTuple -__all__ = ['LinearInterpolator', 'PrecomputedInterpolator'] +__all__ = ['LinearInterpolator', 'PrecomputedInterpolator', 'SincInterpolator'] def check_radius(func): @@ -22,7 +29,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 @@ -125,6 +132,8 @@ class GenericInterpolator(ABC): Abstract base class defining the interface for an interpolator. """ + _name = "generic" + @abstractmethod def inject(self, *args, **kwargs): pass @@ -133,6 +142,13 @@ def inject(self, *args, **kwargs): def interpolate(self, *args, **kwargs): pass + @property + def name(self): + return self._name + + def _arg_defaults(self, **args): + return {} + class WeightedInterpolator(GenericInterpolator): @@ -142,6 +158,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 @@ -204,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): + def _interp_idx(self, variables, implicit_dims=None, pos_only=()): """ Generate interpolation indices for the DiscreteFunctions in ``variables``. """ @@ -223,6 +241,12 @@ 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] + 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 @@ -347,10 +371,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(fields, implicit_dims=implicit_dims, + pos_only=variables) # Substitute coordinate base symbols into the interpolation coefficients eqns = [Inc(_field.xreplace(idx_subs), @@ -370,6 +393,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 @@ -403,11 +429,13 @@ 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) # 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 @@ -421,3 +449,70 @@ 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 + + """ + + _name = 'sinc' + + # Table 1 + _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} + + 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 = {} + 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 = 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.spacing) + 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 - self.r + 1 - coords[:, j] + num = i0(b*np.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..85910ff955 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': 4} + + class AbstractSparseFunction(DiscreteFunction): """ @@ -799,22 +804,36 @@ 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', 'interpolation') def __init_finalize__(self, *args, **kwargs): + # Interpolation method + self.__interp_setup__(**kwargs) + + # Initialization super().__init_finalize__(*args, **kwargs) - self.interpolator = LinearInterpolator(self) # Set up sparse point coordinates coordinates = kwargs.get('coordinates', kwargs.get('coordinates_data')) self._coordinates = self.__subfunc_setup__(coordinates, 'coords') self._dist_origin = {self._coordinates: self.grid.origin_offset} + 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 interpolation == 'linear' and self._radius != 1: + self._radius = 1 + @cached_property def _coordinate_symbols(self): """Symbol representing the coordinate values in each Dimension.""" @@ -827,6 +846,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.get(key.coordinates.name, key.coordinates.data) + defaults.update(key.interpolator._arg_defaults(coords=coords, + sfunc=key)) + return defaults + class SparseTimeFunction(AbstractSparseTimeFunction, SparseFunction): """ 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 7e922cbb6e..bbcc4a7c41 100644 --- a/examples/seismic/acoustic/acoustic_example.py +++ b/examples/seismic/acoustic/acoustic_example.py @@ -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, @@ -60,21 +60,26 @@ 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)]) @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) @@ -98,4 +103,4 @@ def test_isoacoustic(fs, normrec, dtype): 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/acoustic/operators.py b/examples/seismic/acoustic/operators.py index 3623915a28..c8a2fdec5c 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 +from devito.symbolics import retrieve_functions, INT, retrieve_derivatives def freesurface(model, eq): @@ -15,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 + dzs = {d for d in retrieve_derivatives(rhs) if z in d.dims} + # Remove inner duplicate + 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(dzs.values())} + mapper = {} # Antisymmetric mirror at negative indices # TODO: Make a proper "mirror_indices" tool function @@ -30,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 dzs.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): @@ -119,11 +133,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 +171,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 +215,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 +252,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/utils.py b/examples/seismic/utils.py index e6c0648a83..835961fa02 100644 --- a/examples/seismic/utils.py +++ b/examples/seismic/utils.py @@ -3,23 +3,28 @@ 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 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) + 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 +63,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 +90,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 +147,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 +162,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 +173,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 +189,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} @@ -234,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 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 506a3b88de..34bd39bbe0 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 @@ -752,14 +753,88 @@ def test_inject_function(): assert u.data[1, i, j] == 0 -def test_interpolation_radius(): +@pytest.mark.parametrize('r, interp', [(2, 'linear'), (4, 'sinc')]) +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) + src = SparseTimeFunction(name="src", grid=grid, nt=nt, npoint=1, + 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 == 4 + + 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 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,))