diff --git a/.travis.yml b/.travis.yml index 42cd3a7620..009143ba2e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -29,7 +29,7 @@ matrix: packages: - gcc-5 - g++-5 - env: DEVITO_ARCH=gcc-5 DEVITO_OPENMP=0 RUN_EXAMPLES=True INSTALL_TYPE=conda MPI_INSTALL=1 + env: DEVITO_ARCH=gcc-5 DEVITO_OPENMP=0 RUN_EXAMPLES=False INSTALL_TYPE=conda MPI_INSTALL=0 - os: linux python: "3.6" addons: @@ -47,9 +47,9 @@ matrix: sources: - ubuntu-toolchain-r-test # For gcc 4.9, 5 and 7 packages: - - gcc-4.9 - - g++-4.9 - env: DEVITO_ARCH=gcc-4.9 DEVITO_OPENMP=1 OMP_NUM_THREADS=2 RUN_EXAMPLES=True INSTALL_TYPE=conda MPI_INSTALL=0 + - gcc-7 + - g++-7 + env: DEVITO_ARCH=gcc-7 DEVITO_OPENMP=0 DEVITO_BACKEND=yask YC_CXX=g++-7 INSTALL_TYPE=conda RUN_EXAMPLES=False MPI_INSTALL=1 - os: linux python: "3.6" addons: @@ -59,7 +59,7 @@ matrix: packages: - gcc-7 - g++-7 - env: DEVITO_ARCH=gcc-7 DEVITO_OPENMP=0 DEVITO_BACKEND=yask YC_CXX=g++-7 INSTALL_TYPE=conda RUN_EXAMPLES=False MPI_INSTALL=1 + env: DEVITO_ARCH=gcc-7 DEVITO_OPENMP=0 DEVITO_BACKEND=ops YC_CXX=g++-7 INSTALL_TYPE=conda RUN_EXAMPLES=False MPI_INSTALL=0 allow_failures: - os: linux python: "2.7" @@ -87,14 +87,16 @@ before_install: - conda update -q conda # Useful for debugging any issues with conda - conda info -a - - sudo apt-get install -y -q mpich libmpich-dev + - if [[ "$MPI_INSTALL" == '1' ]]; then + sudo apt-get install -y -q mpich libmpich-dev; + fi install: # Install devito with conda - if [[ $INSTALL_TYPE == 'conda' ]]; then conda env create -q -f environment.yml python=$TRAVIS_PYTHON_VERSION; source activate devito; - pip install -e .; + if [[ $MPI_INSTALL == '1' ]]; then pip install -e .[extras]; else pip install -e .; fi conda list; fi - if [[ "$MPI_INSTALL" == '1' ]]; then diff --git a/Jenkinsfile b/Jenkinsfile index 555efeb324..c565e1c783 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -19,6 +19,7 @@ pipeline { environment { HOME="${WORKSPACE}" DEVITO_OPENMP=0 + INSTALL_MPI=1 PYTHONPATH="${WORKSPACE}/lib/python3.6/site-packages/" } steps { @@ -27,13 +28,14 @@ pipeline { runPipTests() } } - stage('Build and test gcc-4.9 OpenMP container') { + stage('Build and test gcc-4.9 OpenMP container no MPI') { agent { dockerfile { label 'azure-linux-8core' filename 'Dockerfile.jenkins' additionalBuildArgs "--build-arg gccvers=4.9" } } environment { HOME="${WORKSPACE}" DEVITO_OPENMP=1 + INSTALL_MPI=0 OMP_NUM_THREADS=2 } steps { @@ -52,17 +54,19 @@ pipeline { environment { HOME="${WORKSPACE}" DEVITO_OPENMP=0 + INSTALL_MPI=1 } steps { cleanWorkspace() condaInstallDevito() + condaInstallMPI() runCondaTests() runExamples() runCodecov() buildDocs() } } - stage('Build and test gcc-7 container') { + stage('Build and test gcc-7 YASK container') { agent { dockerfile { label 'azure-linux' filename 'Dockerfile.jenkins' additionalBuildArgs "--build-arg gccvers=7" } } @@ -70,17 +74,35 @@ pipeline { HOME="${WORKSPACE}" DEVITO_BACKEND="yask" DEVITO_OPENMP="0" + INSTALL_MPI="1" YC_CXX="g++-7" } steps { cleanWorkspace() condaInstallDevito() + condaInstallMPI() installYask() runCondaTests() runCodecov() buildDocs() } } + stage('Build and test gcc-7 OPS container') { + agent { dockerfile { label 'azure-linux' + filename 'Dockerfile.jenkins' + additionalBuildArgs "--build-arg gccvers=7" } } + environment { + HOME="${WORKSPACE}" + DEVITO_BACKEND="ops" + } + steps { + cleanWorkspace() + condaInstallDevito() + runCondaTests() + runCodecov() + buildDocs() + } + } stage('Build and test gcc-8 container') { agent { dockerfile { label 'azure-linux-8core' filename 'Dockerfile.jenkins' @@ -88,10 +110,12 @@ pipeline { environment { HOME="${WORKSPACE}" DEVITO_OPENMP=0 + INSTALL_MPI=1 } steps { cleanWorkspace() condaInstallDevito() + condaInstallMPI() runCondaTests() runExamples() runCodecov() @@ -114,6 +138,10 @@ def condaInstallDevito () { sh 'source activate devito ; flake8 --exclude .conda,.git --builtins=ArgumentError .' } +def condaInstallMPI () { + sh 'source activate devito ; pip install -r requirements-optional.txt' +} + def pipInstallDevito () { sh "mkdir -p ${WORKSPACE}/lib/python3.6/site-packages/" sh "python setup.py install --prefix=${WORKSPACE}" diff --git a/README.md b/README.md index dedd01381d..933fd0090f 100644 --- a/README.md +++ b/README.md @@ -53,6 +53,9 @@ docker-compose run devito /tests # start a jupyter notebook server on port 8888 docker-compose up devito + +# start a bash shell with devito +docker-compose run devito /bin/bash ``` ## Examples diff --git a/devito/__init__.py b/devito/__init__.py index 1dbdcd7a5d..c4c9380109 100644 --- a/devito/__init__.py +++ b/devito/__init__.py @@ -8,6 +8,7 @@ import cpuinfo from devito.base import * # noqa +from devito.builtins import * # noqa from devito.data.allocators import * # noqa from devito.dimension import * # noqa from devito.equation import * # noqa @@ -27,10 +28,18 @@ __version__ = get_versions()['version'] del get_versions +# Setup compiler and backend configuration.add('compiler', 'custom', list(compiler_registry), callback=lambda i: compiler_registry[i]()) configuration.add('backend', 'core', list(backends_registry), callback=init_backend) +# Should Devito run a first-touch Operator upon allocating data? +configuration.add('first-touch', 0, [0, 1], lambda i: bool(i), False) + +# Should Devito ignore any unknown runtime arguments supplied to Operator.apply(), +# or rather raise an exception (the default behaviour)? +configuration.add('ignore-unknowns', 0, [0, 1], lambda i: bool(i), False) + # Execution mode setup def _reinit_compiler(val): # noqa # Force re-build the compiler @@ -43,7 +52,7 @@ def _reinit_compiler(val): # noqa # Autotuning setup AT_LEVELs = ['off', 'basic', 'aggressive'] AT_MODEs = ['preemptive', 'runtime'] -at_default_mode = {'core': 'preemptive', 'yask': 'runtime'} +at_default_mode = {'core': 'preemptive', 'yask': 'runtime', 'ops': 'runtime'} at_setup = namedtuple('at_setup', 'level mode') at_accepted = AT_LEVELs + [list(i) for i in product(AT_LEVELs, AT_MODEs)] def _at_callback(val): # noqa @@ -58,10 +67,11 @@ def _at_callback(val): # noqa return at_setup(level, 'preemptive') else: return at_setup(level, mode) -configuration.add('autotuning', 'off', at_accepted, callback=_at_callback) # noqa +configuration.add('autotuning', 'off', at_accepted, callback=_at_callback, # noqa + impacts_jit=False) # Should Devito emit the JIT compilation commands? -configuration.add('debug_compiler', 0, [0, 1], lambda i: bool(i)) +configuration.add('debug-compiler', 0, [0, 1], lambda i: bool(i), False) # Set the Instruction Set Architecture (ISA) ISAs = ['cpp', 'avx', 'avx2', 'avx512'] diff --git a/devito/backends.py b/devito/backends.py index 2b7c065cc4..d979cfac0a 100644 --- a/devito/backends.py +++ b/devito/backends.py @@ -14,7 +14,7 @@ backends = {} -backends_registry = ('core', 'yask', 'void') +backends_registry = ('core', 'yask', 'void', 'ops') class void(object): diff --git a/devito/builtins.py b/devito/builtins.py new file mode 100644 index 0000000000..1dd02186a8 --- /dev/null +++ b/devito/builtins.py @@ -0,0 +1,193 @@ +""" +Built-in :class:`Operator`s provided by Devito. +""" + +from sympy import Abs, Pow +import numpy as np + +import devito as dv + +__all__ = ['assign', 'smooth', 'norm', 'sumall', 'inner'] + + +def assign(f, v=0): + """ + Assign a value to a :class:`Function`. + + Parameters + ---------- + f : Function + The left-hand side of the assignment. + v : scalar, optional + The right-hand side of the assignment. + """ + dv.Operator(dv.Eq(f, v), name='assign')() + + +def smooth(f, g, axis=None): + """ + Smooth a :class:`Function` through simple moving average. + + Parameters + ---------- + f : Function + The left-hand side of the smoothing kernel, that is the smoothed Function. + g : Function + The right-hand side of the smoothing kernel, that is the Function being smoothed. + axis : Dimension or list of Dimensions, optional + The :class:`Dimension` along which the smoothing operation is performed. + Defaults to ``f``'s innermost Dimension. + + Notes + ----- + More info about simple moving average available at: :: + + https://en.wikipedia.org/wiki/Moving_average#Simple_moving_average + """ + if g.is_Constant: + # Return a scaled version of the input if it's a Constant + f.data[:] = .9 * g.data + else: + if axis is None: + axis = g.dimensions[-1] + dv.Operator(dv.Eq(f, g.avg(dims=axis)), name='smoother')() + + +# Reduction-inducing builtins + +class MPIReduction(object): + """ + A context manager to build MPI-aware reduction Operators. + """ + + def __init__(self, *functions): + grids = {f.grid for f in functions} + if len(grids) == 0: + self.grid = None + elif len(grids) == 1: + self.grid = grids.pop() + else: + raise ValueError("Multiple Grids found") + dtype = {f.dtype for f in functions} + if len(dtype) == 1: + self.dtype = dtype.pop() + else: + raise ValueError("Illegal mixed data types") + self.v = None + + def __enter__(self): + i = dv.Dimension(name='i',) + self.n = dv.Function(name='n', shape=(1,), dimensions=(i,), + grid=self.grid, dtype=self.dtype) + self.n.data[0] = 0 + return self + + def __exit__(self, exc_type, exc_value, traceback): + if self.grid is None or not dv.configuration['mpi']: + assert self.n.data.size == 1 + self.v = self.n.data[0] + else: + comm = self.grid.distributor.comm + self.v = comm.allreduce(np.asarray(self.n.data))[0] + + +def norm(f, order=2): + """ + Compute the norm of a :class:`Function`. + + Parameters + ---------- + f : Function + Input Function. + order : int, optional + The order of the norm. Defaults to 2. + """ + kwargs = {} + if f.is_TimeFunction and f._time_buffering: + kwargs[f.time_dim.max_name] = f._time_size - 1 + + # Protect SparseFunctions from accessing duplicated (out-of-domain) data, + # otherwise we would eventually be summing more than expected + p, eqns = f.guard() if f.is_SparseFunction else (f, []) + + with MPIReduction(f) as mr: + op = dv.Operator(eqns + [dv.Inc(mr.n[0], Abs(Pow(p, order)))], + name='norm%d' % order) + op.apply(**kwargs) + + v = Pow(mr.v, 1/order) + + return np.float(v) + + +def sumall(f): + """ + Compute the sum of the values in a :class:`Function`. + + Parameters + ---------- + f : Function + Input Function. + """ + kwargs = {} + if f.is_TimeFunction and f._time_buffering: + kwargs[f.time_dim.max_name] = f._time_size - 1 + + # Protect SparseFunctions from accessing duplicated (out-of-domain) data, + # otherwise we would eventually be summing more than expected + p, eqns = f.guard() if f.is_SparseFunction else (f, []) + + with MPIReduction(f) as mr: + op = dv.Operator(eqns + [dv.Inc(mr.n[0], p)], name='sum') + op.apply(**kwargs) + + return np.float(mr.v) + + +def inner(f, g): + """ + Inner product of two :class:`Function`s. + + Parameters + ---------- + f : Function + First input operand + g : Function + Second input operand + + Raises + ------ + ValueError + If the two input Functions are defined over different grids, or have + different dimensionality, or their dimension-wise sizes don't match. + If in input are two SparseFunctions and their coordinates don't match, + the exception is raised. + + Notes + ----- + The inner product is the sum of all dimension-wise products. For 1D Functions, + the inner product corresponds to the dot product. + """ + # Input check + if f.is_TimeFunction and f._time_buffering != g._time_buffering: + raise ValueError("Cannot compute `inner` between save/nosave TimeFunctions") + if f.shape != g.shape: + raise ValueError("`f` and `g` must have same shape") + if f._data is None or g._data is None: + raise ValueError("Uninitialized input") + if f.is_SparseFunction and not np.all(f.coordinates_data == g.coordinates_data): + raise ValueError("Non-matching coordinates") + + kwargs = {} + if f.is_TimeFunction and f._time_buffering: + kwargs[f.time_dim.max_name] = f._time_size - 1 + + # Protect SparseFunctions from accessing duplicated (out-of-domain) data, + # otherwise we would eventually be summing more than expected + rhs, eqns = f.guard(f*g) if f.is_SparseFunction else (f*g, []) + + with MPIReduction(f, g) as mr: + op = dv.Operator(eqns + [dv.Inc(mr.n[0], rhs)], name='inner') + op.apply(**kwargs) + + return np.float(mr.v) diff --git a/devito/compiler.py b/devito/compiler.py index c5e132159e..8e5ffcaf56 100644 --- a/devito/compiler.py +++ b/devito/compiler.py @@ -332,7 +332,7 @@ def save(soname, binary, compiler): debug("%s: `%s` was not saved in `%s` as it already exists" % (compiler, sofile.name, get_jit_dir())) else: - with open(str(path), 'wb') as f: + with open(str(sofile), 'wb') as f: f.write(binary) debug("%s: `%s` successfully saved in `%s`" % (compiler, sofile.name, get_jit_dir())) @@ -353,16 +353,25 @@ def jit_compile(soname, code, compiler): target = str(get_jit_dir().joinpath(soname)) src_file = "%s.%s" % (target, compiler.src_ext) + # This makes a suite of cache directories based on the soname + cache_dir = get_codepy_dir().joinpath(soname[:7]) + cache_dir.mkdir(parents=True, exist_ok=True) + # `catch_warnings` suppresses codepy complaining that it's taking # too long to acquire the cache lock. This warning can only appear # in a multiprocess session, typically (but not necessarily) when # many processes are frequently attempting jit-compilation (e.g., # when running the test suite in parallel) with warnings.catch_warnings(): + warnings.simplefilter('ignore') + tic = time() + # Spinlock in case of MPI + sleep_delay = 0 if configuration['mpi'] else 1 _, _, _, recompiled = compile_from_string(compiler, target, code, src_file, - cache_dir=get_codepy_dir(), - debug=configuration['debug_compiler']) + cache_dir=cache_dir, + debug=configuration['debug-compiler'], + sleep_delay=sleep_delay) toc = time() if recompiled: diff --git a/devito/data/allocators.py b/devito/data/allocators.py index 5fffd63b67..8632fad4b7 100644 --- a/devito/data/allocators.py +++ b/devito/data/allocators.py @@ -7,15 +7,13 @@ import ctypes from ctypes.util import find_library -from devito.equation import Eq from devito.logger import logger from devito.parameters import configuration from devito.tools import numpy_to_ctypes -import devito __all__ = ['ALLOC_FLAT', 'ALLOC_NUMA_LOCAL', 'ALLOC_NUMA_ANY', 'ALLOC_KNL_MCDRAM', 'ALLOC_KNL_DRAM', 'ALLOC_GUARD', - 'default_allocator', 'first_touch'] + 'default_allocator'] class MemoryAllocator(object): @@ -322,11 +320,3 @@ def default_allocator(): return ALLOC_NUMA_LOCAL else: return ALLOC_FLAT - - -def first_touch(array): - """ - Uses an Operator to initialize the given array in the same pattern that - would later be used to access it. - """ - devito.Operator(Eq(array, 0.))() diff --git a/devito/data/data.py b/devito/data/data.py index d7b6b7b414..6611757b9d 100644 --- a/devito/data/data.py +++ b/devito/data/data.py @@ -12,56 +12,45 @@ class Data(np.ndarray): """ - A special :class:`numpy.ndarray` allowing logical indexing. - - The type :class:`numpy.ndarray` is subclassed as indicated at: :: + A :class:`numpy.ndarray` supporting distributed dimensions. + + Parameters + ---------- + shape : tuple of ints + Shape of created array. + dtype : numpy.dtype + The data type of the raw data. + decomposition : tuple of :class:`Decomposition`, optional + The data decomposition, for each dimension. + modulo : tuple of bool, optional + If the i-th entry is True, then the i-th array dimension uses modulo indexing. + allocator : :class:`MemoryAllocator`, optional + Used to allocate memory. Defaults to ``ALLOC_FLAT``. + + Notes + ----- + NumPy array subclassing is described at: :: https://docs.scipy.org/doc/numpy-1.13.0/user/basics.subclassing.html - :param shape: Shape of the array in grid points. - :param dimensions: The array :class:`Dimension`s. - :param dtype: A ``numpy.dtype`` for the raw data. - :param decomposition: (Optional) a mapper from :class:`Dimension`s in - ``dimensions`` to :class:`Decomposition`s, which - describe how the Data is distributed over a set - of processes. The Decompositions will be used to - translate global array indices into local indices. - The local indices are relative to the calling process. - This is only relevant in the case of distributed - memory execution (via MPI). - :param allocator: (Optional) a :class:`MemoryAllocator` to specialize memory - allocation. Defaults to ``ALLOC_FLAT``. - - .. note:: - - This type supports logical indexing over modulo buffered dimensions. - - .. note:: - - Any view or copy ``A`` created starting from ``self``, for instance via - a slice operation or a universal function ("ufunc" in NumPy jargon), will - still be of type :class:`Data`. However, if ``A``'s rank is different than - ``self``'s rank, namely if ``A.ndim != self.ndim``, then the capability of - performing logical indexing is lost. + Any view or copy created from ``self``, for instance via a slice operation + or a universal function ("ufunc" in NumPy jargon), will still be of type + :class:`Data`. """ - def __new__(cls, shape, dimensions, dtype, decomposition=None, allocator=ALLOC_FLAT): - assert len(shape) == len(dimensions) + def __new__(cls, shape, dtype, decomposition=None, modulo=None, allocator=ALLOC_FLAT): + assert len(shape) == len(modulo) ndarray, memfree_args = allocator.alloc(shape, dtype) obj = np.asarray(ndarray).view(cls) obj._allocator = allocator obj._memfree_args = memfree_args - obj._decomposition = tuple((decomposition or {}).get(i) for i in dimensions) - obj._modulo = tuple(True if i.is_Stepping else False for i in dimensions) - - # By default, the indices used to access array values are interpreted as - # local indices - obj._glb_indexing = False + obj._decomposition = decomposition or (None,)*len(shape) + obj._modulo = modulo or (False,)*len(shape) # This cannot be a property, as Data objects constructed from this # object might not have any `decomposition`, but they would still be # distributed. Hence, in `__array_finalize__` we must copy this value - obj._is_decomposed = any(i is not None for i in obj._decomposition) + obj._is_distributed = any(i is not None for i in obj._decomposition) # Saves the last index used in `__getitem__`. This allows `__array_finalize__` # to reconstruct information about the computed view (e.g., `decomposition`) @@ -96,18 +85,17 @@ def __array_finalize__(self, obj): if type(obj) != Data: # Definitely from view casting - self._glb_indexing = False - self._is_decomposed = False + self._is_distributed = False self._modulo = tuple(False for i in range(self.ndim)) self._decomposition = (None,)*self.ndim elif obj._index_stash is not None: # From `__getitem__` - self._glb_indexing = obj._glb_indexing - self._is_decomposed = obj._is_decomposed - idx = obj._normalize_index(obj._index_stash) - self._modulo = tuple(m for i, m in zip(idx, obj._modulo) if not is_integer(i)) + self._is_distributed = obj._is_distributed + glb_idx = obj._normalize_index(obj._index_stash) + self._modulo = tuple(m for i, m in zip(glb_idx, obj._modulo) + if not is_integer(i)) decomposition = [] - for i, dec in zip(idx, obj._decomposition): + for i, dec in zip(glb_idx, obj._decomposition): if is_integer(i): continue elif dec is None: @@ -116,8 +104,7 @@ def __array_finalize__(self, obj): decomposition.append(dec.reshape(i)) self._decomposition = tuple(decomposition) else: - self._glb_indexing = obj._glb_indexing - self._is_decomposed = obj._is_decomposed + self._is_distributed = obj._is_distributed if self.ndim == obj.ndim: # E.g., from a ufunc, such as `np.add` self._modulo = obj._modulo @@ -129,21 +116,25 @@ def __array_finalize__(self, obj): @property def _local(self): - """Return a view of ``self`` with disabled global indexing.""" + """A view of ``self`` with global indexing disabled.""" ret = self.view() - ret._glb_indexing = False + ret._is_distributed = False return ret - @property - def _global(self): - """Return a view of ``self`` with enabled global indexing.""" - ret = self.view() - ret._glb_indexing = True + def _global(self, glb_idx, decomposition): + """A "global" view of ``self`` over a given :class:`Decomposition`.""" + if self._is_distributed: + raise ValueError("Cannot derive a decomposed view from a decomposed Data") + if len(decomposition) != self.ndim: + raise ValueError("`decomposition` should have ndim=%d entries" % self.ndim) + ret = self[glb_idx] + ret._decomposition = decomposition + ret._is_distributed = any(i is not None for i in decomposition) return ret @property def _is_mpi_distributed(self): - return self._is_decomposed and configuration['mpi'] + return self._is_distributed and configuration['mpi'] def __repr__(self): return super(Data, self._local).__repr__() @@ -155,7 +146,7 @@ def __getitem__(self, glb_idx): # self's data partition, so None is returned return None else: - self._index_stash = glb_idx # Will be popped in `__array_finalize__` + self._index_stash = glb_idx retval = super(Data, self).__getitem__(loc_idx) self._index_stash = None return retval @@ -166,45 +157,50 @@ def __setitem__(self, glb_idx, val): # no-op return elif np.isscalar(val): - pass - elif isinstance(val, Data) and val._is_decomposed: - if self._is_decomposed: + if index_is_basic(loc_idx): + # Won't go through `__getitem__` as it's basic indexing mode, + # so we should just propage `loc_idx` + super(Data, self).__setitem__(loc_idx, val) + else: + super(Data, self).__setitem__(glb_idx, val) + elif isinstance(val, Data) and val._is_distributed: + if self._is_distributed: # `val` is decomposed, `self` is decomposed -> local set - pass + super(Data, self).__setitem__(glb_idx, val) else: # `val` is decomposed, `self` is replicated -> gatherall-like raise NotImplementedError - elif isinstance(val, Iterable): - if self._is_decomposed: + elif isinstance(val, np.ndarray): + if self._is_distributed: # `val` is replicated, `self` is decomposed -> `val` gets decomposed - if self._glb_indexing: - val_idx = self._normalize_index(glb_idx) - val_idx = [index_dist_to_repl(i, dec) for i, dec in - zip(val_idx, self._decomposition)] - if NONLOCAL in val_idx: - # no-op - return - val_idx = [i for i in val_idx if i is not PROJECTED] - # NumPy broadcasting note: - # When operating on two arrays, NumPy compares their shapes - # element-wise. It starts with the trailing dimensions, and works - # its way forward. Two dimensions are compatible when - # * they are equal, or - # * one of them is 1 - # Conceptually, below we apply the same rule - val_idx = val_idx[len(val_idx)-val.ndim:] - val = val[val_idx] + val_idx = self._normalize_index(glb_idx) + val_idx = [index_dist_to_repl(i, dec) for i, dec in + zip(val_idx, self._decomposition)] + if NONLOCAL in val_idx: + # no-op + return + val_idx = [i for i in val_idx if i is not PROJECTED] + # NumPy broadcasting note: + # When operating on two arrays, NumPy compares their shapes + # element-wise. It starts with the trailing dimensions, and works + # its way forward. Two dimensions are compatible when + # * they are equal, or + # * one of them is 1 + # Conceptually, below we apply the same rule + val_idx = val_idx[len(val_idx)-val.ndim:] + val = val[val_idx] else: # `val` is replicated`, `self` is replicated -> plain ndarray.__setitem__ pass + super(Data, self).__setitem__(glb_idx, val) + elif isinstance(val, Iterable): + if self._is_mpi_distributed: + raise NotImplementedError("With MPI data can only be set " + "via scalars or numpy arrays") + super(Data, self).__setitem__(glb_idx, val) else: raise ValueError("Cannot insert obj of type `%s` into a Data" % type(val)) - # Finally, perform the `__setitem__` - # Note: we pass `glb_idx`, rather than `loc_idx`, as `__setitem__` calls - # `__getitem__`, which in turn expects a global index - super(Data, self).__setitem__(glb_idx, val) - def _normalize_index(self, idx): if isinstance(idx, np.ndarray): # Advanced indexing mode @@ -235,7 +231,7 @@ def _convert_index(self, glb_idx): if mod is True: # Need to wrap index based on modulo v = index_apply_modulo(i, s) - elif self._glb_indexing is True and dec is not None: + elif self._is_distributed is True and dec is not None: # Need to convert the user-provided global indices into local indices. # Obviously this will have no effect if MPI is not used try: @@ -278,7 +274,12 @@ class Index(Tag): def index_is_basic(idx): - return all(is_integer(i) or (i is NONLOCAL) for i in idx) + if is_integer(idx): + return True + elif isinstance(idx, (slice, np.ndarray)): + return False + else: + return all(is_integer(i) or (i is NONLOCAL) for i in idx) def index_apply_modulo(idx, modulo): @@ -310,11 +311,8 @@ def index_dist_to_repl(idx, decomposition): """ Convert a distributed array index a replicated array index. """ - if is_integer(idx): - return PROJECTED - if decomposition is None: - return idx + return PROJECTED if is_integer(idx) else idx # Derive shift value value = idx.start if isinstance(idx, slice) else idx @@ -326,7 +324,9 @@ def index_dist_to_repl(idx, decomposition): # Convert into absolute local index idx = decomposition.convert_index(idx, rel=False) - if idx is None: + if is_integer(idx): + return PROJECTED + elif idx is None: return NONLOCAL elif isinstance(idx, (tuple, list)): return [i - value for i in idx] diff --git a/devito/data/decomposition.py b/devito/data/decomposition.py index 7323c16947..afc789f715 100644 --- a/devito/data/decomposition.py +++ b/devito/data/decomposition.py @@ -41,7 +41,7 @@ class Decomposition(tuple): >>> d = Decomposition([[0, 1, 2], [3, 4], [5, 6, 7]], 1) >>> d - Decomposition([0 1 2], <<[3 4]>>, [5 6 7]) + Decomposition([0,2], <<[3,4]>>, [5,7]) >>> d.loc_abs_min 3 >>> d.loc_abs_max @@ -112,9 +112,14 @@ def __eq__(self, o): all(np.all(i == j) for i, j in zip(self, o)) def __repr__(self): - items = ', '.join('<<%s>>' % str(v) if self.local == i else str(v) - for i, v in enumerate(self)) - return "Decomposition(%s)" % items + ret = [] + for i, v in enumerate(self): + bounds = (min(v, default=None), max(v, default=None)) + item = '[]' if bounds == (None, None) else '[%d,%d]' % bounds + if self.local == i: + item = "<<%s>>" % item + ret.append(item) + return 'Decomposition(%s)' % ', '.join(ret) def __call__(self, *args, rel=True): """Alias for ``self.convert_index``.""" @@ -157,7 +162,7 @@ def convert_index(self, *args, rel=True): >>> d = Decomposition([[0, 1, 2], [3, 4], [5, 6, 7], [8, 9, 10, 11]], 2) >>> d - Decomposition([0 1 2], [3 4], <<[5 6 7]>>, [ 8 9 10 11]) + Decomposition([0,2], [3,4], <<[5,7]>>, [8,11]) A global index as single argument: @@ -292,42 +297,43 @@ def reshape(self, *args): -------- >>> d = Decomposition([[0, 1, 2], [3, 4], [5, 6, 7], [8, 9, 10, 11]], 2) >>> d - Decomposition([0 1 2], [3 4], <<[5 6 7]>>, [ 8 9 10 11]) + Decomposition([0,2], [3,4], <<[5,7]>>, [8,11]) Providing explicit values >>> d.reshape(1, 1) - Decomposition([0 1 2 3], [4 5], <<[6 7 8]>>, [ 9 10 11 12 13]) + Decomposition([0,3], [4,5], <<[6,8]>>, [9,13]) >>> d.reshape(-2, 2) - Decomposition([0], [1 2], <<[3 4 5]>>, [ 6 7 8 9 10 11]) + Decomposition([0,0], [1,2], <<[3,5]>>, [6,11]) Providing a slice >>> d.reshape(slice(2, 9)) - Decomposition([0], [1 2], <<[3 4 5]>>, [6]) + Decomposition([0,0], [1,2], <<[3,5]>>, [6,6]) >>> d.reshape(slice(2, -2)) - Decomposition([0], [1 2], <<[3 4 5]>>, [6 7]) + Decomposition([0,0], [1,2], <<[3,5]>>, [6,7]) >>> d.reshape(slice(4)) - Decomposition([0 1 2], [3], <<[]>>, []) + Decomposition([0,2], [3,3], <<[]>>, []) """ if len(args) == 1: - if isinstance(args[0], slice): - if args[0].start is None: + arg = args[0] + if isinstance(arg, slice): + if arg.start is None or self.glb_min is None: nleft = 0 else: - nleft = self.glb_min - args[0].start - if args[0].stop is None: + nleft = self.glb_min - arg.start + if arg.stop is None or self.glb_max is None: nright = 0 - elif args[0].stop < 0: - nright = args[0].stop + elif arg.stop < 0: + nright = arg.stop else: - nright = args[0].stop - self.glb_max - 1 - elif isinstance(args[0], Iterable): - items = [np.array([j for j in i if j in args[0]]) for i in self] + nright = arg.stop - self.glb_max - 1 + elif isinstance(arg, Iterable): + items = [np.array([j for j in i if j in arg]) for i in self] for i, arr in enumerate(list(items)): items[i] = np.arange(arr.size) + sum(j.size for j in items[:i]) return Decomposition(items, self.local) diff --git a/devito/dle/transformer.py b/devito/dle/transformer.py index 4d6b2e3876..b2cec9c99b 100644 --- a/devito/dle/transformer.py +++ b/devito/dle/transformer.py @@ -25,7 +25,7 @@ This dictionary may be modified at backend-initialization time.""" configuration.add('dle', 'advanced', list(default_modes)) -configuration.add('dle_options', +configuration.add('dle-options', ';'.join('%s:%s' % (k, v) for k, v in default_options.items()), list(default_options)) @@ -75,7 +75,7 @@ def transform(node, mode='basic', options=None): if i not in default_options: dle_warning("Illegal DLE parameter '%s'" % i) params.pop(i) - params.update({k: v for k, v in configuration['dle_options'].items() + params.update({k: v for k, v in configuration['dle-options'].items() if k not in params}) params.update({k: v for k, v in default_options.items() if k not in params}) params['compiler'] = configuration['compiler'] diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index 1711400ca4..4009829a3f 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -380,12 +380,13 @@ def generate_fd_shortcuts(function): derivatives[name_fd] = (deriv, desciption) else: # Left - deriv = partial(first_derivative, order=space_fd_order, dim=d, side=left) + dim_order = time_fd_order if d.is_Time else space_fd_order + deriv = partial(first_derivative, order=dim_order, dim=d, side=left) name_fd = 'd%sl' % name desciption = 'left first order derivative w.r.t dimension %s' % d derivatives[name_fd] = (deriv, desciption) # Right - deriv = partial(first_derivative, order=space_fd_order, dim=d, side=right) + deriv = partial(first_derivative, order=dim_order, dim=d, side=right) name_fd = 'd%sr' % name desciption = 'right first order derivative w.r.t dimension %s' % d derivatives[name_fd] = (deriv, desciption) diff --git a/devito/function.py b/devito/function.py index e38c190680..a01a4eee67 100644 --- a/devito/function.py +++ b/devito/function.py @@ -6,8 +6,9 @@ from psutil import virtual_memory from cached_property import cached_property +from devito.builtins import assign from devito.cgen_utils import INT, cast_mapper -from devito.data import Data, default_allocator, first_touch +from devito.data import Data, default_allocator from devito.dimension import Dimension, ConditionalDimension, DefaultDimension from devito.equation import Eq, Inc from devito.exceptions import InvalidArgument @@ -45,7 +46,7 @@ class Constant(AbstractCachedSymbol, ArgProvider): def __init__(self, *args, **kwargs): if not self._cached(): - self._value = kwargs.get('value') + self._value = kwargs.get('value', 0) @classmethod def __dtype_setup__(cls, **kwargs): @@ -143,7 +144,7 @@ def __init__(self, *args, **kwargs): # Data-related properties and data initialization self._data = None - self._first_touch = kwargs.get('first_touch', configuration['first_touch']) + self._first_touch = kwargs.get('first_touch', configuration['first-touch']) self._allocator = kwargs.get('allocator', default_allocator()) initializer = kwargs.get('initializer') if initializer is None or callable(initializer): @@ -154,7 +155,7 @@ def __init__(self, *args, **kwargs): # a reference to the user-provided buffer self._initializer = None if len(initializer) > 0: - self.data_allocated[:] = initializer + self.data_with_halo[:] = initializer else: # This is a corner case -- we might get here, for example, when # running with MPI and some processes get 0-size arrays after @@ -170,21 +171,21 @@ def _allocate_memory(func): def wrapper(self): if self._data is None: debug("Allocating memory for %s%s" % (self.name, self.shape_allocated)) - self._data = Data(self.shape_allocated, self.indices, self.dtype, - self._decomposition, self._allocator) + self._data = Data(self.shape_allocated, self.dtype, + modulo=self._mask_modulo, allocator=self._allocator) if self._first_touch: - first_touch(self) + assign(self, 0) if callable(self._initializer): if self._first_touch: warning("`first touch` together with `initializer` causing " "redundant data initialization") try: - self._initializer(self._data) + self._initializer(self.data_with_halo) except ValueError: # Perhaps user only wants to initialise the physical domain - self._initializer(self._data[self._mask_domain]) + self._initializer(self.data) else: - self._data.fill(0) + self.data_with_halo.fill(0) return func(self) return wrapper @@ -236,9 +237,9 @@ def __distributor_setup__(self, **kwargs): @property def _data_buffer(self): - """Reference to the data. Unlike ``data, data_with_halo, data_allocated``, + """Reference to the data. Unlike :attr:`data` and :attr:`data_with_halo`, this *never* returns a view of the data. This method is for internal use only.""" - return self.data_allocated + return self._data_allocated @property def _mem_external(self): @@ -255,41 +256,94 @@ def staggered(self): @cached_property def shape(self): """ - Shape of the domain associated with this :class:`TensorFunction`. - The domain constitutes the area of the data written to in a - stencil update. + Shape of the domain region. The domain constitutes the area of the + data written to by an :class:`Operator`. + + Notes + ----- + In an MPI context, this is the *local* domain region shape. """ return self.shape_domain @cached_property def shape_domain(self): """ - Shape of the domain associated with this :class:`TensorFunction`. - The domain constitutes the area of the data written to in a - stencil update. + Shape of the domain region. The domain constitutes the area of the + data written to by an :class:`Operator`. - .. note:: + Notes + ----- + In an MPI context, this is the *local* domain region shape. - Alias to ``self.shape``. + Alias to ``self.shape``. """ return tuple(i - j for i, j in zip(self._shape, self.staggered)) @cached_property def shape_with_halo(self): """ - Shape of the domain plus the read-only stencil boundary associated - with this :class:`Function`. + Shape of the domain+outhalo region. The outhalo is the region + surrounding the domain that may be read by an :class:`Operator`. + + Notes + ----- + In an MPI context, this is the *local* with_halo region shape. + + Further, note that the outhalo of inner ranks is typically empty, while + the outhalo of boundary ranks contains a number of elements depending + on the rank position in the decomposed grid (corner, side, ...). + """ + return tuple(j + i + k for i, (j, k) in zip(self.shape_domain, + self._extent_outhalo)) + + _shape_with_outhalo = shape_with_halo + + @cached_property + def _shape_with_inhalo(self): + """ + Shape of the domain+inhalo region. The inhalo region comprises the + outhalo as well as any additional "ghost" layers for MPI halo + exchanges. Data in the inhalo region are exchanged when running + :class:`Operator`s to maintain consistent values as in sequential runs. + + Notes + ----- + Typically, this property won't be used in user code, but it may come + in handy for testing or debugging """ return tuple(j + i + k for i, (j, k) in zip(self.shape_domain, self._halo)) @cached_property def shape_allocated(self): """ - Shape of the allocated data associated with this :class:`Function`. - It includes the domain and halo regions, as well as any additional - padding outside of the halo. + Shape of the allocated data. It includes the domain and inhalo regions, + as well as any additional padding surrounding the halo. + + Notes + ----- + In an MPI context, this is the *local* with_halo region shape. + """ + return tuple(j + i + k for i, (j, k) in zip(self._shape_with_inhalo, + self._padding)) + + @cached_property + def shape_global(self): + """ + Global shape of the domain region. The domain constitutes the area of + the data written to by an :class:`Operator`. + + Notes + ----- + In an MPI context, this is the *global* domain region shape, which is + therefore identical on all MPI ranks. """ - return tuple(j + i + k for i, (j, k) in zip(self.shape_with_halo, self._padding)) + if self.grid is None: + return self.shape + retval = [] + for d, s in zip(self.dimensions, self.shape): + size = self.grid.dimension_map.get(d) + retval.append(size.glb if size is not None else s) + return tuple(retval) _offset_inhalo = AbstractCachedFunction._offset_halo _extent_inhalo = AbstractCachedFunction._extent_halo @@ -312,6 +366,13 @@ def _extent_outhalo(self): return EnrichedTuple(*extents, getters=self.dimensions, left=left, right=right) + @cached_property + def _mask_modulo(self): + """ + A boolean mask telling which :class:`Dimension`s support modulo-indexing. + """ + return tuple(True if i.is_Stepping else False for i in self.dimensions) + @cached_property def _mask_domain(self): """ @@ -336,6 +397,30 @@ def _mask_outhalo(self): return tuple(slice(i.start - j.left, i.stop and i.stop + j.right or None) for i, j in zip(self._mask_domain, self._extent_outhalo)) + @cached_property + def _decomposition(self): + """ + A tuple of :class:`Decomposition`s, representing the domain + decomposition. None is used as a placeholder for non-decomposed + Dimensions. + """ + if self._distributor is None: + return (None,)*self.ndim + mapper = {d: self._distributor.decomposition[d] for d in self._dist_dimensions} + return tuple(mapper.get(d) for d in self.dimensions) + + @cached_property + def _decomposition_outhalo(self): + """ + A tuple of :class:`Decomposition`s, representing the domain+outhalo + decomposition. None is used as a placeholder for non-decomposed + Dimensions. + """ + if self._distributor is None: + return (None,)*self.ndim + return tuple(v.reshape(*self._extent_inhalo[d]) if v is not None else v + for d, v in zip(self.dimensions, self._decomposition)) + @property def data(self): """ @@ -343,11 +428,11 @@ def data(self): Elements are stored in row-major format. - .. note:: - - With this accessor you are claiming that you will modify - the values you get back. If you only need to look at the - values, use :meth:`data_ro` instead. + Notes + ----- + With this accessor you are claiming that you will modify the values you + get back. If you only need to look at the values, use :meth:`data_ro` + instead. """ return self.data_domain @@ -359,18 +444,16 @@ def data_domain(self): Elements are stored in row-major format. - .. note:: - - With this accessor you are claiming that you will modify - the values you get back. If you only need to look at the - values, use :meth:`data_ro_domain` instead. - - .. note:: + Notes + ----- + Alias to ``self.data``. - Alias to ``self.data``. + With this accessor you are claiming that you will modify the values you + get back. If you only need to look at the values, use + :meth:`data_ro_domain` instead. """ self._is_halo_dirty = True - return self._data[self._mask_domain]._global + return self._data._global(self._mask_domain, self._decomposition) @property @_allocate_memory @@ -380,15 +463,16 @@ def data_with_halo(self): Elements are stored in row-major format. - .. note:: + Notes + ----- - With this accessor you are claiming that you will modify - the values you get back. If you only need to look at the - values, use :meth:`data_ro_with_halo` instead. + With this accessor you are claiming that you will modify the values you + get back. If you only need to look at the values, use + :meth:`data_ro_with_halo` instead. """ self._is_halo_dirty = True self._halo_exchange() - return self._data[self._mask_outhalo]._global + return self._data._global(self._mask_outhalo, self._decomposition_outhalo) _data_with_outhalo = data_with_halo @@ -400,39 +484,43 @@ def _data_with_inhalo(self): Elements are stored in row-major format. - .. note:: + Notes + ----- + This accessor does *not* support global indexing. - With this accessor you are claiming that you will modify - the values you get back. If you only need to look at the - values, use :meth:`data_ro_with_inhalo` instead. + With this accessor you are claiming that you will modify the values you + get back. If you only need to look at the values, use + :meth:`data_ro_with_inhalo` instead. - .. note:: - - Typically, this property won't be used in user code to set - or read data values. Instead, it may come in handy for - testing or debugging + Typically, this accessor won't be used in user code to set or read data + values. Instead, it may come in handy for testing or debugging """ self._is_halo_dirty = True self._halo_exchange() - return self._data[self._mask_inhalo]._global + return np.asarray(self._data[self._mask_inhalo]) @property @_allocate_memory - def data_allocated(self): + def _data_allocated(self): """ - The allocated data values, that is domain+halo+padding. + The allocated data values, that is domain+inhalo+padding. Elements are stored in row-major format. - .. note:: + Notes + ----- + This accessor does *not* support global indexing. + + With this accessor you are claiming that you will modify the values you + get back. If you only need to look at the values, use + :meth:`data_ro_allocated` instead. - With this accessor you are claiming that you will modify - the values you get back. If you only need to look at the - values, use :meth:`data_ro_allocated` instead. + Typically, this accessor won't be used in user code to set or read data + values. Instead, it may come in handy for testing or debugging """ self._is_halo_dirty = True self._halo_exchange() - return self._data._global + return np.asarray(self._data) @property @_allocate_memory @@ -440,15 +528,17 @@ def data_ro_domain(self): """ A read-only view of the domain data values. """ - view = self._data[self._mask_domain]._global + view = self._data._global(self._mask_domain, self._decomposition) view.setflags(write=False) return view @property @_allocate_memory def data_ro_with_halo(self): - """A read-only view of the domain+outhalo data values.""" - view = self._data[self._mask_outhalo]._global + """ + A read-only view of the domain+outhalo data values. + """ + view = self._data._global(self._mask_outhalo, self._decomposition_outhalo) view.setflags(write=False) return view @@ -457,18 +547,30 @@ def data_ro_with_halo(self): @property @_allocate_memory def _data_ro_with_inhalo(self): - """A read-only view of the domain+inhalo data values.""" - view = self._data[self._mask_inhalo]._global + """ + A read-only view of the domain+inhalo data values. + + Notes + ----- + This accessor does *not* support global indexing. + """ + view = self._data[self._mask_inhalo] view.setflags(write=False) - return view + return np.asarray(view) @property @_allocate_memory - def data_ro_allocated(self): - """A read-only view of the domain+halo+padding data values.""" - view = self._data._global + def _data_ro_allocated(self): + """ + A read-only view of the domain+inhalo+padding data values. + + Notes + ----- + This accessor does *not* support global indexing. + """ + view = self._data view.setflags(write=False) - return view + return np.asarray(view) @cached_property def local_indices(self): @@ -476,12 +578,12 @@ def local_indices(self): A tuple of slices representing the global indices that logically belong to the calling MPI rank. - .. note:: - - Given a Function ``f(x, y)`` with shape ``(nx, ny)``, when *not* - using MPI this property will return ``(slice(0, nx-1), slice(0, ny-1))``. - On the other hand, when MPI is used, the local ranges depend on the - domain decomposition, which is carried by ``self.grid``. + Notes + ----- + Given a Function ``f(x, y)`` with shape ``(nx, ny)``, when *not* using + MPI this property will return ``(slice(0, nx-1), slice(0, ny-1))``. On + the other hand, when MPI is used, the local ranges depend on the domain + decomposition, which is carried by ``self.grid``. """ if self._distributor is None: return tuple(slice(0, s) for s in self.shape) @@ -489,15 +591,22 @@ def local_indices(self): return tuple(self._distributor.glb_slices.get(d, slice(0, s)) for s, d in zip(self.shape, self.dimensions)) - @property + @cached_property def space_dimensions(self): """Tuple of :class:`Dimension`s that define physical space.""" return tuple(d for d in self.indices if d.is_Space) + @cached_property + def _dist_dimensions(self): + """Tuple of MPI-distributed :class:`Dimension`s.""" + if self._distributor is None: + return () + return tuple(d for d in self.indices if d in self._distributor.dimensions) + @property def initializer(self): if self._data is not None: - return self._data.view(np.ndarray) + return self.data_with_halo.view(np.ndarray) else: return self._initializer @@ -515,31 +624,13 @@ def symbolic_shape(self): for i, j in zip(symbolic_shape, self.staggered)) return EnrichedTuple(*ret, getters=self.dimensions) - @cached_property - def _decomposition(self): - """ - A mapper from self's distributed :class:`Dimension`s to their - :class:`Decomposition`s. - - .. note:: - - The partitioning encoded in the returned :class:`Decomposition`s - includes the indices falling in the halo+padding regions. - """ - if self._distributor is None: - return {} - mapper = {d: self._distributor.decomposition[d] for d in self.dimensions - if d in self._distributor.dimensions} - # Extend the `Decomposition`s to include the non-domain indices - return {d: v.reshape(*self._offset_domain[d]) for d, v in mapper.items()} - def _halo_exchange(self): """Perform the halo exchange with the neighboring processes.""" if not MPI.Is_initialized() or MPI.COMM_WORLD.size == 1: # Nothing to do return if MPI.COMM_WORLD.size > 1 and self._distributor is None: - raise RuntimeError("`%s` cannot perfom a halo exchange as it has " + raise RuntimeError("`%s` cannot perform a halo exchange as it has " "no Grid attached" % self.name) if self._in_flight: raise RuntimeError("`%s` cannot initiate a halo exchange as previous " @@ -749,30 +840,32 @@ def __indices_setup__(cls, **kwargs): def __shape_setup__(cls, **kwargs): grid = kwargs.get('grid') dimensions = kwargs.get('dimensions') - shape = kwargs.get('shape') + shape = kwargs.get('shape', kwargs.get('shape_global')) if grid is None: if shape is None: raise TypeError("Need either `grid` or `shape`") elif shape is None: - shape = grid.shape_domain + if dimensions is not None and dimensions != grid.dimensions: + raise TypeError("Need `shape` as not all `dimensions` are in `grid`") + shape = grid.shape_local elif dimensions is None: raise TypeError("`dimensions` required if both `grid` and " "`shape` are provided") else: # Got `grid`, `dimensions`, and `shape`. We sanity-check that the - # Dimensions in `dimensions` which also appear in `grid` have - # size (given by `shape`) matching the one in `grid` + # Dimensions in `dimensions` also appearing in `grid` have same size + # (given by `shape`) as that provided in `grid` if len(shape) != len(dimensions): - raise TypeError("`shape` and `dimensions` must have the " - "same number of entries") + raise ValueError("`shape` and `dimensions` must have the " + "same number of entries") loc_shape = [] for d, s in zip(dimensions, shape): if d in grid.dimensions: size = grid.dimension_map[d] - if size.glb != s: - raise TypeError("Dimension `%s` is given size `%d`, " - "while `grid` says `%s` has size `%d` " - % (d, s, d, size.glb)) + if size.glb != s and s is not None: + raise ValueError("Dimension `%s` is given size `%d`, " + "while `grid` says `%s` has size `%d` " + % (d, s, d, size.glb)) else: loc_shape.append(size.loc) else: @@ -848,7 +941,7 @@ def avg(self, p=None, dims=None): # Pickling support _pickle_kwargs = TensorFunction._pickle_kwargs +\ - ['space_order', 'shape', 'dimensions'] + ['space_order', 'shape_global', 'dimensions'] class TimeFunction(Function): @@ -994,7 +1087,7 @@ def __shape_setup__(cls, **kwargs): raise TypeError("Ambiguity detected: provide either `grid` and `save` " "or just `shape` ") elif shape is None: - shape = list(grid.shape_domain) + shape = list(grid.shape_local) if save is None: shape.insert(cls._time_position, time_order + 1) elif isinstance(save, Buffer): @@ -1143,10 +1236,10 @@ def gridpoints(self): """ The *reference* grid point corresponding to each sparse point. - .. note:: - - When using MPI, this property refers to the *physically* owned - sparse points. + Notes + ----- + When using MPI, this property refers to the *physically* owned + sparse points. """ raise NotImplementedError @@ -1166,7 +1259,7 @@ def _support(self): @property def _dist_datamap(self): """ - Return a mapper ``M : MPI rank -> required sparse data``. + Mapper ``M : MPI rank -> required sparse data``. """ ret = {} for i, s in enumerate(self._support): @@ -1178,12 +1271,11 @@ def _dist_datamap(self): @property def _dist_scatter_mask(self): """ - Return a mask to index into ``self.data``, which creates a new - data array that logically contains N consecutive groups of sparse - data values, where N is the number of MPI ranks. The i-th group - contains the sparse data values accessible by the i-th MPI rank. - Thus, sparse data values along the boundary of two or more MPI - ranks are duplicated. + A mask to index into ``self.data``, which creates a new data array that + logically contains N consecutive groups of sparse data values, where N + is the number of MPI ranks. The i-th group contains the sparse data + values accessible by the i-th MPI rank. Thus, sparse data values along + the boundary of two or more MPI ranks are duplicated. """ dmap = self._dist_datamap mask = np.array(flatten(dmap[i] for i in sorted(dmap)), dtype=int) @@ -1203,10 +1295,10 @@ def _dist_subfunc_scatter_mask(self): @property def _dist_gather_mask(self): """ - Return a mask to index into the ``data`` received upon returning - from ``self._dist_alltoall``. This mask creates a new data array - in which duplicate sparse data values have been discarded. The - resulting data array can thus be used to populate ``self.data``. + A mask to index into the ``data`` received upon returning from + ``self._dist_alltoall``. This mask creates a new data array in which + duplicate sparse data values have been discarded. The resulting data + array can thus be used to populate ``self.data``. """ ret = list(self._dist_scatter_mask) mask = ret[self._sparse_position] @@ -1217,9 +1309,8 @@ def _dist_gather_mask(self): @property def _dist_count(self): """ - Return a 2-tuple of comm-sized iterables, which tells how many sparse - points is this MPI rank expected to send/receive to/from each other - MPI rank. + A 2-tuple of comm-sized iterables, which tells how many sparse points + is this MPI rank expected to send/receive to/from each other MPI rank. """ dmap = self._dist_datamap comm = self.grid.distributor.comm @@ -1230,11 +1321,20 @@ def _dist_count(self): return ssparse, rsparse + @cached_property + def _dist_reorder_mask(self): + """ + An ordering mask that puts ``self._sparse_position`` at the front. + """ + ret = (self._sparse_position,) + ret += tuple(i for i, d in enumerate(self.indices) if d is not self._sparse_dim) + return ret + @property def _dist_alltoall(self): """ - Return the metadata necessary to perform an ``MPI_Alltoallv`` distributing - the sparse data values across the MPI ranks needing them. + The metadata necessary to perform an ``MPI_Alltoallv`` distributing the + sparse data values across the MPI ranks needing them. """ ssparse, rsparse = self._dist_count @@ -1265,12 +1365,17 @@ def _dist_alltoall(self): rshape = list(self.shape) rshape[self._sparse_position] = sum(rsparse) + # May have to swap axes, as `MPI_Alltoallv` expects contiguous data, and + # the sparse dimension may not be the outermost + sshape = [sshape[i] for i in self._dist_reorder_mask] + rshape = [rshape[i] for i in self._dist_reorder_mask] + return sshape, scount, sdisp, rshape, rcount, rdisp @property def _dist_subfunc_alltoall(self): """ - Return the metadata necessary to perform an ``MPI_Alltoallv`` distributing + The metadata necessary to perform an ``MPI_Alltoallv`` distributing self's SubFunction values across the MPI ranks needing them. """ raise NotImplementedError @@ -1498,7 +1603,7 @@ def coordinates(self): @property def coordinates_data(self): - return self.coordinates.data + return self.coordinates.data.view(np.ndarray) @property def _coefficients(self): @@ -1588,7 +1693,7 @@ def _index_matrix(self, offset): # A unique symbol for each indirection index indices = filter_ordered(flatten(index_matrix)) - points = OrderedDict([(p, Symbol(name='ii%d' % i)) + points = OrderedDict([(p, Symbol(name='ii_%s_%d' % (self.name, i))) for i, p in enumerate(indices)]) return index_matrix, points @@ -1633,14 +1738,16 @@ def gridpoints(self): return ret def interpolate(self, expr, offset=0, increment=False, self_subs={}): - """Creates a :class:`sympy.Eq` equation for the interpolation - of an expression onto this sparse point collection. + """Generate equations interpolating an arbitrary expression into ``self``. - :param expr: The expression to interpolate. - :param offset: Additional offset from the boundary for - absorbing boundary conditions. - :param increment: (Optional) if True, perform an increment rather - than an assignment. Defaults to False. + Parameters + ---------- + expr : sympy.Expr + Input expression to interpolate. + offset : int, optional + Additional offset from the boundary. + increment: bool, optional + If True, generate increments (Inc) rather than assignments (Eq). """ variables = list(retrieve_function_carriers(expr)) @@ -1662,12 +1769,16 @@ def interpolate(self, expr, offset=0, increment=False, self_subs={}): return eqns + summands + last def inject(self, field, expr, offset=0): - """Symbol for injection of an expression onto a grid + """Generate equations injecting an arbitrary expression into a field. - :param field: The grid field into which we inject. - :param expr: The expression to inject. - :param offset: Additional offset from the boundary for - absorbing boundary conditions. + Parameters + ---------- + field : Function + Input field into which the injection is performed. + expr : sympy.Expr + Injected expression. + offset : int, optional + Additional offset from the boundary. """ variables = list(retrieve_function_carriers(expr)) + [field] @@ -1681,9 +1792,51 @@ def inject(self, field, expr, offset=0): return eqns + def guard(self, expr=None, offset=0): + """ + Generate a guarded expression, that is expressions that are + evaluated by an Operator only if certain conditions are met. + The introduced condition, here, is that all grid points in the + support of a sparse value must fall within the grid domain (i.e., + *not* on the halo). + + Parameters + ---------- + expr : sympy.Expr, optional + Input expression, from which the guarded expression is derived. + If not specified, defaults to ``self``. + offset : int, optional + Relax the guard condition by introducing a tolerance offset. + """ + _, points = self._index_matrix(offset) + + # Guard through ConditionalDimension + conditions = {} + for d, idx in zip(self.grid.dimensions, self._coordinate_indices): + p = points[idx] + lb = sympy.And(p >= d.symbolic_start - offset, evaluate=False) + ub = sympy.And(p <= d.symbolic_end + offset, evaluate=False) + conditions[p] = sympy.And(lb, ub, evaluate=False) + condition = sympy.And(*conditions.values(), evaluate=False) + cd = ConditionalDimension("%s_g" % self._sparse_dim, self._sparse_dim, + condition=condition) + + if expr is None: + out = self.indexify().xreplace({self._sparse_dim: cd}) + else: + functions = {f for f in retrieve_function_carriers(expr) + if f.is_SparseFunction} + out = indexify(expr).xreplace({f._sparse_dim: cd for f in functions}) + + # Equations for the indirection dimensions + eqns = [Eq(v, k) for k, v in points.items() if v in conditions] + + return out, eqns + @cached_property def _decomposition(self): - return {self._sparse_dim: self._distributor.decomposition[self._sparse_dim]} + mapper = {self._sparse_dim: self._distributor.decomposition[self._sparse_dim]} + return tuple(mapper.get(d) for d in self.dimensions) @property def _dist_subfunc_alltoall(self): @@ -1721,14 +1874,17 @@ def _dist_scatter(self, data=None): comm = distributor.comm mpitype = MPI._typedict[np.dtype(self.dtype).char] - # Pack (reordered) data values so that they can be sent out via an Alltoallv + # Pack sparse data values so that they can be sent out via an Alltoallv data = data[self._dist_scatter_mask] + data = np.ascontiguousarray(np.transpose(data, self._dist_reorder_mask)) # Send out the sparse point values _, scount, sdisp, rshape, rcount, rdisp = self._dist_alltoall scattered = np.empty(shape=rshape, dtype=self.dtype) comm.Alltoallv([data, scount, sdisp, mpitype], [scattered, rcount, rdisp, mpitype]) data = scattered + # Unpack data values so that they follow the expected storage layout + data = np.ascontiguousarray(np.transpose(data, self._dist_reorder_mask)) # Pack (reordered) coordinates so that they can be sent out via an Alltoallv coords = self.coordinates.data._local[self._dist_subfunc_scatter_mask] @@ -1740,7 +1896,7 @@ def _dist_scatter(self, data=None): coords = scattered # Translate global coordinates into local coordinates - coords = coords - np.array(self.grid.origin_domain, dtype=self.dtype) + coords = coords - np.array(self.grid.origin_offset, dtype=self.dtype) return {self: data, self.coordinates: coords} @@ -1753,6 +1909,8 @@ def _dist_gather(self, data): comm = distributor.comm + # Pack sparse data values so that they can be sent out via an Alltoallv + data = np.ascontiguousarray(np.transpose(data, self._dist_reorder_mask)) # Send back the sparse point values sshape, scount, sdisp, _, rcount, rdisp = self._dist_alltoall gathered = np.empty(shape=sshape, dtype=self.dtype) @@ -1760,8 +1918,8 @@ def _dist_gather(self, data): comm.Alltoallv([data, rcount, rdisp, mpitype], [gathered, scount, sdisp, mpitype]) data = gathered - - # Insert back into `self.data` based on the original (expected) data layout + # Unpack data values so that they follow the expected storage layout + data = np.ascontiguousarray(np.transpose(data, self._dist_reorder_mask)) self._data[:] = data[self._dist_gather_mask] # Note: this method "mirrors" `_dist_scatter`: a sparse point that is sent diff --git a/devito/grid.py b/devito/grid.py index ce311f23a3..d8700ee4ea 100644 --- a/devito/grid.py +++ b/devito/grid.py @@ -187,9 +187,9 @@ def spacing_map(self): return dict(zip(self.spacing_symbols, self.spacing)) @property - def origin_domain(self): + def origin_offset(self): """ - Origin of the local (per-process) physical domain. + Offset of the local (per-process) origin from the domain origin. """ grid_origin = [min(i) for i in self.distributor.glb_numb] assert len(grid_origin) == len(self.spacing) @@ -201,7 +201,7 @@ def shape(self): return self._shape @property - def shape_domain(self): + def shape_local(self): """Shape of the local (per-process) physical domain.""" return self._distributor.shape @@ -212,7 +212,7 @@ def dimension_map(self): local size. """ return {d: namedtuple('Size', 'glb loc')(g, l) - for d, g, l in zip(self.dimensions, self.shape, self.shape_domain)} + for d, g, l in zip(self.dimensions, self.shape, self.shape_local)} @property def distributor(self): @@ -241,6 +241,9 @@ def _arg_defaults(self): """ args = ReducerMap() + for k, v in self.dimension_map.items(): + args.update(k._arg_defaults(start=0, size=v.loc)) + if configuration['mpi']: distributor = self.distributor args[distributor._C_comm.name] = distributor._C_comm.value diff --git a/devito/ir/equations/equation.py b/devito/ir/equations/equation.py index 5a39fb3b73..2c9fac556f 100644 --- a/devito/ir/equations/equation.py +++ b/devito/ir/equations/equation.py @@ -230,3 +230,7 @@ def __new__(cls, *args): else: raise ValueError("Cannot construct DummyEq from args=%s" % str(args)) return ClusterizedEq.__new__(cls, obj, ispace=obj.ispace, dspace=obj.dspace) + + # Pickling support + _pickle_args = ['lhs', 'rhs'] + _pickle_kwargs = [] diff --git a/devito/logger.py b/devito/logger.py index a0f10ede4c..31af5dea39 100644 --- a/devito/logger.py +++ b/devito/logger.py @@ -104,8 +104,8 @@ def set_log_noperf(): logger.setLevel(WARNING) -configuration.add('log_level', 'INFO', list(logger_registry), - lambda i: set_log_level(i)) +configuration.add('log-level', 'INFO', list(logger_registry), + lambda i: set_log_level(i), False) class silencio(object): @@ -120,10 +120,10 @@ def __init__(self, log_level='WARNING'): def __call__(self, func, *args, **kwargs): @wraps(func) def wrapper(*args, **kwargs): - previous = configuration['log_level'] - configuration['log_level'] = self.log_level + previous = configuration['log-level'] + configuration['log-level'] = self.log_level result = func(*args, **kwargs) - configuration['log_level'] = previous + configuration['log-level'] = previous return result return wrapper diff --git a/devito/mpi/distributed.py b/devito/mpi/distributed.py index 288b1a0b63..7206b32051 100644 --- a/devito/mpi/distributed.py +++ b/devito/mpi/distributed.py @@ -12,7 +12,7 @@ from devito.data import Decomposition from devito.parameters import configuration -from devito.types import LEFT, RIGHT +from devito.types import LEFT, RIGHT, CompositeObject, Object from devito.tools import EnrichedTuple, as_tuple, is_integer @@ -26,8 +26,10 @@ except ImportError: class MPI(object): COMM_NULL = None + + @classmethod def Is_initialized(): - return False + return False __all__ = ['Distributor', 'SparseDistributor', 'MPI'] @@ -66,6 +68,15 @@ def mycoords(self): """The coordinates of the calling MPI rank in the Distributor topology.""" return + @abstractmethod + def nprocs(self): + """A shortcut for the number of processes in the MPI communicator.""" + return + + @property + def is_parallel(self): + return self.nprocs > 1 + @cached_property def glb_numb(self): """The global indices owned by the calling MPI rank.""" @@ -300,26 +311,12 @@ def neighbours(self): @cached_property def _C_comm(self): - """ - A :class:`Object` wrapping an MPI communicator. - - Extracted from: :: - - https://github.com/mpi4py/mpi4py/blob/master/demo/wrap-ctypes/helloworld.py - """ - from devito.types import Object - if MPI._sizeof(self._comm) == sizeof(c_int): - ctype = type('MPI_Comm', (c_int,), {}) - else: - ctype = type('MPI_Comm', (c_void_p,), {}) - comm_ptr = MPI._addressof(self._comm) - comm_val = ctype.from_address(comm_ptr) - return Object(name='comm', dtype=ctype, value=comm_val) + """An :class:`Object` wrapping an MPI communicator.""" + return MPICommObject(self.comm) @cached_property def _C_neighbours(self): - """A ctypes Struct to access the neighborhood of a given rank.""" - from devito.types import CompositeObject + """A :class:`ctypes.Struct` to access the neighborhood of a given rank.""" entries = list(product(self.dimensions, [LEFT, RIGHT])) fields = [('%s%s' % (d, i), c_int) for d, i in entries] obj = CompositeObject('nb', 'neighbours', Structure, fields) @@ -399,6 +396,28 @@ def nprocs(self): return self.distributor.nprocs +class MPICommObject(Object): + + name = 'comm' + + # See https://github.com/mpi4py/mpi4py/blob/master/demo/wrap-ctypes/helloworld.py + if MPI._sizeof(MPI.Comm) == sizeof(c_int): + dtype = type('MPI_Comm', (c_int,), {}) + else: + dtype = type('MPI_Comm', (c_void_p,), {}) + + def __init__(self, comm=None): + if comm is None: + # Should only end up here upon unpickling + comm = MPI.COMM_WORLD + comm_ptr = MPI._addressof(comm) + comm_val = self.dtype.from_address(comm_ptr) + self.value = comm_val + + # Pickling support + _pickle_args = [] + + def compute_dims(nprocs, ndim): # We don't do anything clever here. In fact, we do something very basic -- # we just try to distribute `nprocs` evenly over the number of dimensions, diff --git a/devito/mpi/routines.py b/devito/mpi/routines.py index 935f57949e..74263c67df 100644 --- a/devito/mpi/routines.py +++ b/devito/mpi/routines.py @@ -85,12 +85,9 @@ def sendrecv(f, fixed): # the domain boundary, where the sender is actually MPI.PROC_NULL scatter = Conditional(CondNe(fromrank, Macro('MPI_PROC_NULL')), scatter) - MPI_Status = type('MPI_Status', (c_void_p,), {}) - srecv = LocalObject(name='srecv', dtype=MPI_Status) - - MPI_Request = type('MPI_Request', (c_void_p,), {}) - rrecv = LocalObject(name='rrecv', dtype=MPI_Request) - rsend = LocalObject(name='rsend', dtype=MPI_Request) + srecv = MPIStatusObject(name='srecv') + rrecv = MPIRequestObject(name='rrecv') + rsend = MPIRequestObject(name='rsend') count = reduce(mul, bufs.shape, 1) recv = Call('MPI_Irecv', [bufs, count, Macro(numpy_to_mpitypes(f.dtype)), @@ -161,3 +158,25 @@ def update_halo(f, fixed): parameters = ([f] + masks + [comm, nb] + list(fixed.values()) + [d.symbolic_size for d in f.dimensions]) return Callable('halo_exchange_%s' % f.name, iet, 'void', parameters, ('static',)) + + +class MPIStatusObject(LocalObject): + + dtype = type('MPI_Status', (c_void_p,), {}) + + def __init__(self, name): + self.name = name + + # Pickling support + _pickle_args = ['name'] + + +class MPIRequestObject(LocalObject): + + dtype = type('MPI_Request', (c_void_p,), {}) + + def __init__(self, name): + self.name = name + + # Pickling support + _pickle_args = ['name'] diff --git a/devito/operator.py b/devito/operator.py index 9b2d22f3d0..f1f92923ee 100644 --- a/devito/operator.py +++ b/devito/operator.py @@ -12,7 +12,7 @@ from devito.dle import transform from devito.dse import rewrite from devito.exceptions import InvalidOperator -from devito.logger import info, perf +from devito.logger import info, perf, warning from devito.ir.equations import LoweredEq from devito.ir.clusters import clusterize from devito.ir.iet import (Callable, List, MetaCall, iet_build, iet_insert_C_decls, @@ -171,9 +171,11 @@ def _prepare_arguments(self, **kwargs): args = self._autotune(args) # Check all user-provided keywords are known to the Operator - for k, v in kwargs.items(): - if k not in self._known_arguments: - raise ValueError("Unrecognized argument %s=%s passed to `apply`" % (k, v)) + if not configuration['ignore-unknowns']: + for k, v in kwargs.items(): + if k not in self._known_arguments: + raise ValueError("Unrecognized argument %s=%s passed to `apply`" + % (k, v)) return args @@ -316,6 +318,7 @@ def _build_casts(self, iet): def __getstate__(self): if self._lib: state = dict(self.__dict__) + state.pop('_soname') # The compiled shared-object will be pickled; upon unpickling, it # will be restored into a potentially different temporary directory, # so the entire process during which the shared-object is loaded and @@ -329,10 +332,26 @@ def __getstate__(self): return self.__dict__ def __setstate__(self, state): - binary = state.pop('binary') + soname = state.pop('_soname', None) + binary = state.pop('binary', None) for k, v in state.items(): setattr(self, k, v) - save(self._soname, binary, self._compiler) + # If the `sonames` don't match, there *might* be a hidden bug as the + # unpickled Operator might be generating code that differs from that + # generated by the pickled Operator. For example, a stupid bug that we + # had to fix was due to rebuilding SymPy expressions which weren't + # automatically getting the flag `evaluate=False`, thus producing x+2 + # on the unpickler instead of x+1+1). However, different `sonames` + # doesn't necessarily means there's a bug: if the unpickler and the + # pickler are two distinct processes and the unpickler runs with a + # different `configuration` dictionary, then the `sonames` might indeed + # be different, depending on which entries in `configuration` differ. + if soname is not None: + if soname != self._soname: + warning("The pickled and unpickled Operators have different .sonames; " + "this might be a bug, or simply a harmless difference in " + "`configuration`. You may check they produce the same code.") + save(self._soname, binary, self._compiler) class OperatorRunnable(Operator): diff --git a/devito/ops/__init__.py b/devito/ops/__init__.py new file mode 100644 index 0000000000..a022095401 --- /dev/null +++ b/devito/ops/__init__.py @@ -0,0 +1,25 @@ +""" +The ``ops`` Devito backend uses the OPS library to generate, +JIT-compile, and run kernels on multiple architectures. +""" + +from devito.dle import (BasicRewriter, AdvancedRewriter, AdvancedRewriterSafeMath, + SpeculativeRewriter, init_dle) +from devito.parameters import Parameters, add_sub_configuration + +ops_configuration = Parameters('ops') +env_vars_mapper = {} +add_sub_configuration(ops_configuration, env_vars_mapper) + +# Initialize the DLE +modes = {'basic': BasicRewriter, + 'advanced': AdvancedRewriter, + 'advanced-safemath': AdvancedRewriterSafeMath, + 'speculative': SpeculativeRewriter} +init_dle(modes) + +# The following used by backends.backendSelector +from devito.function import * # noqa +from devito.grid import Grid # noqa +from devito.ops.operator import Operator # noqa +from devito.types import CacheManager # noqa diff --git a/devito/ops/operator.py b/devito/ops/operator.py new file mode 100644 index 0000000000..79735e88d3 --- /dev/null +++ b/devito/ops/operator.py @@ -0,0 +1,12 @@ +from devito.operator import OperatorRunnable + +__all__ = ['Operator'] + + +class Operator(OperatorRunnable): + """ + A special :class:`OperatorCore` to JIT-compile and run operators through OPS. + """ + + def _specialize_iet(self, iet, **kwargs): + raise NotImplementedError diff --git a/devito/parameters.py b/devito/parameters.py index ef8f0be163..be5dbf8489 100644 --- a/devito/parameters.py +++ b/devito/parameters.py @@ -27,6 +27,7 @@ def __init__(self, name=None, **kwargs): self._name = name self._accepted = {} self._defaults = {} + self._impact_jit = {} self._update_functions = {} if kwargs is not None: for key, value in kwargs.items(): @@ -66,7 +67,7 @@ def update(self, key, value): """ super(Parameters, self).__setitem__(key, value) - def add(self, key, value, accepted=None, callback=None): + def add(self, key, value, accepted=None, callback=None, impacts_jit=True): """ Add a new parameter ``key`` with default value ``value``. @@ -74,10 +75,15 @@ def add(self, key, value, accepted=None, callback=None): If provided, make sure ``callback`` is executed when the value of ``key`` changes. + + If ``impacts_jit`` is False (defaults to True), then it can be assumed + that the parameter doesn't affect code generation, so it can be excluded + from the construction of the hash key. """ super(Parameters, self).__setitem__(key, value) self._accepted[key] = accepted self._defaults[key] = value + self._impact_jit[key] = impacts_jit if callable(callback): self._update_functions[key] = callback @@ -94,10 +100,9 @@ def name(self): return self._name def _signature_items(self): - # Note: we are discarding some vars that do not affect the c level + # Note: we are discarding some vars that do not affect the C level # code in order to avoid recompiling when such vars are modified - discard = ['profiling', 'autotuning', 'log_level', 'first_touch'] - items = sorted((k, v) for k, v in self.items() if k not in discard) + items = sorted((k, v) for k, v in self.items() if self._impact_jit[k]) return tuple(str(items)) + tuple(str(sorted(self.backend.items()))) @@ -110,13 +115,14 @@ def _signature_items(self): 'DEVITO_DEVELOP': 'develop-mode', 'DEVITO_DSE': 'dse', 'DEVITO_DLE': 'dle', - 'DEVITO_DLE_OPTIONS': 'dle_options', + 'DEVITO_DLE_OPTIONS': 'dle-options', 'DEVITO_OPENMP': 'openmp', 'DEVITO_MPI': 'mpi', 'DEVITO_AUTOTUNING': 'autotuning', - 'DEVITO_LOGGING': 'log_level', - 'DEVITO_FIRST_TOUCH': 'first_touch', - 'DEVITO_DEBUG_COMPILER': 'debug_compiler', + 'DEVITO_LOGGING': 'log-level', + 'DEVITO_FIRST_TOUCH': 'first-touch', + 'DEVITO_DEBUG_COMPILER': 'debug-compiler', + 'DEVITO_IGNORE_UNKNOWN_PARAMS': 'ignore-unknowns' } diff --git a/devito/profiling.py b/devito/profiling.py index 617770a7d6..81bb68d692 100644 --- a/devito/profiling.py +++ b/devito/profiling.py @@ -204,12 +204,16 @@ def __init__(self, name, sections): [(i, c_double) for i in sections]) def reset(self): - for i in self.pfields: + for i, _ in self.pfields: setattr(self.value._obj, i, 0.0) return self.value + @property + def sections(self): + return [i for i, _ in self.pfields] + # Pickling support - _pickle_args = ['name', 'pfields'] + _pickle_args = ['name', 'sections'] _pickle_kwargs = [] @@ -247,7 +251,7 @@ def create_profile(name): """ Create a new :class:`Profiler`. """ - if configuration['log_level'] == 'DEBUG': + if configuration['log-level'] == 'DEBUG': # Enforce performance profiling in DEBUG mode level = 'advanced' else: @@ -269,7 +273,7 @@ def create_profile(name): 'advanced': AdvancedProfiler, 'advisor': AdvisorProfiler } -configuration.add('profiling', 'basic', list(profiler_registry)) +configuration.add('profiling', 'basic', list(profiler_registry), impacts_jit=False) def locate_intel_advisor(): diff --git a/devito/symbolics/extended_sympy.py b/devito/symbolics/extended_sympy.py index 51ef9231ff..5f1cfa25e9 100644 --- a/devito/symbolics/extended_sympy.py +++ b/devito/symbolics/extended_sympy.py @@ -43,6 +43,9 @@ def xreplace(self, rule): return self.func(*args, evaluate=False) return self + def evalf(self, *args, **kwargs): + return self + class Eq(sympy.Eq, FrozenExpr): @@ -77,15 +80,18 @@ def canonical(self): class Mul(sympy.Mul, FrozenExpr): - pass + def __new__(cls, *args, **kwargs): + return sympy.Mul.__new__(cls, *args, evaluate=False) class Add(sympy.Add, FrozenExpr): - pass + def __new__(cls, *args, **kwargs): + return sympy.Add.__new__(cls, *args, evaluate=False) class Pow(sympy.Pow, FrozenExpr): - pass + def __new__(cls, *args, **kwargs): + return sympy.Pow.__new__(cls, *args, evaluate=False) class IntDiv(sympy.Expr): diff --git a/devito/types.py b/devito/types.py index 01a47da979..1939af197d 100644 --- a/devito/types.py +++ b/devito/types.py @@ -11,13 +11,10 @@ import numpy as np import sympy -from devito.parameters import configuration from devito.tools import ArgProvider, EnrichedTuple, Pickable, Tag, ctypes_to_C __all__ = ['Symbol', 'Indexed'] -configuration.add('first_touch', 0, [0, 1], lambda i: bool(i)) - # This cache stores a reference to each created data object # so that we may re-create equivalent symbols during symbolic # manipulation with the correct shapes, pointers, etc. @@ -436,11 +433,11 @@ def symbolic_shape(self): halo, and domain regions. While halo and padding are known quantities (integers), the domain size is represented by a symbol. """ - halo_sizes = [sympy.Add(*i, evaluate=False) for i in self._extent_halo] - padding_sizes = [sympy.Add(*i, evaluate=False) for i in self._extent_padding] + from devito.symbolics.extended_sympy import Add + halo_sizes = [Add(*i) for i in self._extent_halo] + pad_sizes = [Add(*i) for i in self._extent_padding] domain_sizes = [i.symbolic_size for i in self.indices] - ret = tuple(sympy.Add(i, j, k, evaluate=False) - for i, j, k in zip(domain_sizes, halo_sizes, padding_sizes)) + ret = tuple(Add(i, j, k) for i, j, k in zip(domain_sizes, halo_sizes, pad_sizes)) return EnrichedTuple(*ret, getters=self.dimensions) @property @@ -792,7 +789,7 @@ def __init__(self, name, pname, ptype, pfields, value=None): @property def pfields(self): - return tuple(i for i, _ in self.dtype._type_._fields_) + return tuple(self.dtype._type_._fields_) @property def ptype(self): diff --git a/devito/yask/__init__.py b/devito/yask/__init__.py index 2abc5a2b36..b09c2b5691 100644 --- a/devito/yask/__init__.py +++ b/devito/yask/__init__.py @@ -72,11 +72,11 @@ def __init__(self, *args, **kwargs): yask_configuration = Parameters('yask') yask_configuration.add('compiler', YaskCompiler()) callback = lambda i: eval(i) if i else () -yask_configuration.add('folding', (), callback=callback) -yask_configuration.add('blockshape', (), callback=callback) -yask_configuration.add('clustering', (), callback=callback) -yask_configuration.add('options', None) -yask_configuration.add('dump', None) +yask_configuration.add('folding', (), callback=callback, impacts_jit=False) +yask_configuration.add('blockshape', (), callback=callback, impacts_jit=False) +yask_configuration.add('clustering', (), callback=callback, impacts_jit=False) +yask_configuration.add('options', None, impacts_jit=False) +yask_configuration.add('dump', None, impacts_jit=False) env_vars_mapper = { 'DEVITO_YASK_FOLDING': 'folding', diff --git a/devito/yask/function.py b/devito/yask/function.py index 757fca09eb..00f55b8dc5 100644 --- a/devito/yask/function.py +++ b/devito/yask/function.py @@ -155,7 +155,7 @@ def data_with_halo(self): @cached_property @_allocate_memory - def data_allocated(self): + def _data_allocated(self): return Data(self._data.grid, self.shape_allocated, self.indices, self.dtype) def _arg_defaults(self, alias=None): diff --git a/devito/yask/operator.py b/devito/yask/operator.py index 419db9bb4e..759e295538 100644 --- a/devito/yask/operator.py +++ b/devito/yask/operator.py @@ -28,8 +28,7 @@ class Operator(OperatorRunnable): A special :class:`OperatorCore` to JIT-compile and run operators through YASK. """ - _default_headers = OperatorRunnable._default_headers - _default_headers += ['#define restrict __restrict'] + _default_headers = OperatorRunnable._default_headers + ['#define restrict __restrict'] _default_includes = OperatorRunnable._default_includes + ['yask_kernel_api.hpp'] def __init__(self, expressions, **kwargs): @@ -164,7 +163,7 @@ def apply(self, **kwargs): return self._profile_output(args) def __getstate__(self): - state = super(Operator, self).__getstate__() + state = dict(super(Operator, self).__getstate__()) # A YASK solution object needs to be recreated upon unpickling. Steps: # 1) upon pickling: serialise all files generated by this Operator via YASK # 2) upon unpickling: deserialise and explicitly recreate the YASK solution diff --git a/docker-compose.yml b/docker-compose.yml index af958f6ef8..9ea9759e53 100644 --- a/docker-compose.yml +++ b/docker-compose.yml @@ -9,7 +9,7 @@ services: dockerfile: docker/Dockerfile environment: - DEVITO_ARCH: gcc-4.9 + DEVITO_ARCH: gcc DEVITO_OPENMP: 0 volumes: diff --git a/docker/Dockerfile b/docker/Dockerfile index 4df1990b31..0f0e4e5012 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -28,7 +28,7 @@ RUN chmod +x \ WORKDIR /app -ENV DEVITO_ARCH="gcc-4.9" +ENV DEVITO_ARCH="gcc" ENV DEVITO_OPENMP="0" EXPOSE 8888 diff --git a/docker/entrypoint.sh b/docker/entrypoint.sh index 51eddd40f5..71c4ed9d65 100644 --- a/docker/entrypoint.sh +++ b/docker/entrypoint.sh @@ -2,4 +2,7 @@ find /app -type f -name '*.pyc' -delete +export PATH=/venv/bin:$PATH +export PYTHONPATH=$PYTHONPATH:/app + exec "$@" diff --git a/environment.yml b/environment.yml index 890e49f6d9..3c2e003721 100644 --- a/environment.yml +++ b/environment.yml @@ -8,7 +8,6 @@ dependencies: - sympy==1.2 - cgen - scipy - - mpi4py - matplotlib - pytest - flake8>=2.1.0 diff --git a/examples/cfd/05_laplace.ipynb b/examples/cfd/05_laplace.ipynb index 5f83619847..28550c141b 100644 --- a/examples/cfd/05_laplace.ipynb +++ b/examples/cfd/05_laplace.ipynb @@ -253,7 +253,7 @@ "\n", "# Silence the runtime performance logging\n", "from devito import configuration\n", - "configuration['log_level'] = 'ERROR'\n", + "configuration['log-level'] = 'ERROR'\n", "\n", "# Initialise the two buffer fields\n", "p.data[:] = 0.\n", diff --git a/examples/cfd/06_poisson.ipynb b/examples/cfd/06_poisson.ipynb index d3fbedc356..54be602bb5 100644 --- a/examples/cfd/06_poisson.ipynb +++ b/examples/cfd/06_poisson.ipynb @@ -119,7 +119,7 @@ "from devito import Grid, Function, TimeFunction, Operator, configuration, Eq, solve\n", "\n", "# Silence the runtime performance logging\n", - "configuration['log_level'] = 'ERROR'\n", + "configuration['log-level'] = 'ERROR'\n", "\n", "# Now with Devito we will turn `p` into `TimeFunction` object \n", "# to make all the buffer switching implicit\n", @@ -233,7 +233,7 @@ "\n", "# We can even switch performance logging back on,\n", "# since we only require a single kernel invocation.\n", - "configuration['log_level'] = 'INFO'\n", + "configuration['log-level'] = 'INFO'\n", " \n", "# Create and execute the operator for a number of timesteps\n", "op = Operator([eq_stencil] + bc)\n", diff --git a/examples/cfd/tools.py b/examples/cfd/tools.py index f194c13283..0e6cef093a 100644 --- a/examples/cfd/tools.py +++ b/examples/cfd/tools.py @@ -16,7 +16,7 @@ def plot_field(field, xmax=2., ymax=2., zmax=None, view=None, linewidth=0): y_coord = np.linspace(0, ymax, field.shape[1]) fig = pyplot.figure(figsize=(11, 7), dpi=100) ax = fig.gca(projection='3d') - X, Y = np.meshgrid(x_coord, y_coord) + X, Y = np.meshgrid(x_coord, y_coord, indexing='ij') ax.plot_surface(X, Y, field[:], cmap=cm.viridis, rstride=1, cstride=1, linewidth=linewidth, antialiased=False) diff --git a/examples/compiler/01_iet.ipynb b/examples/compiler/01_iet.ipynb index dd3f2b3698..2ffac51f5f 100644 --- a/examples/compiler/01_iet.ipynb +++ b/examples/compiler/01_iet.ipynb @@ -49,7 +49,9 @@ "source": [ "Here, we just defined a time-varying grid with 3x3 points in each of the space `Dimension`s _x_ and _y_. `u.data` gives us access to the actual data values. In particular, `u.data[0, :, :]` holds the data values at timestep `t=0`, whereas `u.data[1, :, :]` holds the data values at timestep `t=1`.\n", "\n", - "We now create an `Operator` that increments by 1 all points in the computational domain." + "We now create an `Operator` that increments by 1 all points in the computational domain.", + "\n", + "Note: User should not change to `configuration['openmp'] = 1` as the IET structure changes with OpenMP enabled." ] }, { diff --git a/examples/compiler/03_data_regions.ipynb b/examples/compiler/03_data_regions.ipynb new file mode 100644 index 0000000000..189ecf6c3f --- /dev/null +++ b/examples/compiler/03_data_regions.ipynb @@ -0,0 +1,555 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Domain, Halo and Padding regions" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this tutorial we will learn about data regions and how these impact the Operator construction. We will use the simple time marching example shown in the [01_iet](https://github.com/opesci/devito/blob/master/examples/compiler/01_iet.ipynb) tutorial. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from devito import Eq, Grid, TimeFunction, Operator\n", + "\n", + "grid = Grid(shape=(3, 3))\n", + "u = TimeFunction(name='u', grid=grid)\n", + "u.data[:] = 1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At this point, we have a time-varying 3x3 grid filled with `1's`. Below, we can see the `domain` data values:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[1. 1. 1.]\n", + " [1. 1. 1.]\n", + " [1. 1. 1.]]\n", + "\n", + " [[1. 1. 1.]\n", + " [1. 1. 1.]\n", + " [1. 1. 1.]]]\n" + ] + } + ], + "source": [ + "print(u.data)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now create an `Operator` that, at each timestep, increments by `2` all points in the computational domain." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from devito import configuration\n", + "configuration['openmp'] = 0\n", + "\n", + "eq = Eq(u.forward, u+2)\n", + "op = Operator(eq)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can print `op` to get the generated code." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "#define _POSIX_C_SOURCE 200809L\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(float *restrict u_vec, const int time_M, const int time_m, struct profiler* timers, const int x_M, const int x_m, const int x_size, const int y_M, const int y_m, const int y_size)\n", + "{\n", + " float (*restrict u)[x_size + 1 + 1][y_size + 1 + 1] __attribute__((aligned(64))) = (float (*)[x_size + 1 + 1][y_size + 1 + 1]) u_vec;\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + " for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))\n", + " {\n", + " struct timeval start_section0, end_section0;\n", + " gettimeofday(&start_section0, NULL);\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " #pragma omp simd\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " u[t1][x + 1][y + 1] = u[t0][x + 1][y + 1] + 2;\n", + " }\n", + " }\n", + " gettimeofday(&end_section0, NULL);\n", + " timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;\n", + " }\n", + " return 0;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "print(op)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When we take a look at the constructed expression, `u[t1][x + 1][y + 1] = u[t0][x + 1][y + 1] + 2`, we see several `+1` were added to the `u`'s spatial indexes." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is because the domain region is actually surrounded by 'ghost' points, which can be accessed via a stencil when iterating in proximity of the domain boundary. The ghost points define the `halo` region. The halo region can be accessed through the `data_with_halo` data accessor; as we see below, the halo points correspond to the zeros surrounding the domain region." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]]\n" + ] + } + ], + "source": [ + "print(u.data_with_halo)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By adding the '+1' offsets, the Devito compiler ensures the array accesses are logically aligned to the equation’s physical domain. For instance, the `TimeFunction` `u(t, x, y)` used in the example above has one point on each side of the `x` and `y` halo regions; if the user writes an expression including `u(t, x, y)` and `u(t, x + 2, y + 2)`, the compiler will ultimately generate `u[t, x + 1, y + 1]` and `u[t, x + 3, y + 3]`. When `x = y = 0`, therefore, the values `u[t, 1, 1]` and `u[t, 3, 3]` are fetched, representing the first and third points in the physical domain. \n", + "\n", + "By default, the halo region has `space_order` points on each side of the space dimensions. Sometimes, these points may be unnecessary, or, depending on the partial differential equation being approximated, extra points may be needed." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[1. 1. 1.]\n", + " [1. 1. 1.]\n", + " [1. 1. 1.]]\n", + "\n", + " [[1. 1. 1.]\n", + " [1. 1. 1.]\n", + " [1. 1. 1.]]]\n" + ] + } + ], + "source": [ + "u0 = TimeFunction(name='u0', grid=grid, space_order=0)\n", + "u0.data[:] = 1\n", + "print(u0.data_with_halo)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]]]\n" + ] + } + ], + "source": [ + "u2 = TimeFunction(name='u2', grid=grid, space_order=2)\n", + "u2.data[:] = 1\n", + "print(u2.data_with_halo)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "One can also pass a 3-tuple `(o, lp, rp)` instead of a single integer representing the discretization order. Here, `o` is the discretization order, while `lp` and `rp` indicate how many points are expected on left (lp) and right (rp) sides of a point of interest." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "u_new = TimeFunction(name='u_new', grid=grid, space_order=(4, 3, 1))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]]]\n" + ] + } + ], + "source": [ + "u_new.data[:] = 1\n", + "print(u_new.data_with_halo)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's have a look at the Operator generated code when using `u_new`." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "#define _POSIX_C_SOURCE 200809L\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(float *restrict u_new_vec, const int time_M, const int time_m, struct profiler* timers, const int x_M, const int x_m, const int x_size, const int y_M, const int y_m, const int y_size)\n", + "{\n", + " float (*restrict u_new)[x_size + 1 + 3][y_size + 1 + 3] __attribute__((aligned(64))) = (float (*)[x_size + 1 + 3][y_size + 1 + 3]) u_new_vec;\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + " for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))\n", + " {\n", + " struct timeval start_section0, end_section0;\n", + " gettimeofday(&start_section0, NULL);\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " #pragma omp simd\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " u_new[t1][x + 3][y + 3] = u_new[t0][x + 3][y + 3] + 2;\n", + " }\n", + " }\n", + " gettimeofday(&end_section0, NULL);\n", + " timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;\n", + " }\n", + " return 0;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "equation = Eq(u_new.forward, u_new + 2)\n", + "op = Operator(equation)\n", + "print(op)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And finally, let's run it, to convince ourselves that only the domain region values will be incremented at each timestep." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` run in 0.00 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 5. 5. 5. 0.]\n", + " [0. 0. 0. 5. 5. 5. 0.]\n", + " [0. 0. 0. 5. 5. 5. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 7. 7. 7. 0.]\n", + " [0. 0. 0. 7. 7. 7. 0.]\n", + " [0. 0. 0. 7. 7. 7. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0.]]]\n" + ] + } + ], + "source": [ + "op.apply(time_M=2)\n", + "print(u_new.data_with_halo)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The halo region, in turn, is surrounded by the `padding` region, which can be used for data alignment. By default, there is no padding. This can be changed by passing a suitable value to 'padding', as shown below:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "#define _POSIX_C_SOURCE 200809L\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(float *restrict u_pad_vec, const int time_M, const int time_m, struct profiler* timers, const int x_M, const int x_m, const int x_size, const int y_M, const int y_m, const int y_size)\n", + "{\n", + " float (*restrict u_pad)[x_size + 2 + 2 + 2 + 2][y_size + 2 + 2 + 2 + 2] __attribute__((aligned(64))) = (float (*)[x_size + 2 + 2 + 2 + 2][y_size + 2 + 2 + 2 + 2]) u_pad_vec;\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + " for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))\n", + " {\n", + " struct timeval start_section0, end_section0;\n", + " gettimeofday(&start_section0, NULL);\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " #pragma omp simd\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " u_pad[t1][x + 4][y + 4] = u_pad[t0][x + 4][y + 4] + 2;\n", + " }\n", + " }\n", + " gettimeofday(&end_section0, NULL);\n", + " timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;\n", + " }\n", + " return 0;\n", + "}\n", + "\n" + ] + } + ], + "source": [ + "u_pad = TimeFunction(name='u_pad', grid=grid, space_order=2, padding=(0,2,2))\n", + "u_pad.data_with_halo[:] = 1\n", + "u_pad.data[:] = 2\n", + "equation = Eq(u_pad.forward, u_pad + 2)\n", + "op = Operator(equation)\n", + "print(op)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that `+4` was added to each space dimension index. Of this '+4', '+2' is due to the halo (`space_order=2`), and another +2 to the padding.\n", + "Although in practice not very useful, with the (private) `_data_allocated` accessor one can see the entire domain+halo+padding region." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 2. 2. 2. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 2. 2. 2. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 2. 2. 2. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 2. 2. 2. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 2. 2. 2. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 2. 2. 2. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]]\n" + ] + } + ], + "source": [ + "print(u_pad._data_allocated)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Above, the domain is filled with 2's, the halo is filled with 1's, and the padding is filled with 0's." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py index 337f5619f9..284f02f243 100644 --- a/examples/seismic/acoustic/acoustic_example.py +++ b/examples/seismic/acoustic/acoustic_example.py @@ -2,10 +2,9 @@ from argparse import ArgumentParser from devito.logger import info -from devito import Constant, Function +from devito import Constant, Function, smooth from examples.seismic.acoustic import AcousticWaveSolver from examples.seismic import demo_model, TimeAxis, RickerSource, Receiver -from examples.seismic.utils import smooth def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), @@ -30,7 +29,8 @@ def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), rec = Receiver(name='rec', grid=model.grid, time_range=time_range, npoint=nrec) rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec) if len(shape) > 1: - rec.coordinates.data[:, 1:] = src.coordinates.data[0, 1:] + rec.coordinates.data[:, 1] = np.array(model.domain_size)[1] * .5 + rec.coordinates.data[:, -1] = model.origin[-1] + 2 * spacing[-1] # Create solver object to provide relevant operators solver = AcousticWaveSolver(model, source=src, receiver=rec, kernel=kernel, diff --git a/examples/seismic/source.py b/examples/seismic/source.py index fd6bee1d9e..158d11f8ec 100644 --- a/examples/seismic/source.py +++ b/examples/seismic/source.py @@ -93,7 +93,7 @@ def __new__(cls, **kwargs): time_order = kwargs.pop('time_order', 2) p_dim = kwargs.pop('dimension', Dimension(name='p_%s' % name)) - coordinates = kwargs.pop('coordinates', None) + coordinates = kwargs.pop('coordinates', kwargs.pop('coordinates_data', None)) # Either `npoint` or `coordinates` must be provided npoint = kwargs.pop('npoint', None) if npoint is None: @@ -117,10 +117,6 @@ def __new__(cls, **kwargs): return obj - def __init__(self, *args, **kwargs): - if not self._cached(): - super(PointSource, self).__init__(*args, **kwargs) - @cached_property def time_values(self): return self._time_range.time_values diff --git a/examples/seismic/tutorials/02_rtm.ipynb b/examples/seismic/tutorials/02_rtm.ipynb index e863bb3d2b..3a0bafaac5 100644 --- a/examples/seismic/tutorials/02_rtm.ipynb +++ b/examples/seismic/tutorials/02_rtm.ipynb @@ -44,7 +44,7 @@ "%matplotlib inline\n", "\n", "from devito import configuration\n", - "configuration['log_level'] = 'WARNING'" + "configuration['log-level'] = 'WARNING'" ] }, { diff --git a/examples/seismic/tutorials/03_fwi.ipynb b/examples/seismic/tutorials/03_fwi.ipynb index 9f46c84e9a..098073f27e 100644 --- a/examples/seismic/tutorials/03_fwi.ipynb +++ b/examples/seismic/tutorials/03_fwi.ipynb @@ -37,7 +37,7 @@ "%matplotlib inline\n", "\n", "from devito import configuration\n", - "configuration['log_level'] = 'WARNING'" + "configuration['log-level'] = 'WARNING'" ] }, { diff --git a/examples/seismic/tutorials/05_skimage_tv.ipynb b/examples/seismic/tutorials/05_skimage_tv.ipynb index 352a4f27db..5442957cb0 100644 --- a/examples/seismic/tutorials/05_skimage_tv.ipynb +++ b/examples/seismic/tutorials/05_skimage_tv.ipynb @@ -276,7 +276,7 @@ "# Change to the WARNING log level to reduce log output\n", "# as compared to the default DEBUG\n", "from devito import configuration\n", - "configuration['log_level'] = 'WARNING'\n", + "configuration['log-level'] = 'WARNING'\n", "\n", "# Set up inversion parameters.\n", "param = {'t0': 0.,\n", diff --git a/examples/seismic/utils.py b/examples/seismic/utils.py index 25966b38da..898b4d7a68 100644 --- a/examples/seismic/utils.py +++ b/examples/seismic/utils.py @@ -1,20 +1,6 @@ from scipy import ndimage -from devito import Operator, Eq -__all__ = ['smooth', 'scipy_smooth'] - - -# Velocity models -def smooth(dest, f): - """ - Run an n-point moving average kernel over ``f`` and store the result - into ``dest``. The average is computed along the innermost ``f`` dimension. - """ - if f.is_Constant: - # Return a scaled version of the input if it's a Constant - dest.data[:] = .9 * f.data - else: - Operator(Eq(dest, f.avg(dims=f.dimensions[-1])), name='smoother').apply() +__all__ = ['scipy_smooth'] def scipy_smooth(img, sigma=5): diff --git a/requirements-optional.txt b/requirements-optional.txt new file mode 100644 index 0000000000..c501e9666d --- /dev/null +++ b/requirements-optional.txt @@ -0,0 +1,2 @@ +mpi4py +matplotlib \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index e74e898b75..6c459125fe 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,6 @@ numpy==1.14 sympy>=1.2 scipy -matplotlib pytest-runner flake8>=2.1.0 jedi diff --git a/setup.py b/setup.py index c4b93a9962..32ec38c625 100644 --- a/setup.py +++ b/setup.py @@ -6,6 +6,9 @@ with open('requirements.txt') as f: required = f.read().splitlines() +with open('requirements-optional.txt') as f: + optionals = f.read().splitlines() + reqs = [] links = [] for ir in required: @@ -15,8 +18,14 @@ else: reqs += [ir] -if os.environ.get('INSTALL_MPI') == '1': - reqs += ['mpi4py'] +opt_reqs = [] +opt_links = [] +for ir in optionals: + if ir[0:3] == 'git': + opt_links += [ir + '#egg=' + ir.split('/')[-1] + '-0'] + opt_reqs += [ir.split('/')[-1]] + else: + opt_reqs += [ir] setup(name='devito', version=versioneer.get_version(), @@ -29,11 +38,12 @@ equations defined in SymPy to create and execute highly optimised Finite Difference kernels on multiple computer platforms.""", - url='http://www.opesci.org/devito', + url='http://www.devitoproject.org', author="Imperial College London", author_email='opesci@imperial.ac.uk', license='MIT', packages=find_packages(exclude=['docs', 'tests', 'examples']), install_requires=reqs, + extras_require={'extras': opt_reqs}, dependency_links=links, test_suite='tests') diff --git a/tests/conftest.py b/tests/conftest.py index 85eeb24be9..8066bf5839 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -18,16 +18,20 @@ from devito.tools import as_tuple try: - import mpi4py - no_mpi = False + from mpi4py import MPI # noqa except ImportError: - no_mpi = True pass -skipif_yask = pytest.mark.skipif(configuration['backend'] == 'yask', - reason="YASK testing is currently restricted") -skipif_mpi = pytest.mark.skipif(no_mpi, reason="mpi not installed") +def skipif_backend(backends): + conditions = [] + for b in backends: + conditions.append(b == configuration['backend']) + return pytest.mark.skipif(any(conditions), + reason="{} testing is currently restricted".format(b)) + + +skipif_nompi = pytest.mark.skipif(MPI is None, reason="mpi not installed") # Testing dimensions for space and time @@ -57,7 +61,6 @@ def timefunction(name, space_order=1): return TimeFunction(name=name, grid=grid, space_order=space_order) -@pytest.fixture(scope="session") def unit_box(name='a', shape=(11, 11), grid=None): """Create a field with value 0. to 1. in each dimension""" grid = grid or Grid(shape=shape) @@ -67,7 +70,6 @@ def unit_box(name='a', shape=(11, 11), grid=None): return a -@pytest.fixture(scope="session") def unit_box_time(name='a', shape=(11, 11)): """Create a field with value 0. to 1. in each dimension""" grid = Grid(shape=shape) @@ -78,7 +80,6 @@ def unit_box_time(name='a', shape=(11, 11)): return a -@pytest.fixture(scope="session") def points(grid, ranges, npoints, name='points'): """Create a set of sparse points from a set of coordinate ranges for each spatial dimension. @@ -89,7 +90,6 @@ def points(grid, ranges, npoints, name='points'): return points -@pytest.fixture(scope="session") def time_points(grid, ranges, npoints, name='points', nt=10): """Create a set of sparse points from a set of coordinate ranges for each spatial dimension. @@ -287,19 +287,18 @@ def wrapper(*args, **kwargs): return dec -# Support to run MPI tests -# This is partly extracted from: -# `https://github.com/firedrakeproject/firedrake/blob/master/tests/conftest.py` - -mpi_exec = 'mpiexec' -mpi_distro = sniff_mpi_distro(mpi_exec) - - def parallel(item): """Run a test in parallel. :parameter item: The test item to run. """ + # Support to run MPI tests + # This is partly extracted from: + # `https://github.com/firedrakeproject/firedrake/blob/master/tests/conftest.py` + + mpi_exec = 'mpiexec' + mpi_distro = sniff_mpi_distro(mpi_exec) + marker = item.get_closest_marker("parallel") nprocs = as_tuple(marker.kwargs.get("nprocs", 2)) for i in nprocs: diff --git a/tests/test_adjoint.py b/tests/test_adjoint.py index 5e4aaff7ab..b0b8155de5 100644 --- a/tests/test_adjoint.py +++ b/tests/test_adjoint.py @@ -1,13 +1,15 @@ import numpy as np import pytest from numpy import linalg -from conftest import skipif_yask, unit_box, points - -from devito import clear_cache, Operator +from conftest import unit_box, points +from devito import clear_cache, Operator, configuration from devito.logger import info from examples.seismic import demo_model, TimeAxis, RickerSource, Receiver from examples.seismic.acoustic import AcousticWaveSolver +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") presets = { 'constant': {'preset': 'constant-isotropic'}, @@ -15,7 +17,6 @@ } -@skipif_yask class TestAdjoint(object): def setup_method(self, method): diff --git a/tests/test_autotuner.py b/tests/test_autotuner.py index fc99e5a5ac..99c0ba0c4d 100644 --- a/tests/test_autotuner.py +++ b/tests/test_autotuner.py @@ -1,5 +1,4 @@ from __future__ import absolute_import - from functools import reduce from operator import mul try: @@ -9,16 +8,16 @@ from io import StringIO import pytest -from conftest import skipif_yask - import numpy as np - from devito import Grid, Function, TimeFunction, Eq, Operator, configuration, silencio from devito.logger import logger, logging +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + @silencio(log_level='DEBUG') -@skipif_yask @pytest.mark.parametrize("shape,expected", [ ((30, 30), 17), ((30, 30, 30), 21) @@ -61,7 +60,6 @@ def test_at_is_actually_working(shape, expected): @silencio(log_level='DEBUG') -@skipif_yask def test_timesteps_per_at_run(): """ Check that each autotuning run (ie with a given block shape) takes diff --git a/tests/test_checkpointing.py b/tests/test_checkpointing.py index 1919231145..ffc0c4295a 100644 --- a/tests/test_checkpointing.py +++ b/tests/test_checkpointing.py @@ -3,15 +3,17 @@ from examples.seismic import Receiver from pyrevolve import Revolver import numpy as np -from conftest import skipif_yask import pytest from functools import reduce -from devito import Grid, TimeFunction, Operator, Function, Eq, silencio +from devito import Grid, TimeFunction, Operator, Function, Eq, silencio, configuration + +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") @silencio(log_level='WARNING') -@skipif_yask def test_segmented_incremment(): """ Test for segmented operator execution of a one-sided first order @@ -41,7 +43,6 @@ def test_segmented_incremment(): @silencio(log_level='WARNING') -@skipif_yask def test_segmented_fibonacci(): """ Test for segmented operator execution of a one-sided second order @@ -78,7 +79,6 @@ def test_segmented_fibonacci(): @silencio(log_level='WARNING') -@skipif_yask def test_segmented_averaging(): """ Test for segmented operator execution of a two-sided, second order @@ -95,24 +95,23 @@ def test_segmented_averaging(): # We add the average to the point itself, so the grid "interior" # (domain) is updated only. f_ref = TimeFunction(name='f', grid=grid) - f_ref.data_allocated[:] = 1. + f_ref.data_with_halo[:] = 1. op(f=f_ref, time=1) assert (f_ref.data[1, :] == 2.).all() - assert (f_ref.data_allocated[1, 0] == 1.).all() - assert (f_ref.data_allocated[1, -1] == 1.).all() + assert (f_ref.data_with_halo[1, 0] == 1.).all() + assert (f_ref.data_with_halo[1, -1] == 1.).all() # Now we sweep the x direction in 4 segmented steps of 5 iterations each nsteps = 5 - f.data_allocated[:] = 1. + f.data_with_halo[:] = 1. for i in range(4): op(f=f, time=1, x_m=i*nsteps, x_M=(i+1)*nsteps-1) assert (f_ref.data[1, :] == 2.).all() - assert (f_ref.data_allocated[1, 0] == 1.).all() - assert (f_ref.data_allocated[1, -1] == 1.).all() + assert (f_ref.data_with_halo[1, 0] == 1.).all() + assert (f_ref.data_with_halo[1, -1] == 1.).all() @silencio(log_level='WARNING') -@skipif_yask @pytest.mark.parametrize('space_order', [4]) @pytest.mark.parametrize('kernel', ['OT2']) @pytest.mark.parametrize('shape', [(70, 80), (50, 50, 50)]) @@ -151,7 +150,6 @@ def test_forward_with_breaks(shape, kernel, space_order): @silencio(log_level='WARNING') -@skipif_yask def test_acoustic_save_and_nosave(shape=(50, 50), spacing=(15.0, 15.0), tn=500., time_order=2, space_order=4, nbpml=10): """ Run the acoustic example with and without save=True. Make sure the result is the @@ -169,7 +167,6 @@ def test_acoustic_save_and_nosave(shape=(50, 50), spacing=(15.0, 15.0), tn=500., assert(np.allclose(rec.data, rec_bk)) -@skipif_yask def test_index_alignment(const): """ A much simpler test meant to ensure that the forward and reverse indices are correctly aligned (i.e. u * v , where u is the forward field and v the reverse field diff --git a/tests/test_data.py b/tests/test_data.py index 2e8881415b..3965ddc849 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -1,11 +1,15 @@ -from conftest import skipif_yask +from conftest import skipif_backend, skipif_nompi import pytest import numpy as np -from devito import Grid, Function, TimeFunction, Eq, Operator, ALLOC_GUARD, ALLOC_FLAT +from devito import (Grid, Function, TimeFunction, Dimension, Eq, Operator, # noqa + configuration, ALLOC_GUARD, ALLOC_FLAT) from devito.data import Decomposition from devito.types import LEFT, RIGHT +pytestmark = pytest.mark.skipif(configuration['backend'] == 'ops', + reason="testing is currently restricted") + class TestDataBasic(object): @@ -76,7 +80,8 @@ def test_halo_indexing(self): u = Function(name='yu3D', grid=grid, space_order=2) assert u.shape == u.data.shape == domain_shape - assert u.shape_with_halo == u.data_with_halo.shape == (20, 20, 20) + assert u._shape_with_inhalo == u.data_with_halo.shape == (20, 20, 20) + assert u.shape_with_halo == u._shape_with_inhalo # W/o MPI, these two coincide # Test simple insertion and extraction u.data_with_halo[0, 0, 0] = 1. @@ -90,7 +95,7 @@ def test_halo_indexing(self): assert u.data[-1, -1, -1] == 0. assert u.data_with_halo[-1, -1, -1] == 3. - def test_data_arithmetic(self): + def test_arithmetic(self): """ Tests arithmetic operations between :class:`Data` objects and values. """ @@ -119,7 +124,7 @@ def test_data_arithmetic(self): arr.fill(2.) assert np.all(arr - u.data == 1.) - @skipif_yask # YASK not throwing excpetions yet + @skipif_backend(['yask']) # YASK and OPS backends do not support MPI yet def test_illegal_indexing(self): """ Tests that indexing into illegal entries throws an exception. @@ -167,7 +172,8 @@ def test_domain_vs_halo(self): u0 = Function(name='u0', grid=grid, space_order=0) u2 = Function(name='u2', grid=grid, space_order=2) - assert u0.shape == u0.shape_with_halo == u0.shape_allocated + assert u0.shape == u0._shape_with_inhalo == u0.shape_allocated + assert u0.shape_with_halo == u0._shape_with_inhalo # W/o MPI, these two coincide assert len(u2.shape) == len(u2._extent_halo.left) assert tuple(i + j*2 for i, j in zip(u2.shape, u2._extent_halo.left)) ==\ u2.shape_with_halo @@ -193,7 +199,7 @@ def test_domain_vs_halo(self): assert v._extent_padding.left == v._extent_padding.right == (1, 3, 4) -@skipif_yask +@skipif_backend(['yask']) class TestDecomposition(object): """ @@ -309,7 +315,8 @@ def test_reshape_iterable(self): assert d.reshape((1, 3, 10, 11, 14)) == Decomposition([[0], [1], [], [2, 3]], 2) -@skipif_yask # YASK backend does not support MPI yet +@skipif_nompi +@skipif_backend(['yask']) # YASK and OPS backends do not support MPI yet class TestDataDistributed(object): """ @@ -317,7 +324,7 @@ class TestDataDistributed(object): """ @pytest.mark.parallel(nprocs=4) - def test_data_localviews(self): + def test_localviews(self): grid = Grid(shape=(4, 4)) x, y = grid.dimensions glb_pos_map = grid.distributor.glb_pos_map @@ -353,14 +360,12 @@ def test_trivial_insertion(self): u.data[:] = 1. assert np.all(u.data == 1.) assert np.all(u.data._local == 1.) - assert np.all(u.data._global == 1.) v.data_with_halo[:] = 1. assert v.data_with_halo[:].shape == (3, 3) assert np.all(v.data_with_halo == 1.) assert np.all(v.data_with_halo[:] == 1.) assert np.all(v.data_with_halo._local == 1.) - assert np.all(v.data_with_halo._global == 1.) @pytest.mark.parallel(nprocs=4) def test_indexing(self): @@ -407,7 +412,6 @@ def test_slicing(self): # to the local rank domain, so it must be == myrank assert np.all(u.data == myrank) assert np.all(u.data._local == myrank) - assert np.all(u.data._global == myrank) if LEFT in glb_pos_map[x] and LEFT in glb_pos_map[y]: assert np.all(u.data[:2, :2] == myrank) assert u.data[:2, 2:].size == u.data[2:, :2].size == u.data[2:, 2:].size == 0 @@ -534,6 +538,86 @@ def test_from_replicated_to_distributed(self): except: assert False + @pytest.mark.parallel(nprocs=4) + def test_misc_setup(self): + """Test setup of Functions with mixed distributed/replicated Dimensions.""" + grid = Grid(shape=(4, 4)) + _, y = grid.dimensions + dy = Dimension(name='dy') + + # Note: `grid` must be passed to `c` since `x` is a distributed dimension, + # and `grid` carries the `x` decomposition + c = Function(name='c', grid=grid, dimensions=(y, dy), shape=(4, 5)) + + # The following should be identical to `c` in everything but the name + c2 = Function(name='c2', grid=grid, dimensions=(y, dy), shape=(None, 5)) + assert c.shape == c2.shape == (2, 5) + assert c.shape_with_halo == c2.shape_with_halo + assert c._decomposition == c2._decomposition + + # The following should all raise an exception as illegal + try: + Function(name='c3', grid=grid, dimensions=(y, dy)) + assert False + except TypeError: + # Missing `shape` + assert True + + # The following should all raise an exception as illegal + try: + Function(name='c4', grid=grid, dimensions=(y, dy), shape=(3, 5)) + assert False + except ValueError: + # The provided y-size, 3, doesn't match the y-size in grid (4) + assert True + + # The following should all raise an exception as illegal + try: + Function(name='c4', grid=grid, dimensions=(y, dy), shape=(4,)) + assert False + except ValueError: + # Too few entries for `shape` (two expected, for `y` and `dy`) + assert True + + @pytest.mark.parallel(nprocs=4) + def test_misc_data(self): + """ + Test data insertion/indexing for Functions with mixed + distributed/replicated Dimensions. + """ + dx = Dimension(name='dx') + grid = Grid(shape=(4, 4)) + x, y = grid.dimensions + glb_pos_map = grid.distributor.glb_pos_map + + # Note: `grid` must be passed to `c` since `x` is a distributed dimension, + # and `grid` carries the `x` decomposition + c = Function(name='c', grid=grid, dimensions=(x, dx), shape=(4, 5)) + + # Data insertion + for i in range(4): + c.data[i, 0] = 1.0+i + c.data[i, 1] = 1.0+i + c.data[i, 2] = 3.0+i + c.data[i, 3] = 6.0+i + c.data[i, 4] = 5.0+i + + # Data indexing + if LEFT in glb_pos_map[x]: + assert(np.all(c.data[0] == [1., 1., 3., 6., 5.])) + assert(np.all(c.data[1] == [2., 2., 4., 7., 6.])) + else: + assert(np.all(c.data[2] == [3., 3., 5., 8., 7.])) + assert(np.all(c.data[3] == [4., 4., 6., 9., 8.])) + + # Same as before, but with negative indices and non-trivial slices + if LEFT in glb_pos_map[x]: + assert(np.all(c.data[0:-3] == [1., 1., 3., 6., 5.])) + assert(np.all(c.data[-3:-2] == [2., 2., 4., 7., 6.])) + else: + assert(np.all(c.data[-2:-1] == [3., 3., 5., 8., 7.])) + assert(np.all(c.data[-1] == [4., 4., 6., 9., 8.])) + def test_scalar_arg_substitution(t0, t1): """ @@ -573,4 +657,4 @@ def test_oob_guard(): if __name__ == "__main__": from devito import configuration configuration['mpi'] = True - TestDataDistributed().test_indexing_in_views() + TestDataDistributed().test_misc_data() diff --git a/tests/test_dependency_bugs.py b/tests/test_dependency_bugs.py index 2fdd2f01f9..e7141a1fdc 100644 --- a/tests/test_dependency_bugs.py +++ b/tests/test_dependency_bugs.py @@ -1,6 +1,12 @@ import numpy as np +import pytest from numpy.random import rand +from devito import configuration + +pytestmark = pytest.mark.skipif(configuration['backend'] == 'ops', + reason="testing is currently restricted") + def test_numpy_dot(): # Checking for bug in numpy.dot diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index ed0bd7de78..34ecaaf734 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -1,267 +1,250 @@ +from conftest import skipif_backend + import numpy as np import pytest -from conftest import skipif_yask from sympy import Derivative, simplify, diff -from devito import (Grid, Function, TimeFunction, Eq, Operator, - clear_cache, ConditionalDimension, left, right, centered, - staggered_diff) +from devito import (Grid, Function, TimeFunction, Eq, Operator, clear_cache, + ConditionalDimension, left, right, centered, staggered_diff) from devito.finite_differences import Differentiable _PRECISION = 9 -@pytest.fixture -def shape(xdim=20, ydim=30, zdim=20): - return (xdim, ydim, zdim) - - -@pytest.fixture -def grid(shape): - return Grid(shape=shape) - - -@pytest.fixture def x(grid): return grid.dimensions[0] -@pytest.fixture def y(grid): return grid.dimensions[1] -@pytest.fixture def z(grid): return grid.dimensions[2] -@pytest.fixture def t(grid): return grid.stepping_dim -@skipif_yask -@pytest.mark.parametrize('SymbolType, dim', [ - (Function, x), (Function, y), - (TimeFunction, x), (TimeFunction, y), (TimeFunction, t), -]) -def test_stencil_derivative(grid, shape, SymbolType, dim): - """Test symbolic behaviour when expanding stencil derivatives""" - i = dim(grid) # issue fixtures+parametrize: github.com/pytest-dev/pytest/issues/349 - u = SymbolType(name='u', grid=grid) - u.data[:] = 66.6 - di = u.diff(i) - dii = u.diff(i, i) - # Check for sympy Derivative objects - assert(isinstance(di, Derivative) and isinstance(dii, Derivative)) - s_di = di.as_finite_difference([i - i.spacing, i]) - s_dii = dii.as_finite_difference([i - i.spacing, i, i + i.spacing]) - # Check stencil length of first and second derivatives - assert(len(s_di.args) == 2 and len(s_dii.args) == 3) - u_di = s_di.args[0].args[1] - u_dii = s_di.args[0].args[1] - # Ensure that devito meta-data survived symbolic transformation - assert(u_di.grid.shape == shape and u_dii.grid.shape == shape) - assert(u_di.shape == u.shape and u_dii.shape == u.shape) - assert(np.allclose(u_di.data, 66.6)) - assert(np.allclose(u_dii.data, 66.6)) - - -@skipif_yask -@pytest.mark.parametrize('SymbolType, derivative, dim', [ - (Function, 'dx2', 3), (Function, 'dy2', 3), - (TimeFunction, 'dx2', 3), (TimeFunction, 'dy2', 3), (TimeFunction, 'dt', 2) -]) -def test_preformed_derivatives(grid, SymbolType, derivative, dim): - """Test the stencil expressions provided by devito objects""" - u = SymbolType(name='u', grid=grid, time_order=2, space_order=2) - expr = getattr(u, derivative) - assert(len(expr.args) == dim) - - -@skipif_yask -@pytest.mark.parametrize('derivative, dim', [ - ('dx', x), ('dy', y), ('dz', z) -]) -@pytest.mark.parametrize('order', [1, 2, 4, 6, 8, 10, 12, 14, 16]) -def test_derivatives_space(grid, derivative, dim, order): - """Test first derivative expressions against native sympy""" - dim = dim(grid) # issue fixtures+parametrize: github.com/pytest-dev/pytest/issues/349 - u = TimeFunction(name='u', grid=grid, time_order=2, space_order=order) - expr = getattr(u, derivative) - # Establish native sympy derivative expression - width = int(order / 2) - if order == 1: - indices = [dim, dim + dim.spacing] - else: - indices = [(dim + i * dim.spacing) for i in range(-width, width + 1)] - s_expr = u.diff(dim).as_finite_difference(indices).evalf(_PRECISION) - assert(simplify(expr - s_expr) == 0) # Symbolic equality - assert(expr == s_expr) # Exact equailty - - -@skipif_yask -@pytest.mark.parametrize('derivative, dim', [ - ('dx2', x), ('dy2', y), ('dz2', z) -]) -@pytest.mark.parametrize('order', [2, 4, 6, 8, 10, 12, 14, 16]) -def test_second_derivatives_space(grid, derivative, dim, order): - """Test second derivative expressions against native sympy""" - dim = dim(grid) # issue fixtures+parametrize: github.com/pytest-dev/pytest/issues/349 - u = TimeFunction(name='u', grid=grid, time_order=2, space_order=order) - expr = getattr(u, derivative) - # Establish native sympy derivative expression - width = int(order / 2) - indices = [(dim + i * dim.spacing) for i in range(-width, width + 1)] - s_expr = u.diff(dim, dim).as_finite_difference(indices).evalf(_PRECISION) - assert(simplify(expr - s_expr) == 0) # Symbolic equality - assert(expr == s_expr) # Exact equailty - - -@skipif_yask -@pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]) -# Only test x and t as y and z are the same as x -@pytest.mark.parametrize('derivative', ['dx', 'dxl', 'dxr', 'dx2']) -def test_fd_space(derivative, space_order): +@skipif_backend(['yask', 'ops']) +class TestFD(object): """ - This test compares the discrete finite-difference scheme against polynomials - For a given order p, the finite difference scheme should - be exact for polynomials of order p - :param derivative: name of the derivative to be tested - :param space_order: space order of the finite difference stencil + Class for finite difference testing + Tests the accuracy w.r.t polynomials + Test that the shortcut produce the same answer as the FD functions """ - clear_cache() - # dummy axis dimension - nx = 100 - xx = np.linspace(-1, 1, nx) - dx = xx[1] - xx[0] - # Symbolic data - grid = Grid(shape=(nx,), dtype=np.float32) - x = grid.dimensions[0] - u = Function(name="u", grid=grid, space_order=space_order) - du = Function(name="du", grid=grid, space_order=space_order) - # Define polynomial with exact fd - coeffs = np.ones((space_order,), dtype=np.float32) - polynome = sum([coeffs[i]*x**i for i in range(0, space_order)]) - polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32) - # Fill original data with the polynomial values - u.data[:] = polyvalues - # True derivative of the polynome - Dpolynome = diff(diff(polynome)) if derivative == 'dx2' else diff(polynome) - Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx], np.float32) - # FD derivative, symbolic - u_deriv = getattr(u, derivative) - # Compute numerical FD - stencil = Eq(du, u_deriv) - op = Operator(stencil, subs={x.spacing: dx}) - op.apply() - - # Check exactness of the numerical derivative except inside space_brd - space_border = space_order - error = abs(du.data[space_border:-space_border] - - Dpolyvalues[space_border:-space_border]) - assert np.isclose(np.mean(error), 0., atol=1e-3) - -@skipif_yask -@pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]) -@pytest.mark.parametrize('stagger', [centered, right, left]) -# Only test x and t as y and z are the same as x -def test_fd_space_staggered(space_order, stagger): - """ - This test compares the discrete finite-difference scheme against polynomials - For a given order p, the finite difference scheme should - be exact for polynomials of order p - :param derivative: name of the derivative to be tested - :param space_order: space order of the finite difference stencil - """ - clear_cache() - if stagger == left: - off = -.5 - elif stagger == right: - off = .5 - else: - off = 0 - # dummy axis dimension - nx = 100 - xx = np.linspace(-1, 1, nx) - dx = xx[1] - xx[0] - # Location of the staggered function - xx2 = xx + off * dx - # Symbolic data - grid = Grid(shape=(nx,), dtype=np.float32) - x = grid.dimensions[0] - u = Function(name="u", grid=grid, space_order=space_order, stagger=(1,)) - du = Function(name="du", grid=grid, space_order=space_order) - # Define polynomial with exact fd - coeffs = np.ones((space_order,), dtype=np.float32) - polynome = sum([coeffs[i]*x**i for i in range(0, space_order)]) - polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32) - # Fill original data with the polynomial values - u.data[:] = polyvalues - # True derivative of the polynome - Dpolynome = diff(polynome) - Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx2], np.float32) - # FD derivative, symbolic - u_deriv = staggered_diff(u, deriv_order=1, fd_order=space_order, - dim=x, stagger=stagger) - # Compute numerical FD - stencil = Eq(du, u_deriv) - op = Operator(stencil, subs={x.spacing: dx}) - op.apply() - - # Check exactness of the numerical derivative except inside space_brd - space_border = space_order - error = abs(du.data[space_border:-space_border] - - Dpolyvalues[space_border:-space_border]) - - assert np.isclose(np.mean(error), 0., atol=1e-3) - - -@skipif_yask -def test_subsampled_fd(): - """ - Test that the symbolic interface is working for space subsampled - functions. - """ - nt = 19 - grid = Grid(shape=(12, 12), extent=(11, 11)) - - u = TimeFunction(name='u', grid=grid, save=nt, space_order=2) - assert(grid.time_dim in u.indices) - - # Creates subsampled spatial dimensions and according grid - dims = tuple([ConditionalDimension(d.name+'sub', parent=d, factor=2) - for d in u.grid.dimensions]) - grid2 = Grid((6, 6), dimensions=dims) - u2 = TimeFunction(name='u2', grid=grid2, save=nt, space_order=1) - for i in range(nt): - for j in range(u2.data_allocated.shape[2]): - u2.data_allocated[i, :, j] = np.arange(u2.data_allocated.shape[2]) - - eqns = [Eq(u.forward, u + 1.), Eq(u2.forward, u2.dx)] - op = Operator(eqns, dse="advanced") - op.apply(time_M=nt-2) - # Verify that u2[1, x,y]= du2/dx[0, x, y] - - assert np.allclose(u.data[-1], nt-1) - assert np.allclose(u2.data[1], 0.5) - - -@skipif_yask -@pytest.mark.parametrize('expr,expected', [ - ('f.dx', '-f(x)/h_x + f(x + h_x)/h_x'), - ('f.dx + g.dx', '-f(x)/h_x + f(x + h_x)/h_x - g(x)/h_x + g(x + h_x)/h_x'), - ('-f', '-f(x)'), - ('-(f + g)', '-f(x) - g(x)') -]) -def test_shortcuts(expr, expected): - grid = Grid(shape=(1,)) - f = Function(name='f', grid=grid) # noqa - g = Function(name='g', grid=grid) # noqa - - expr = eval(expr) - - assert isinstance(expr, Differentiable) - assert expected == str(expr) + def setup_method(self): + self.shape = (20, 20, 20) + self.grid = Grid(self.shape) + + @pytest.mark.parametrize('SymbolType, dim', [ + (Function, x), (Function, y), + (TimeFunction, x), (TimeFunction, y), (TimeFunction, t), + ]) + def test_stencil_derivative(self, SymbolType, dim): + """Test symbolic behaviour when expanding stencil derivatives""" + i = dim(self.grid) + u = SymbolType(name='u', grid=self.grid) + u.data[:] = 66.6 + di = u.diff(i) + dii = u.diff(i, i) + # Check for sympy Derivative objects + assert(isinstance(di, Derivative) and isinstance(dii, Derivative)) + s_di = di.as_finite_difference([i - i.spacing, i]) + s_dii = dii.as_finite_difference([i - i.spacing, i, i + i.spacing]) + # Check stencil length of first and second derivatives + assert(len(s_di.args) == 2 and len(s_dii.args) == 3) + u_di = s_di.args[0].args[1] + u_dii = s_di.args[0].args[1] + # Ensure that devito meta-data survived symbolic transformation + assert(u_di.grid.shape == self.shape and u_dii.grid.shape == self.shape) + assert(u_di.shape == u.shape and u_dii.shape == u.shape) + assert(np.allclose(u_di.data, 66.6)) + assert(np.allclose(u_dii.data, 66.6)) + + @pytest.mark.parametrize('SymbolType, derivative, dim', [ + (Function, 'dx2', 3), (Function, 'dy2', 3), + (TimeFunction, 'dx2', 3), (TimeFunction, 'dy2', 3), (TimeFunction, 'dt', 2) + ]) + def test_preformed_derivatives(self, SymbolType, derivative, dim): + """Test the stencil expressions provided by devito objects""" + u = SymbolType(name='u', grid=self.grid, time_order=2, space_order=2) + expr = getattr(u, derivative) + assert(len(expr.args) == dim) + + @pytest.mark.parametrize('derivative, dim', [ + ('dx', x), ('dy', y), ('dz', z) + ]) + @pytest.mark.parametrize('order', [1, 2, 4, 6, 8, 10, 12, 14, 16]) + def test_derivatives_space(self, derivative, dim, order): + """Test first derivative expressions against native sympy""" + dim = dim(self.grid) + u = TimeFunction(name='u', grid=self.grid, time_order=2, space_order=order) + expr = getattr(u, derivative) + # Establish native sympy derivative expression + width = int(order / 2) + if order == 1: + indices = [dim, dim + dim.spacing] + else: + indices = [(dim + i * dim.spacing) for i in range(-width, width + 1)] + s_expr = u.diff(dim).as_finite_difference(indices).evalf(_PRECISION) + assert(simplify(expr - s_expr) == 0) # Symbolic equality + assert(expr == s_expr) # Exact equailty + + @pytest.mark.parametrize('derivative, dim', [ + ('dx2', x), ('dy2', y), ('dz2', z) + ]) + @pytest.mark.parametrize('order', [2, 4, 6, 8, 10, 12, 14, 16]) + def test_second_derivatives_space(self, derivative, dim, order): + """Test second derivative expressions against native sympy""" + dim = dim(self.grid) + u = TimeFunction(name='u', grid=self.grid, time_order=2, space_order=order) + expr = getattr(u, derivative) + # Establish native sympy derivative expression + width = int(order / 2) + indices = [(dim + i * dim.spacing) for i in range(-width, width + 1)] + s_expr = u.diff(dim, dim).as_finite_difference(indices).evalf(_PRECISION) + assert(simplify(expr - s_expr) == 0) # Symbolic equality + assert(expr == s_expr) # Exact equailty + + @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]) + # Only test x and t as y and z are the same as x + @pytest.mark.parametrize('derivative', ['dx', 'dxl', 'dxr', 'dx2']) + def test_fd_space(self, derivative, space_order): + """ + This test compares the discrete finite-difference scheme against polynomials + For a given order p, the finite difference scheme should + be exact for polynomials of order p + :param derivative: name of the derivative to be tested + :param space_order: space order of the finite difference stencil + """ + clear_cache() + # dummy axis dimension + nx = 100 + xx = np.linspace(-1, 1, nx) + dx = xx[1] - xx[0] + # Symbolic data + grid = Grid(shape=(nx,), dtype=np.float32) + x = grid.dimensions[0] + u = Function(name="u", grid=grid, space_order=space_order) + du = Function(name="du", grid=grid, space_order=space_order) + # Define polynomial with exact fd + coeffs = np.ones((space_order,), dtype=np.float32) + polynome = sum([coeffs[i]*x**i for i in range(0, space_order)]) + polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32) + # Fill original data with the polynomial values + u.data[:] = polyvalues + # True derivative of the polynome + Dpolynome = diff(diff(polynome)) if derivative == 'dx2' else diff(polynome) + Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx], np.float32) + # FD derivative, symbolic + u_deriv = getattr(u, derivative) + # Compute numerical FD + stencil = Eq(du, u_deriv) + op = Operator(stencil, subs={x.spacing: dx}) + op.apply() + + # Check exactness of the numerical derivative except inside space_brd + space_border = space_order + error = abs(du.data[space_border:-space_border] - + Dpolyvalues[space_border:-space_border]) + assert np.isclose(np.mean(error), 0., atol=1e-3) + + @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]) + @pytest.mark.parametrize('stagger', [centered, right, left]) + # Only test x and t as y and z are the same as x + def test_fd_space_staggered(self, space_order, stagger): + """ + This test compares the discrete finite-difference scheme against polynomials + For a given order p, the finite difference scheme should + be exact for polynomials of order p + :param derivative: name of the derivative to be tested + :param space_order: space order of the finite difference stencil + """ + clear_cache() + if stagger == left: + off = -.5 + elif stagger == right: + off = .5 + else: + off = 0 + # dummy axis dimension + nx = 100 + xx = np.linspace(-1, 1, nx) + dx = xx[1] - xx[0] + # Location of the staggered function + xx2 = xx + off * dx + # Symbolic data + grid = Grid(shape=(nx,), dtype=np.float32) + x = grid.dimensions[0] + u = Function(name="u", grid=grid, space_order=space_order, stagger=(1,)) + du = Function(name="du", grid=grid, space_order=space_order) + # Define polynomial with exact fd + coeffs = np.ones((space_order,), dtype=np.float32) + polynome = sum([coeffs[i]*x**i for i in range(0, space_order)]) + polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32) + # Fill original data with the polynomial values + u.data[:] = polyvalues + # True derivative of the polynome + Dpolynome = diff(polynome) + Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx2], np.float32) + # FD derivative, symbolic + u_deriv = staggered_diff(u, deriv_order=1, fd_order=space_order, + dim=x, stagger=stagger) + # Compute numerical FD + stencil = Eq(du, u_deriv) + op = Operator(stencil, subs={x.spacing: dx}) + op.apply() + + # Check exactness of the numerical derivative except inside space_brd + space_border = space_order + error = abs(du.data[space_border:-space_border] - + Dpolyvalues[space_border:-space_border]) + + assert np.isclose(np.mean(error), 0., atol=1e-3) + + def test_subsampled_fd(self): + """ + Test that the symbolic interface is working for space subsampled + functions. + """ + nt = 19 + grid = Grid(shape=(12, 12), extent=(11, 11)) + + u = TimeFunction(name='u', grid=grid, save=nt, space_order=2) + assert(grid.time_dim in u.indices) + + # Creates subsampled spatial dimensions and according grid + dims = tuple([ConditionalDimension(d.name+'sub', parent=d, factor=2) + for d in u.grid.dimensions]) + grid2 = Grid((6, 6), dimensions=dims) + u2 = TimeFunction(name='u2', grid=grid2, save=nt, space_order=1) + for i in range(nt): + for j in range(u2.data_with_halo.shape[2]): + u2.data_with_halo[i, :, j] = np.arange(u2.data_with_halo.shape[2]) + + eqns = [Eq(u.forward, u + 1.), Eq(u2.forward, u2.dx)] + op = Operator(eqns, dse="advanced") + op.apply(time_M=nt-2) + # Verify that u2[1, x,y]= du2/dx[0, x, y] + + assert np.allclose(u.data[-1], nt-1) + assert np.allclose(u2.data[1], 0.5) + + @pytest.mark.parametrize('expr,expected', [ + ('f.dx', '-f(x)/h_x + f(x + h_x)/h_x'), + ('f.dx + g.dx', '-f(x)/h_x + f(x + h_x)/h_x - g(x)/h_x + g(x + h_x)/h_x'), + ('-f', '-f(x)'), + ('-(f + g)', '-f(x) - g(x)') + ]) + def test_shortcuts(self, expr, expected): + grid = Grid(shape=(10,)) + f = Function(name='f', grid=grid) # noqa + g = Function(name='g', grid=grid) # noqa + + expr = eval(expr) + + assert isinstance(expr, Differentiable) + assert expected == str(expr) diff --git a/tests/test_dimension.py b/tests/test_dimension.py index b0e899c22e..f5058e8698 100644 --- a/tests/test_dimension.py +++ b/tests/test_dimension.py @@ -4,13 +4,14 @@ from sympy import And import pytest -from conftest import skipif_yask, configuration_override +from conftest import skipif_backend, configuration_override from devito import (ConditionalDimension, Grid, Function, TimeFunction, SparseFunction, # noqa Eq, Operator, Constant, SubDimension) from devito.ir.iet import Iteration, FindNodes, retrieve_iteration_tree +@skipif_backend(['ops']) class TestSubDimension(object): def test_interior(self): @@ -35,7 +36,7 @@ def test_interior(self): assert np.all(u.data[1, :, :, 0] == 0.) assert np.all(u.data[1, :, :, -1] == 0.) - @skipif_yask + @skipif_backend(['yask']) def test_domain_vs_interior(self): """ Tests application of an Operator consisting of two equations, one @@ -125,7 +126,7 @@ def test_bcs(self): for i in range(1, thickness + 1)) assert np.all(u.data[0, thickness:-thickness, thickness:-thickness] == 1.) - @skipif_yask + @skipif_backend(['yask']) def test_flow_detection_interior(self): """ Test detection of flow directions when :class:`SubDimension`s are used @@ -175,7 +176,7 @@ def test_flow_detection_interior(self): assert np.all(u.data[1, :, 0:5] == 0) assert np.all(u.data[1, :, 6:] == 0) - @skipif_yask + @skipif_backend(['yask']) @pytest.mark.parametrize('exprs,expected,', [ # Carried dependence in both /t/ and /x/ (['Eq(u[t+1, x, y], u[t+1, x-1, y] + u[t, x, y])'], 'y'), @@ -208,7 +209,7 @@ def test_iteration_property_parallel(self, exprs, expected): assert all(i.is_Sequential for i in iterations if i.dim.name != expected) assert all(i.is_Parallel for i in iterations if i.dim.name == expected) - @skipif_yask + @skipif_backend(['yask']) @pytest.mark.parametrize('exprs,expected,', [ (['Eq(u[t, x, yleft], u[t, x, yleft] + 1.)'], ['yleft']), # All outers are parallel, carried dependence in `yleft`, so no SIMD in `yleft` @@ -412,7 +413,7 @@ def test_subdim_fd(self): assert np.all(u.data[1, 1:18, 1:18] == 0.) -@skipif_yask +@skipif_backend(['yask', 'ops']) class TestConditionalDimension(object): """A collection of tests to check the correct functioning of diff --git a/tests/test_dle.py b/tests/test_dle.py index 9b7566663e..8382f08ab1 100644 --- a/tests/test_dle.py +++ b/tests/test_dle.py @@ -1,18 +1,20 @@ from __future__ import absolute_import - from functools import reduce from operator import mul import numpy as np import pytest -from conftest import skipif_yask, EVAL - +from conftest import EVAL from devito.dle import transform -from devito import Grid, Function, TimeFunction, Eq, Operator, solve +from devito import Grid, Function, TimeFunction, Eq, Operator, solve, configuration from devito.ir.equations import DummyEq from devito.ir.iet import (ELEMENTAL, Expression, Callable, Iteration, List, tagger, Transformer, FindNodes, iet_analyze, retrieve_iteration_tree) from unittest.mock import patch +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + @pytest.fixture(scope="module") def exprs(a, b, c, d, a_dense, b_dense): @@ -125,7 +127,6 @@ def _new_operator3(shape, **kwargs): return u.data[1, :], op -@skipif_yask def test_create_elemental_functions_simple(simple_function): roots = [i[-1] for i in retrieve_iteration_tree(simple_function)] retagged = [i._rebuild(properties=tagger(0)) for i in roots] @@ -168,7 +169,6 @@ def test_create_elemental_functions_simple(simple_function): }""") -@skipif_yask def test_create_elemental_functions_complex(complex_function): roots = [i[-1] for i in retrieve_iteration_tree(complex_function)] retagged = [j._rebuild(properties=tagger(i)) for i, j in enumerate(roots)] @@ -233,7 +233,6 @@ def test_create_elemental_functions_complex(complex_function): }""") -@skipif_yask @pytest.mark.parametrize("blockinner,expected", [ (False, 4), (True, 8) @@ -264,7 +263,6 @@ def test_cache_blocking_structure(blockinner, expected): assert 'omp for' in outermost.pragmas[0].value -@skipif_yask @pytest.mark.parametrize("shape", [(10,), (10, 45), (10, 31, 45)]) @pytest.mark.parametrize("blockshape", [2, 7, (3, 3), (2, 9, 1)]) @pytest.mark.parametrize("blockinner", [False, True]) @@ -277,7 +275,6 @@ def test_cache_blocking_no_time_loop(shape, blockshape, blockinner): assert np.equal(wo_blocking.data, w_blocking.data).all() -@skipif_yask @pytest.mark.parametrize("shape", [(20, 33), (45, 31, 45)]) @pytest.mark.parametrize("time_order", [2]) @pytest.mark.parametrize("blockshape", [2, (13, 20), (11, 15, 23)]) @@ -291,7 +288,6 @@ def test_cache_blocking_time_loop(shape, time_order, blockshape, blockinner): assert np.equal(wo_blocking.data, w_blocking.data).all() -@skipif_yask @pytest.mark.parametrize("shape,blockshape", [ ((25, 25, 46), (None, None, None)), ((25, 25, 46), (7, None, None)), @@ -314,7 +310,6 @@ def test_cache_blocking_edge_cases(shape, blockshape): assert np.equal(wo_blocking.data, w_blocking.data).all() -@skipif_yask @pytest.mark.parametrize("shape,blockshape", [ ((3, 3), (3, 4)), ((4, 4), (3, 4)), @@ -338,7 +333,6 @@ def test_cache_blocking_edge_cases_highorder(shape, blockshape): assert np.equal(wo_blocking.data, w_blocking.data).all() -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ # trivial 1D (['Eq(fa[x], fa[x] + fb[x])'], @@ -399,7 +393,6 @@ def test_loops_ompized(fa, fb, fc, fd, t0, t1, t2, t3, exprs, expected, iters): assert 'omp for' not in k.value -@skipif_yask @pytest.mark.parametrize("shape", [(41,), (20, 33), (45, 31, 45)]) def test_composite_transformation(shape): wo_blocking, _ = _new_operator1(shape, dle='noop') @@ -408,7 +401,6 @@ def test_composite_transformation(shape): assert np.equal(wo_blocking.data, w_blocking.data).all() -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ # trivial 1D (['Eq(fe[x,y,z], fe[x,y,z] + fe[x,y,z])'], diff --git a/tests/test_dse.py b/tests/test_dse.py index d8bd200f7e..f86b33d054 100644 --- a/tests/test_dse.py +++ b/tests/test_dse.py @@ -1,12 +1,11 @@ from conftest import EVAL - from sympy import sin # noqa import numpy as np import pytest -from conftest import x, y, z, skipif_yask # noqa +from conftest import x, y, z, skipif_backend # noqa from devito import (Eq, Inc, Constant, Function, TimeFunction, SparseFunction, # noqa - Grid, Operator, ruido) + Grid, Operator, ruido, configuration) from devito.ir import Stencil, FlowGraph, retrieve_iteration_tree from devito.dse import common_subexprs_elimination, collect from devito.symbolics import (xreplace_constrained, iq_timeinvariant, iq_timevarying, @@ -17,6 +16,10 @@ from examples.seismic import demo_model, TimeAxis, RickerSource, GaborSource, Receiver from examples.seismic.tti import AnisotropicWaveSolver +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + # Acoustic @@ -52,7 +55,6 @@ def run_acoustic_forward(dse=None): return u, rec -@skipif_yask def test_acoustic_rewrite_basic(): ret1 = run_acoustic_forward(dse=None) ret2 = run_acoustic_forward(dse='basic') @@ -101,7 +103,6 @@ def tti_nodse(): return v, rec -@skipif_yask def test_tti_rewrite_basic(tti_nodse): operator = tti_operator(dse='basic') rec, u, v, _ = operator.forward() @@ -110,7 +111,6 @@ def test_tti_rewrite_basic(tti_nodse): assert np.allclose(tti_nodse[1].data, rec.data, atol=10e-3) -@skipif_yask def test_tti_rewrite_advanced(tti_nodse): operator = tti_operator(dse='advanced') rec, u, v, _ = operator.forward() @@ -119,7 +119,6 @@ def test_tti_rewrite_advanced(tti_nodse): assert np.allclose(tti_nodse[1].data, rec.data, atol=10e-1) -@skipif_yask def test_tti_rewrite_speculative(tti_nodse): operator = tti_operator(dse='speculative') rec, u, v, _ = operator.forward() @@ -128,7 +127,6 @@ def test_tti_rewrite_speculative(tti_nodse): assert np.allclose(tti_nodse[1].data, rec.data, atol=10e-1) -@skipif_yask def test_tti_rewrite_aggressive(tti_nodse): operator = tti_operator(dse='aggressive') rec, u, v, _ = operator.forward() @@ -138,7 +136,7 @@ def test_tti_rewrite_aggressive(tti_nodse): @ruido(profiling='advanced') -@skipif_yask +@skipif_backend(['yask', 'ops']) @pytest.mark.parametrize('kernel,space_order,expected', [ ('shifted', 8, 355), ('shifted', 16, 622), ('centered', 8, 168), ('centered', 16, 300) @@ -152,7 +150,6 @@ def test_tti_rewrite_aggressive_opcounts(kernel, space_order, expected): # DSE manipulation -@skipif_yask def test_scheduling_after_rewrite(): """Tests loop scheduling after DSE-induced expression hoisting.""" grid = Grid((10, 10)) @@ -175,7 +172,6 @@ def test_scheduling_after_rewrite(): assert trees[1][0].dim == trees[2][0].dim == trees[3][0].dim == grid.time_dim -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ # simple (['Eq(ti1, 4.)', 'Eq(ti0, 3.)', 'Eq(tu, ti0 + ti1 + 5.)'], @@ -201,7 +197,6 @@ def test_xreplace_constrained_time_invariants(tu, tv, tw, ti0, ti1, t0, t1, assert all(str(i.rhs) == j for i, j in zip(found, expected)) -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ # simple (['Eq(ti0, 3.)', 'Eq(tv, 2.4)', 'Eq(tu, tv + 5. + ti0)'], @@ -228,7 +223,6 @@ def test_xreplace_constrained_time_varying(tu, tv, tw, ti0, ti1, t0, t1, assert all(str(i.rhs) == j for i, j in zip(found, expected)) -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ # simple (['Eq(tu, (tv + tw + 5.)*(ti0 + ti1) + (t0 + t1)*(ti0 + ti1))'], @@ -251,7 +245,6 @@ def test_common_subexprs_elimination(tu, tv, tw, ti0, ti1, t0, t1, exprs, expect assert all(str(i.rhs) == j for i, j in zip(processed, expected)) -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ (['Eq(t0, 3.)', 'Eq(t1, 7.)', 'Eq(ti0, t0*3. + 2.)', 'Eq(ti1, t1 + t0 + 1.5)', 'Eq(tv, (ti0 + ti1)*t0)', 'Eq(tw, (ti0 + ti1)*t1)', @@ -266,7 +259,6 @@ def test_graph_trace(tu, tv, tw, ti0, ti1, t0, t1, exprs, expected): assert set([j.lhs for j in g.trace(i)]) == mapper[i] -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ # trivial (['Eq(t0, 1.)', 'Eq(t1, fa[x] + fb[x])'], @@ -292,7 +284,6 @@ def test_graph_isindex(fa, fb, fc, t0, t1, t2, exprs, expected): assert g.is_index(k) == v -@skipif_yask @pytest.mark.parametrize('expr,expected', [ ('2*fa[x] + fb[x]', '2*fa[x] + fb[x]'), ('fa[x]**2', 'fa[x]*fa[x]'), @@ -307,7 +298,6 @@ def test_pow_to_mul(fa, fb, expr, expected): assert str(pow_to_mul(eval(expr))) == expected -@skipif_yask @pytest.mark.parametrize('exprs,expected', [ # none (different distance) (['Eq(t0, fa[x] + fb[x])', 'Eq(t1, fa[x+1] + fb[x])'], @@ -342,7 +332,6 @@ def test_collect_aliases(fa, fb, fc, fd, t0, t1, t2, t3, exprs, expected): assert (len(v.aliased) == 1 and mapper[k] is None) or v.anti_stencil == mapper[k] -@skipif_yask @pytest.mark.parametrize('expr,expected', [ ('Eq(t0, t1)', 0), ('Eq(t0, fa[x] + fb[x])', 1), @@ -361,7 +350,6 @@ def test_estimate_cost(fa, fb, fc, t0, t1, t2, expr, expected): assert estimate_cost(EVAL(expr, fa, fb, fc, t0, t1, t2)) == expected -@skipif_yask @pytest.mark.parametrize('exprs,exp_u,exp_v', [ (['Eq(s, 0)', 'Eq(s, s + 4)', 'Eq(u, s)'], 4, 0), (['Eq(s, 0)', 'Eq(s, s + s + 4)', 'Eq(s, s + 4)', 'Eq(u, s)'], 8, 0), diff --git a/tests/test_gradient.py b/tests/test_gradient.py index 380b8fdc61..eb70d8cd02 100644 --- a/tests/test_gradient.py +++ b/tests/test_gradient.py @@ -1,14 +1,15 @@ import numpy as np import pytest from numpy import linalg -from conftest import skipif_yask - -from devito import Function, info, clear_cache +from devito import Function, info, clear_cache, configuration from examples.seismic.acoustic.acoustic_example import smooth, acoustic_setup as setup from examples.seismic import Receiver +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + -@skipif_yask class TestGradient(object): def setup_method(self, method): @@ -64,7 +65,6 @@ def test_gradient_checkpointing(self, shape, kernel, space_order): gradient2, _ = wave.gradient(residual, u0, m=m0, checkpointing=False) assert np.allclose(gradient.data, gradient2.data) - @skipif_yask @pytest.mark.parametrize('space_order', [4]) @pytest.mark.parametrize('kernel', ['OT2']) @pytest.mark.parametrize('shape', [(70, 80)]) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 9b6d3b453f..a8335d25ce 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -1,17 +1,19 @@ import numpy as np import pytest -from conftest import skipif_yask, unit_box, points, unit_box_time, time_points +from conftest import unit_box, points, unit_box_time, time_points from math import sin, floor - from devito.cgen_utils import FLOAT from devito import (Grid, Operator, Function, SparseFunction, Dimension, TimeFunction, PrecomputedSparseFunction, - PrecomputedSparseTimeFunction) + PrecomputedSparseTimeFunction, configuration) from examples.seismic import demo_model, TimeAxis, RickerSource, Receiver from examples.seismic.acoustic import AcousticWaveSolver +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + -@pytest.fixture def a(shape=(11, 11)): grid = Grid(shape=shape) a = Function(name='a', grid=grid) @@ -21,7 +23,6 @@ def a(shape=(11, 11)): return a -@pytest.fixture def at(shape=(11, 11)): grid = Grid(shape=shape) a = TimeFunction(name='a', grid=grid) @@ -62,7 +63,6 @@ def precompute_linear_interpolation(points, grid, origin): return gridpoints, coefficients -@skipif_yask def test_precomputed_interpolation(): """ Test interpolation with PrecomputedSparseFunction which accepts precomputed values for coefficients @@ -94,7 +94,6 @@ def init(data): assert(all(np.isclose(sf.data, expected_values, rtol=1e-6))) -@skipif_yask def test_precomputed_interpolation_time(): """ Test interpolation with PrecomputedSparseFunction which accepts precomputed values for coefficients, but this time with a TimeFunction @@ -127,7 +126,6 @@ def test_precomputed_interpolation_time(): assert all(np.isclose(sf.data[it, :], it)) -@skipif_yask @pytest.mark.parametrize('shape, coords', [ ((11, 11), [(.05, .9), (.01, .8)]), ((11, 11, 11), [(.05, .9), (.01, .8), (0.07, 0.84)]) @@ -146,7 +144,6 @@ def test_interpolate(shape, coords, npoints=20): assert np.allclose(p.data[:], xcoords, rtol=1e-6) -@skipif_yask @pytest.mark.parametrize('shape, coords', [ ((11, 11), [(.05, .9), (.01, .8)]), ((11, 11, 11), [(.05, .9), (.01, .8), (0.07, 0.84)]) @@ -166,7 +163,6 @@ def test_interpolate_cumm(shape, coords, npoints=20): assert np.allclose(p.data[:], xcoords + 1., rtol=1e-6) -@skipif_yask @pytest.mark.parametrize('shape, coords', [ ((11, 11), [(.05, .9), (.01, .8)]), ((11, 11, 11), [(.05, .9), (.01, .8), (0.07, 0.84)]) @@ -200,7 +196,6 @@ def test_interpolate_time_shift(shape, coords, npoints=20): assert np.allclose(p.data[1, :], xcoords, rtol=1e-6) -@skipif_yask @pytest.mark.parametrize('shape, coords', [ ((11, 11), [(.05, .9), (.01, .8)]), ((11, 11, 11), [(.05, .9), (.01, .8), (0.07, 0.84)]) @@ -219,7 +214,6 @@ def test_interpolate_array(shape, coords, npoints=20): assert np.allclose(p.data[:], xcoords, rtol=1e-6) -@skipif_yask @pytest.mark.parametrize('shape, coords', [ ((11, 11), [(.05, .9), (.01, .8)]), ((11, 11, 11), [(.05, .9), (.01, .8), (0.07, 0.84)]) @@ -241,7 +235,6 @@ def test_interpolate_custom(shape, coords, npoints=20): assert np.allclose(p.data[2, :], 2.0 * xcoords, rtol=1e-6) -@skipif_yask @pytest.mark.parametrize('shape, coords', [ ((11, 11), [(.05, .9), (.01, .8)]), ((11, 11, 11), [(.05, .9), (.01, .8), (0.07, 0.84)]) @@ -264,7 +257,6 @@ def test_interpolate_indexed(shape, coords, npoints=20): assert np.allclose(p.data[2, :], 2.0 * xcoords, rtol=1e-6) -@skipif_yask @pytest.mark.parametrize('shape, coords, result', [ ((11, 11), [(.05, .95), (.45, .45)], 1.), ((11, 11, 11), [(.05, .95), (.45, .45), (.45, .45)], 0.5) @@ -286,7 +278,6 @@ def test_inject(shape, coords, result, npoints=19): assert np.allclose(a.data[indices], result, rtol=1.e-5) -@skipif_yask @pytest.mark.parametrize('shape, coords, result', [ ((11, 11), [(.05, .95), (.45, .45)], 1.), ((11, 11, 11), [(.05, .95), (.45, .45), (.45, .45)], 0.5) @@ -327,7 +318,6 @@ def test_inject_time_shift(shape, coords, result, npoints=19): assert np.allclose(a.data[indices], result, rtol=1.e-5) -@skipif_yask @pytest.mark.parametrize('shape, coords, result', [ ((11, 11), [(.05, .95), (.45, .45)], 1.), ((11, 11, 11), [(.05, .95), (.45, .45), (.45, .45)], 0.5) @@ -350,7 +340,6 @@ def test_inject_array(shape, coords, result, npoints=19): assert np.allclose(a.data[indices], result, rtol=1.e-5) -@skipif_yask @pytest.mark.parametrize('shape, coords, result', [ ((11, 11), [(.05, .95), (.45, .45)], 1.), ((11, 11, 11), [(.05, .95), (.45, .45), (.45, .45)], 0.5) @@ -373,7 +362,6 @@ def test_inject_from_field(shape, coords, result, npoints=19): assert np.allclose(a.data[indices], result, rtol=1.e-5) -@skipif_yask @pytest.mark.parametrize('shape', [(50, 50, 50)]) def test_position(shape): t0 = 0.0 # Start time diff --git a/tests/test_ir.py b/tests/test_ir.py index 572956473c..48aaeae39f 100644 --- a/tests/test_ir.py +++ b/tests/test_ir.py @@ -1,10 +1,11 @@ import pytest -from conftest import EVAL, time, x, y, z, skipif_yask # noqa +from conftest import EVAL, time, x, y, z, skipif_backend # noqa import numpy as np -from devito import Eq, Inc, Grid, Function, TimeFunction, Operator, Dimension # noqa +from devito import (Eq, Inc, Grid, Function, TimeFunction, # noqa + Operator, Dimension, configuration) from devito.ir.equations import DummyEq, LoweredEq from devito.ir.equations.algorithms import dimension_sort from devito.ir.iet.nodes import Conditional, Expression, Iteration @@ -15,8 +16,11 @@ from devito.ir.support.utils import detect_flow_directions from devito.symbolics import indexify +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + -@skipif_yask class TestVectorDistanceArithmetic(object): @pytest.fixture @@ -178,7 +182,6 @@ def test_timed_access_cmp(self, ta_literal): assert False -@skipif_yask class TestSpace(object): def test_intervals_intersection(self): @@ -304,7 +307,6 @@ def test_intervals_subtract(self): assert ix3.subtract(ix) == ix2 -@skipif_yask class TestDependenceAnalysis(object): @pytest.mark.parametrize('expr,expected', [ @@ -462,7 +464,6 @@ def test_flow_detection(self): assert all(mapper.get(i) == {Any} for i in grid.dimensions) -@skipif_yask class TestIET(object): def test_nodes_conditional(self, fc): @@ -630,7 +631,6 @@ def test_loop_wrapping(self, exprs, wrappable): assert all(not i.is_Wrappable for i in iters if i is not time_iter) -@skipif_yask class TestEquationAlgorithms(object): @pytest.mark.parametrize('expr,expected', [ diff --git a/tests/test_mpi.py b/tests/test_mpi.py index eb09285db6..16d1902564 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -1,19 +1,24 @@ import numpy as np - import pytest -from conftest import skipif_yask, skipif_mpi +from conftest import skipif_backend, skipif_nompi # noqa from devito import (Grid, Constant, Function, TimeFunction, SparseFunction, SparseTimeFunction, Dimension, ConditionalDimension, - SubDimension, Eq, Inc, Operator) + SubDimension, Eq, Inc, Operator, norm, inner, configuration) from devito.ir.iet import Call, Conditional, FindNodes from devito.mpi import MPI, copy, sendrecv, update_halo -from devito.parameters import configuration from devito.types import LEFT, RIGHT +from examples.seismic import demo_model, TimeAxis, RickerSource, Receiver +from examples.seismic.acoustic import AcousticWaveSolver + +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + -@skipif_yask -@skipif_mpi +@skipif_nompi +@skipif_backend(['yask', 'ops']) class TestDistributor(object): @pytest.mark.parallel(nprocs=[2, 4]) @@ -63,8 +68,8 @@ def test_ctypes_neighbours(self): assert all(getattr(obj.value._obj, k) == v for k, v in mapper.items()) -@skipif_yask -@skipif_mpi +@skipif_nompi +@skipif_backend(['yask', 'ops']) class TestFunction(object): @pytest.mark.parallel(nprocs=9) @@ -132,13 +137,13 @@ def test_halo_exchange_bilateral(self): glb_pos_map = grid.distributor.glb_pos_map if LEFT in glb_pos_map[y]: - assert np.all(f._data_ro_with_inhalo._local[1:-1, -1] == 2.) - assert np.all(f._data_ro_with_inhalo._local[:, 0] == 0.) + assert np.all(f._data_ro_with_inhalo[1:-1, -1] == 2.) + assert np.all(f._data_ro_with_inhalo[:, 0] == 0.) else: - assert np.all(f._data_ro_with_inhalo._local[1:-1, 0] == 1.) - assert np.all(f._data_ro_with_inhalo._local[:, -1] == 0.) - assert np.all(f._data_ro_with_inhalo._local[0] == 0.) - assert np.all(f._data_ro_with_inhalo._local[-1] == 0.) + assert np.all(f._data_ro_with_inhalo[1:-1, 0] == 1.) + assert np.all(f._data_ro_with_inhalo[:, -1] == 0.) + assert np.all(f._data_ro_with_inhalo[0] == 0.) + assert np.all(f._data_ro_with_inhalo[-1] == 0.) @pytest.mark.parallel(nprocs=2) def test_halo_exchange_bilateral_asymmetric(self): @@ -180,13 +185,13 @@ def test_halo_exchange_bilateral_asymmetric(self): glb_pos_map = grid.distributor.glb_pos_map if LEFT in glb_pos_map[y]: - assert np.all(f._data_ro_with_inhalo._local[2:-1, -1] == 2.) - assert np.all(f._data_ro_with_inhalo._local[:, 0:2] == 0.) + assert np.all(f._data_ro_with_inhalo[2:-1, -1] == 2.) + assert np.all(f._data_ro_with_inhalo[:, 0:2] == 0.) else: - assert np.all(f._data_ro_with_inhalo._local[2:-1, 0:2] == 1.) - assert np.all(f._data_ro_with_inhalo._local[:, -1] == 0.) - assert np.all(f._data_ro_with_inhalo._local[0:2] == 0.) - assert np.all(f._data_ro_with_inhalo._local[-1] == 0.) + assert np.all(f._data_ro_with_inhalo[2:-1, 0:2] == 1.) + assert np.all(f._data_ro_with_inhalo[:, -1] == 0.) + assert np.all(f._data_ro_with_inhalo[0:2] == 0.) + assert np.all(f._data_ro_with_inhalo[-1] == 0.) @pytest.mark.parallel(nprocs=4) def test_halo_exchange_quadrilateral(self): @@ -240,30 +245,31 @@ def test_halo_exchange_quadrilateral(self): glb_pos_map = grid.distributor.glb_pos_map if LEFT in glb_pos_map[x] and LEFT in glb_pos_map[y]: - assert np.all(f._data_ro_with_inhalo._local[0] == 0.) - assert np.all(f._data_ro_with_inhalo._local[:, 0] == 0.) - assert np.all(f._data_ro_with_inhalo._local[1:-1, -1] == 2.) - assert np.all(f._data_ro_with_inhalo._local[-1, 1:-1] == 3.) - assert f._data_ro_with_inhalo._local[-1, -1] == 4. + assert np.all(f._data_ro_with_inhalo[0] == 0.) + assert np.all(f._data_ro_with_inhalo[:, 0] == 0.) + assert np.all(f._data_ro_with_inhalo[1:-1, -1] == 2.) + assert np.all(f._data_ro_with_inhalo[-1, 1:-1] == 3.) + assert f._data_ro_with_inhalo[-1, -1] == 4. elif LEFT in glb_pos_map[x] and RIGHT in glb_pos_map[y]: - assert np.all(f._data_ro_with_inhalo._local[0] == 0.) - assert np.all(f._data_ro_with_inhalo._local[:, -1] == 0.) - assert np.all(f._data_ro_with_inhalo._local[1:-1, 0] == 1.) - assert np.all(f._data_ro_with_inhalo._local[-1, 1:-1] == 4.) - assert f._data_ro_with_inhalo._local[-1, 0] == 3. + assert np.all(f._data_ro_with_inhalo[0] == 0.) + assert np.all(f._data_ro_with_inhalo[:, -1] == 0.) + assert np.all(f._data_ro_with_inhalo[1:-1, 0] == 1.) + assert np.all(f._data_ro_with_inhalo[-1, 1:-1] == 4.) + assert f._data_ro_with_inhalo[-1, 0] == 3. elif RIGHT in glb_pos_map[x] and LEFT in glb_pos_map[y]: - assert np.all(f._data_ro_with_inhalo._local[-1] == 0.) - assert np.all(f._data_ro_with_inhalo._local[:, 0] == 0.) - assert np.all(f._data_ro_with_inhalo._local[1:-1, -1] == 4.) - assert np.all(f._data_ro_with_inhalo._local[0, 1:-1] == 1.) - assert f._data_ro_with_inhalo._local[0, -1] == 2. + assert np.all(f._data_ro_with_inhalo[-1] == 0.) + assert np.all(f._data_ro_with_inhalo[:, 0] == 0.) + assert np.all(f._data_ro_with_inhalo[1:-1, -1] == 4.) + assert np.all(f._data_ro_with_inhalo[0, 1:-1] == 1.) + assert f._data_ro_with_inhalo[0, -1] == 2. else: - assert np.all(f._data_ro_with_inhalo._local[-1] == 0.) - assert np.all(f._data_ro_with_inhalo._local[:, -1] == 0.) - assert np.all(f._data_ro_with_inhalo._local[1:-1, 0] == 3.) - assert np.all(f._data_ro_with_inhalo._local[0, 1:-1] == 2.) - assert f._data_ro_with_inhalo._local[0, 0] == 1. + assert np.all(f._data_ro_with_inhalo[-1] == 0.) + assert np.all(f._data_ro_with_inhalo[:, -1] == 0.) + assert np.all(f._data_ro_with_inhalo[1:-1, 0] == 3.) + assert np.all(f._data_ro_with_inhalo[0, 1:-1] == 2.) + assert f._data_ro_with_inhalo[0, 0] == 1. + @skipif_nompi @pytest.mark.parallel(nprocs=4) @pytest.mark.parametrize('shape,expected', [ ((15, 15), [((0, 8), (0, 8)), ((0, 8), (8, 15)), @@ -277,8 +283,8 @@ def test_local_indices(self, shape, expected): for i, j in zip(f.local_indices, expected[grid.distributor.myrank])) -@skipif_yask -@skipif_mpi +@skipif_nompi +@skipif_backend(['yask', 'ops']) class TestCodeGeneration(object): def test_iet_copy(self): @@ -365,8 +371,8 @@ def test_iet_update_halo(self): }""" -@skipif_yask -@skipif_mpi +@skipif_nompi +@skipif_backend(['yask', 'ops']) class TestSparseFunction(object): @pytest.mark.parallel(nprocs=4) @@ -455,8 +461,8 @@ def test_scatter_gather(self): assert np.all(sf.data == data[sf.local_indices]*2) -@skipif_yask -@skipif_mpi +@skipif_nompi +@skipif_backend(['yask', 'ops']) class TestOperatorSimple(object): @pytest.mark.parallel(nprocs=[2, 4, 8, 16, 32]) @@ -654,9 +660,26 @@ def test_haloupdate_not_requried(self): calls = FindNodes(Call).visit(op) assert len(calls) == 0 + @pytest.mark.parallel(nprocs=2) + def test_reapply_with_different_functions(self): + grid1 = Grid(shape=(30, 30, 30)) + f1 = Function(name='f', grid=grid1, space_order=4) + + op = Operator(Eq(f1, 1.)) + op.apply() + + grid2 = Grid(shape=(40, 40, 40)) + f2 = Function(name='f', grid=grid2, space_order=4) + + # Re-application + op.apply(f=f2) + + assert np.all(f1.data == 1.) + assert np.all(f2.data == 1.) + -@skipif_yask -@skipif_mpi +@skipif_nompi +@skipif_backend(['yask', 'ops']) class TestOperatorAdvanced(object): @pytest.mark.parallel(nprocs=[4]) @@ -947,6 +970,49 @@ def test_bcs_basic(self): assert np.all(u.data_ro_domain[0, :-thickness] == 1.) assert np.all(u.data_ro_domain[0, -thickness:] == range(2, thickness+2)) + @pytest.mark.parallel(nprocs=4) + def test_misc_dims(self): + """ + Test MPI in presence of Functions with mixed distributed/replicated + Dimensions, with only a strict subset of the Grid dimensions used. + """ + dx = Dimension(name='dx') + grid = Grid(shape=(4, 4)) + x, y = grid.dimensions + glb_pos_map = grid.distributor.glb_pos_map + time = grid.time_dim + + u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2, save=4) + c = Function(name='c', grid=grid, dimensions=(x, dx), shape=(4, 5)) + + step = Eq(u.forward, ( + u[time, x-2, y] * c[x, 0] + + u[time, x-1, y] * c[x, 1] + + u[time, x, y] * c[x, 2] + + u[time, x+1, y] * c[x, 3] + + u[time, x+2, y] * c[x, 4])) + + for i in range(4): + c.data[i, 0] = 1.0+i + c.data[i, 1] = 1.0+i + c.data[i, 2] = 3.0+i + c.data[i, 3] = 6.0+i + c.data[i, 4] = 5.0+i + + u.data[:] = 0.0 + u.data[0, 2, :] = 2.0 + + op = Operator(step) + + op(time_m=0, time_M=0) + + if LEFT in glb_pos_map[x]: + assert(np.all(u.data[1, 0, :] == 10.0)) + assert(np.all(u.data[1, 1, :] == 14.0)) + else: + assert(np.all(u.data[1, 2, :] == 10.0)) + assert(np.all(u.data[1, 3, :] == 8.0)) + @pytest.mark.parallel(nprocs=9) def test_nontrivial_operator(self): """ @@ -1029,25 +1095,76 @@ def test_nontrivial_operator(self): assert np.all(u.data_ro_domain[1] == 3) +@skipif_backend(['yask', 'ops']) class TestIsotropicAcoustic(object): """ - Test the acoustic wave model with MPI. + Test the isotropic acoustic wave equation with MPI. """ - # TODO: Cannot mark the following test as `xfail` since this marker - # doesn't cope well with the `parallel` mark. Leaving it commented out - # for the time being... - # @pytest.mark.parametrize('shape, kernel, space_order, nbpml', [ - # # 1 tests with varying time and space orders - # ((60, ), 'OT2', 4, 10), - # ]) - # @pytest.mark.parallel(nprocs=2) - # def test_adjoint_F(self, shape, kernel, space_order, nbpml): - # from test_adjoint import TestAdjoint - # TestAdjoint().test_adjoint_F('layers', shape, kernel, space_order, nbpml) + @pytest.mark.parametrize('shape,kernel,space_order,nbpml,save,Eu,Erec,Ev,Esrca', [ + ((60, ), 'OT2', 4, 10, False, 976.825, 9372.604, 18851836.075, 47002871.882), + ((60, 70), 'OT2', 8, 10, False, 351.217, 867.420, 405805.482, 239444.952), + ((60, 70, 80), 'OT2', 12, 10, False, 153.122, 205.902, 27484.635, 11736.917) + ]) + @pytest.mark.parallel(nprocs=[4, 8]) + def test_adjoint_F(self, shape, kernel, space_order, nbpml, save, + Eu, Erec, Ev, Esrca): + """ + Unlike `test_adjoint_F` in test_adjoint.py, here we explicitly check the norms + of all Operator-evaluated Functions. The numbers we check against are derived + "manually" from sequential runs of test_adjoint::test_adjoint_F + """ + t0 = 0.0 # Start time + tn = 500. # Final time + nrec = 130 # Number of receivers + spacing = 15. # Grid spacing + + # Create model from preset + model = demo_model(spacing=[spacing for _ in shape], dtype=np.float64, + space_order=space_order, shape=shape, nbpml=nbpml, + preset='layers-isotropic', ratio=3) + + # Derive timestepping from model spacing + dt = model.critical_dt * (1.73 if kernel == 'OT4' else 1.0) + time_range = TimeAxis(start=t0, stop=tn, step=dt) + + # Define source geometry (center of domain, just below surface) + src = RickerSource(name='src', grid=model.grid, f0=0.01, time_range=time_range) + src.coordinates.data[0, :] = np.array(model.domain_size) * .5 + src.coordinates.data[0, -1] = 30. + + # Define receiver geometry (same as source, but spread across x) + rec = Receiver(name='rec', grid=model.grid, time_range=time_range, npoint=nrec) + rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec) + if len(shape) > 1: + rec.coordinates.data[:, 1] = np.array(model.domain_size)[1] * .5 + rec.coordinates.data[:, -1] = 30. + + # Create solver object to provide relevant operators + solver = AcousticWaveSolver(model, source=src, receiver=rec, + kernel=kernel, space_order=space_order) + + # Create adjoint receiver symbol + srca = Receiver(name='srca', grid=model.grid, time_range=solver.source.time_range, + coordinates=solver.source.coordinates.data) + + # Run forward operator + rec, u, _ = solver.forward(save=save, rec=rec) + + assert np.isclose(norm(u), Eu, rtol=Eu*1.e-8) + assert np.isclose(norm(rec), Erec, rtol=Erec*1.e-8) + + # Run adjoint operator + srca, v, _ = solver.adjoint(rec=rec, srca=srca) + + assert np.isclose(norm(v), Ev, rtol=Eu*1.e-8) + assert np.isclose(norm(srca), Esrca, rtol=Erec*1.e-8) - pass + # Adjoint test: Verify matches closely + term1 = inner(srca, solver.source) + term2 = norm(rec)**2 + assert np.isclose((term1 - term2)/term1, 0., rtol=1.e-10) if __name__ == "__main__": @@ -1058,5 +1175,7 @@ class TestIsotropicAcoustic(object): # TestSparseFunction().test_ownership(((1., 1.), (1., 3.), (3., 1.), (3., 3.))) # TestSparseFunction().test_local_indices([(0.5, 0.5), (1.5, 2.5), (1.5, 1.5), (2.5, 1.5)], [[0.], [1.], [2.], [3.]]) # noqa # TestSparseFunction().test_scatter_gather() - TestOperatorAdvanced().test_nontrivial_operator() + # TestOperatorAdvanced().test_nontrivial_operator() # TestOperatorAdvanced().test_interpolation_dup() + TestIsotropicAcoustic().test_adjoint_F((60, 70, 80), 'OT2', 12, 10, False, + 153.122, 205.902, 27484.635, 11736.917) diff --git a/tests/test_operator.py b/tests/test_operator.py index daa38f2180..b39860aabb 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -1,19 +1,20 @@ from __future__ import absolute_import - -from conftest import EVAL, time, x, y, z, skipif_yask - +from conftest import EVAL, time, x, y, z import numpy as np import pytest - from devito import (clear_cache, Grid, Eq, Operator, Constant, Function, TimeFunction, SparseFunction, SparseTimeFunction, Dimension, error, SpaceDimension, - NODE, CELL) + NODE, CELL, configuration) from devito.ir.iet import (Expression, Iteration, ArrayCast, FindNodes, IsPerfectIteration, retrieve_iteration_tree) from devito.ir.support import Any, Backward, Forward from devito.symbolics import indexify, retrieve_indexed from devito.tools import flatten +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + def dimify(dimensions): assert isinstance(dimensions, str) @@ -25,11 +26,10 @@ def symbol(name, dimensions, value=0., shape=(3, 5), mode='function'): and "indexed" API.""" assert(mode in ['function', 'indexed']) s = Function(name=name, dimensions=dimensions, shape=shape) - s.data_allocated[:] = value + s.data_with_halo[:] = value return s.indexify() if mode == 'indexed' else s -@skipif_yask class TestCodeGen(object): @classmethod @@ -130,7 +130,6 @@ def test_multiple_steppers(self, expr, exp_uindices, exp_mods): for i in flatten(retrieve_indexed(i) for i in exprs)) -@skipif_yask class TestArithmetic(object): @classmethod @@ -303,7 +302,6 @@ def test_incs_same_lhs(self): assert np.all(u.data[:] == 3) -@skipif_yask class TestAllocation(object): @classmethod @@ -369,7 +367,6 @@ def test_staggered_time(self, stagg, ndim): assert f.data[index] == 2. -@skipif_yask class TestArguments(object): @classmethod @@ -421,7 +418,7 @@ def test_default_functions(self): 'x_size': 5, 'x_m': 0, 'x_M': 4, 'y_size': 6, 'y_m': 0, 'y_M': 5, 'z_size': 7, 'z_m': 0, 'z_M': 6, - 'f': f.data_allocated, 'g': g.data_allocated, + 'f': f.data_with_halo, 'g': g.data_with_halo, } self.verify_arguments(op.arguments(time=4), expected) exp_parameters = ['f', 'g', 'x_m', 'x_M', 'x_size', 'y_m', @@ -473,7 +470,7 @@ def test_override_function_size(self): 'x_size': 5, 'x_m': 0, 'x_M': 3, 'y_size': 6, 'y_m': 0, 'y_M': 4, 'z_size': 7, 'z_m': 0, 'z_M': 5, - 'g': g.data_allocated + 'g': g.data_with_halo } self.verify_arguments(arguments, expected) # Verify execution @@ -497,7 +494,7 @@ def test_override_function_subrange(self): 'x_size': 5, 'x_m': 1, 'x_M': 3, 'y_size': 6, 'y_m': 2, 'y_M': 4, 'z_size': 7, 'z_m': 3, 'z_M': 5, - 'g': g.data_allocated + 'g': g.data_with_halo } self.verify_arguments(arguments, expected) # Verify execution @@ -528,7 +525,7 @@ def test_override_timefunction_subrange(self): 'y_size': 6, 'y_m': 2, 'y_M': 4, 'z_size': 7, 'z_m': 3, 'z_M': 5, 'time_m': 1, 'time_M': 4, - 'f': f.data_allocated + 'f': f.data_with_halo } self.verify_arguments(arguments, expected) # Verify execution @@ -564,7 +561,7 @@ def test_override_function_data(self): assert (a2.data[:] == 6.).all() # Override with user-allocated numpy data - a3 = np.zeros_like(a.data_allocated) + a3 = np.zeros_like(a.data_with_halo) a3[:] = 4. op(a=a3) assert (a3[[slice(i.left, -i.right) for i in a._offset_domain]] == 7.).all() @@ -597,7 +594,7 @@ def test_override_timefunction_data(self): assert (a2.data[:] == 6.).all() # Override with user-allocated numpy data - a3 = np.zeros_like(a.data_allocated) + a3 = np.zeros_like(a.data_with_halo) a3[:] = 4. op(time_m=0, time=1, a=a3) assert (a3[[slice(i.left, -i.right) for i in a._offset_domain]] == 7.).all() @@ -791,8 +788,30 @@ def test_argument_no_shifting(self): assert (a.data[3:7, :] >= 2.).all() assert (a.data[8:, :] == 1.).all() + def test_argument_unknown(self): + """Check that Operators deal with unknown runtime arguments.""" + grid = Grid(shape=(11, 11)) + a = Function(name='a', grid=grid) + + op = Operator(Eq(a, a + a)) + try: + op.apply(b=3) + assert False + except ValueError: + # `b` means nothing to `op`, so we end up here + assert True + + try: + configuration['ignore-unknowns'] = True + op.apply(b=3) + assert True + except ValueError: + # we should not end up here as we're now ignoring unknown arguments + assert False + finally: + configuration['ignore-unknowns'] = configuration._defaults['ignore-unknowns'] + -@skipif_yask class TestDeclarator(object): @classmethod @@ -909,7 +928,6 @@ def test_stack_vector_temporaries(self, c_stack, e): return 0;""" in str(operator.ccode) -@skipif_yask class TestLoopScheduler(object): def test_consistency_coupled_wo_ofs(self, tu, tv, ti0, t0, t1): @@ -1141,10 +1159,10 @@ def test_equations_mixed_densedata_timedata(self, shape, dimensions): p_aux = Dimension(name='p_aux') b = Function(name='b', shape=shape + (10,), dimensions=dimensions + (p_aux,), space_order=2) - b.data_allocated[:] = 1.0 + b.data_with_halo[:] = 1.0 b2 = Function(name='b2', shape=(10,) + shape, dimensions=(p_aux,) + dimensions, space_order=2) - b2.data_allocated[:] = 1.0 + b2.data_with_halo[:] = 1.0 eqns = [Eq(a.forward, a.laplace + 1.), Eq(b, time*b*a + b)] eqns2 = [Eq(a.forward, a.laplace + 1.), @@ -1161,7 +1179,7 @@ def test_equations_mixed_densedata_timedata(self, shape, dimensions): # Verify both operators produce the same result op(time=10) - a.data_allocated[:] = 0. + a.data_with_halo[:] = 0. op2(time=10) for i in range(10): diff --git a/tests/test_ops.py b/tests/test_ops.py new file mode 100644 index 0000000000..de91c90cbf --- /dev/null +++ b/tests/test_ops.py @@ -0,0 +1,2 @@ +""" Empty file for now +""" diff --git a/tests/test_pickle.py b/tests/test_pickle.py index 348c4a7e60..5fa5ddbb99 100644 --- a/tests/test_pickle.py +++ b/tests/test_pickle.py @@ -1,15 +1,46 @@ +import pytest +from conftest import skipif_backend + import numpy as np from sympy import Symbol from examples.seismic import demo_model -from examples.seismic.source import TimeAxis, RickerSource +from examples.seismic.source import TimeAxis, RickerSource, Receiver from devito import (Constant, Eq, Function, TimeFunction, SparseFunction, Grid, - TimeDimension, SteppingDimension, Operator) + TimeDimension, SteppingDimension, Operator, configuration) +from devito.mpi.routines import MPIStatusObject, MPIRequestObject +from devito.profiling import Timer from devito.symbolics import IntDiv, ListInitializer, FunctionFromPointer import cloudpickle as pickle +pytestmark = pytest.mark.skipif(configuration['backend'] == 'ops', + reason="testing is currently restricted") + + +@pytest.fixture +def enable_mpi_codegen(request): + configuration['mpi'] = True + + def fin(): + configuration['mpi'] = False + + request.addfinalizer(fin) + + +def test_constant(): + c = Constant(name='c') + assert c.data == 0. + c.data = 1. + + pkl_c = pickle.dumps(c) + new_c = pickle.loads(pkl_c) + + # .data is initialized, so it should have been pickled too + assert np.all(c.data == 1.) + assert np.all(new_c.data == 1.) + def test_function(): grid = Grid(shape=(3, 3, 3)) @@ -49,6 +80,22 @@ def test_sparse_function(): assert sf.npoint == new_sf.npoint +def test_receiver(): + grid = Grid(shape=(3,)) + time_range = TimeAxis(start=0., stop=1000., step=0.1) + nreceivers = 3 + + rec = Receiver(name='rec', grid=grid, time_range=time_range, npoint=nreceivers, + coordinates=[(0.,), (1.,), (2.,)]) + rec.data[:] = 1. + + pkl_rec = pickle.dumps(rec) + new_rec = pickle.loads(pkl_rec) + + assert np.all(new_rec.data == 1) + assert np.all(new_rec.coordinates.data == [[0.], [1.], [2.]]) + + def test_symbolics(): a = Symbol('a') @@ -68,6 +115,17 @@ def test_symbolics(): assert li == new_li +def test_timers(): + """Pickling for Timers used in Operators for C-level profiling.""" + timer = Timer('timer', ['sec0', 'sec1']) + pkl_obj = pickle.dumps(timer) + new_obj = pickle.loads(pkl_obj) + assert new_obj.name == timer.name + assert new_obj.sections == timer.sections + assert new_obj.value._obj.sec0 == timer.value._obj.sec0 == 0.0 + assert new_obj.value._obj.sec1 == timer.value._obj.sec1 == 0.0 + + def test_operator_parameters(): grid = Grid(shape=(3, 3, 3)) f = Function(name='f', grid=grid) @@ -79,6 +137,22 @@ def test_operator_parameters(): pickle.loads(pkl_i) +def test_unjitted_operator(): + grid = Grid(shape=(3, 3, 3)) + f = Function(name='f', grid=grid) + + op = Operator(Eq(f, f + 1)) + + pkl_op = pickle.dumps(op) + new_op = pickle.loads(pkl_op) + + assert str(op) == str(new_op) + + +# With yask, broken padding in the generated code upon pickling, since +# in this test the data is allocated after generating the code. +# This is a symptom we need parametric padding +@skipif_backend(['yask']) def test_operator_function(): grid = Grid(shape=(3, 3, 3)) f = Function(name='f', grid=grid) @@ -89,10 +163,33 @@ def test_operator_function(): pkl_op = pickle.dumps(op) new_op = pickle.loads(pkl_op) + assert str(op) == str(new_op) + + new_op.apply(f=f) + assert np.all(f.data == 2) + + +def test_operator_function_w_preallocation(): + grid = Grid(shape=(3, 3, 3)) + f = Function(name='f', grid=grid) + f.data + + op = Operator(Eq(f, f + 1)) + op.apply() + + pkl_op = pickle.dumps(op) + new_op = pickle.loads(pkl_op) + + assert str(op) == str(new_op) + new_op.apply(f=f) assert np.all(f.data == 2) +# With yask, broken padding in the generated code upon pickling, since +# in this test the data is allocated after generating the code. +# This is a symptom we need parametric padding +@skipif_backend(['yask']) def test_operator_timefunction(): grid = Grid(shape=(3, 3, 3)) f = TimeFunction(name='f', grid=grid, save=3) @@ -103,10 +200,88 @@ def test_operator_timefunction(): pkl_op = pickle.dumps(op) new_op = pickle.loads(pkl_op) + assert str(op) == str(new_op) + + new_op.apply(time_m=1, time_M=1, f=f) + assert np.all(f.data[2] == 2) + + +def test_operator_timefunction_w_preallocation(): + grid = Grid(shape=(3, 3, 3)) + f = TimeFunction(name='f', grid=grid, save=3) + f.data + + op = Operator(Eq(f.forward, f + 1)) + op.apply(time=0) + + pkl_op = pickle.dumps(op) + new_op = pickle.loads(pkl_op) + + assert str(op) == str(new_op) + new_op.apply(time_m=1, time_M=1, f=f) assert np.all(f.data[2] == 2) +@skipif_backend(['yask']) +@pytest.mark.parallel(nprocs=[1]) +def test_mpi_objects(enable_mpi_codegen): + # Neighbours + grid = Grid(shape=(4, 4, 4)) + obj = grid.distributor._C_neighbours.obj + pkl_obj = pickle.dumps(obj) + new_obj = pickle.loads(pkl_obj) + assert obj.name == new_obj.name + assert obj.pname == new_obj.pname + assert obj.pfields == new_obj.pfields + assert obj.ptype == new_obj.ptype + + # Communicator + obj = grid.distributor._C_comm + pkl_obj = pickle.dumps(obj) + new_obj = pickle.loads(pkl_obj) + assert obj.name == new_obj.name + assert obj.dtype == new_obj.dtype + + # Status + obj = MPIStatusObject(name='status') + pkl_obj = pickle.dumps(obj) + new_obj = pickle.loads(pkl_obj) + assert obj.name == new_obj.name + assert obj.dtype == new_obj.dtype + + # Request + obj = MPIRequestObject(name='request') + pkl_obj = pickle.dumps(obj) + new_obj = pickle.loads(pkl_obj) + assert obj.name == new_obj.name + assert obj.dtype == new_obj.dtype + + +@skipif_backend(['yask']) +@pytest.mark.parallel(nprocs=[1]) +def test_mpi_operator(enable_mpi_codegen): + grid = Grid(shape=(4,)) + f = TimeFunction(name='f', grid=grid) + g = TimeFunction(name='g', grid=grid) + + # Using `sum` creates a stencil in `x`, which in turn will + # trigger the generation of code for MPI halo exchange + op = Operator(Eq(f.forward, f.sum() + 1)) + op.apply(time=2) + + pkl_op = pickle.dumps(op) + new_op = pickle.loads(pkl_op) + + assert str(op) == str(new_op) + + new_op.apply(time=2, f=g) + assert np.all(f.data[0] == [2., 3., 3., 3.]) + assert np.all(f.data[0] == [2., 3., 3., 3.]) + assert np.all(g.data[0] == f.data[0]) + assert np.all(g.data[1] == f.data[1]) + + def test_full_model(): shape = (50, 50, 50) diff --git a/tests/test_resample.py b/tests/test_resample.py index 8f19782857..6705f7117f 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -1,6 +1,11 @@ import numpy as np +import pytest +from devito import configuration from examples.seismic import TimeAxis, RickerSource, demo_model +pytestmark = pytest.mark.skipif(configuration['backend'] == 'ops', + reason="testing is currently restricted") + def test_resample(): diff --git a/tests/test_save.py b/tests/test_save.py index 999173373d..a0d4a31b9d 100644 --- a/tests/test_save.py +++ b/tests/test_save.py @@ -1,7 +1,13 @@ +import pytest import numpy as np -from conftest import skipif_yask -from devito import Buffer, Grid, Eq, Operator, TimeFunction, solve +from conftest import skipif_backend + +from devito import Buffer, Grid, Eq, Operator, TimeFunction, solve, configuration + + +pytestmark = pytest.mark.skipif(configuration['backend'] == 'ops', + reason="testing is currently restricted") def initial(nt, nx, ny): @@ -35,7 +41,7 @@ def run_simulation(save=False, dx=0.01, dy=0.01, a=0.5, timesteps=100): return u.data[timesteps - 1] -@skipif_yask +@skipif_backend(['yask']) def test_save(): assert(np.array_equal(run_simulation(True), run_simulation())) diff --git a/tests/test_symbol_caching.py b/tests/test_symbol_caching.py index 100f180d15..bb16cf3472 100644 --- a/tests/test_symbol_caching.py +++ b/tests/test_symbol_caching.py @@ -1,15 +1,15 @@ import weakref - import numpy as np import pytest -from conftest import skipif_yask - from devito import (Grid, Function, TimeFunction, SparseFunction, SparseTimeFunction, - Constant, Operator, Eq, Dimension, clear_cache) + Constant, Operator, Eq, Dimension, clear_cache, configuration) from devito.types import _SymbolCache, Scalar +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + -@skipif_yask @pytest.mark.parametrize('FunctionType', [Function, TimeFunction]) def test_cache_function_new(FunctionType): """Test that new u[x, y] instances don't cache""" @@ -22,7 +22,6 @@ def test_cache_function_new(FunctionType): assert np.allclose(u1.data, 2.) -@skipif_yask @pytest.mark.parametrize('FunctionType', [Function, TimeFunction]) def test_cache_function_same_indices(FunctionType): """Test caching of derived u[x, y] instance from derivative""" @@ -36,7 +35,6 @@ def test_cache_function_same_indices(FunctionType): assert np.allclose(u2.data, 6.) -@skipif_yask @pytest.mark.parametrize('FunctionType', [Function, TimeFunction]) def test_cache_function_different_indices(FunctionType): """Test caching of u[x + h, y] instance from derivative""" @@ -48,7 +46,6 @@ def test_cache_function_different_indices(FunctionType): assert np.allclose(u.data, u0.data) -@skipif_yask def test_cache_constant_new(): """Test that new u[x, y] instances don't cache""" u0 = Constant(name='u') @@ -59,7 +56,6 @@ def test_cache_constant_new(): assert u1.data == 2. -@skipif_yask def test_symbol_cache_aliasing(): """Test to assert that our aliasing cache isn't defeated by sympys non-aliasing symbol cache. @@ -100,7 +96,6 @@ def test_symbol_cache_aliasing(): assert u_ref() is None -@skipif_yask def test_symbol_cache_aliasing_reverse(): """Test to assert that removing he original u[x, y] instance does not impede our alisaing cache or leaks memory. @@ -135,7 +130,6 @@ def test_symbol_cache_aliasing_reverse(): assert u_ref() is None -@skipif_yask def test_clear_cache(nx=1000, ny=1000): grid = Grid(shape=(nx, ny), dtype=np.float64) clear_cache() @@ -151,7 +145,6 @@ def test_clear_cache(nx=1000, ny=1000): clear_cache() -@skipif_yask def test_cache_after_indexification(): """Test to assert that the SymPy cache retrieves the right Devito data object after indexification. @@ -166,7 +159,6 @@ def test_cache_after_indexification(): (i.indexify() + 1.).args[1].base.function.space_order -@skipif_yask def test_constant_hash(): """Test that different Constants have different hash value.""" c0 = Constant(name='c') @@ -175,7 +167,6 @@ def test_constant_hash(): assert hash(c0) != hash(c1) -@skipif_yask @pytest.mark.parametrize('FunctionType', [Function, TimeFunction]) def test_function_hash(FunctionType): """Test that different Functions have different hash value.""" @@ -191,7 +182,6 @@ def test_function_hash(FunctionType): assert hash(u0) != hash(u2) -@skipif_yask @pytest.mark.parametrize('FunctionType', [SparseFunction, SparseTimeFunction]) def test_sparse_function_hash(FunctionType): """Test that different Functions have different hash value.""" @@ -207,7 +197,6 @@ def test_sparse_function_hash(FunctionType): assert hash(u0) != hash(u2) -@skipif_yask def test_dimension_cache(): """ Test that :class:`Dimension`s with same name but different attributes do not @@ -231,7 +220,6 @@ def test_dimension_cache(): assert d2 is not d5 -@skipif_yask def test_operator_leakage_function(): """ Test to ensure that :class:`Operator` creation does not cause @@ -259,7 +247,6 @@ def test_operator_leakage_function(): assert w_op() is None -@skipif_yask def test_operator_leakage_sparse(): """ Test to ensure that :class:`Operator` creation does not cause diff --git a/tests/test_timestepping.py b/tests/test_timestepping.py index 4d04395169..e9795b46ce 100644 --- a/tests/test_timestepping.py +++ b/tests/test_timestepping.py @@ -1,9 +1,10 @@ import numpy as np - import pytest -from conftest import skipif_yask +from devito import Grid, Eq, Operator, TimeFunction, configuration -from devito import Grid, Eq, Operator, TimeFunction +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") @pytest.fixture @@ -35,7 +36,6 @@ def d(grid): return TimeFunction(name='d', grid=grid, time_order=2, save=6) -@skipif_yask def test_forward(a): a.data[0, :] = 1. Operator(Eq(a.forward, a + 1.))() @@ -43,7 +43,6 @@ def test_forward(a): assert np.allclose(a.data[i, :], 1. + i, rtol=1.e-12) -@skipif_yask def test_backward(b): b.data[-1, :] = 7. Operator(Eq(b.backward, b - 1.))() @@ -51,7 +50,6 @@ def test_backward(b): assert np.allclose(b.data[i, :], 2. + i, rtol=1.e-12) -@skipif_yask def test_forward_unroll(a, c, nt=5): """Test forward time marching with a buffered and an unrolled t""" a.data[0, :] = 1. @@ -63,7 +61,6 @@ def test_forward_unroll(a, c, nt=5): assert np.allclose(a.data[i, :], 1. + i, rtol=1.e-12) -@skipif_yask def test_forward_backward(a, b, nt=5): """Test a forward operator followed by a backward marching one""" a.data[0, :] = 1. @@ -77,7 +74,6 @@ def test_forward_backward(a, b, nt=5): assert np.allclose(b.data[i, :], 2. + i, rtol=1.e-12) -@skipif_yask def test_forward_backward_overlapping(a, b, nt=5): """ Test a forward operator followed by a backward one, but with @@ -94,7 +90,6 @@ def test_forward_backward_overlapping(a, b, nt=5): assert np.allclose(b.data[i, :], 2. + i, rtol=1.e-12) -@skipif_yask def test_loop_bounds_forward(d): """Test the automatic bound detection for forward time loops""" d.data[:] = 1. @@ -106,7 +101,6 @@ def test_loop_bounds_forward(d): assert np.allclose(d.data[i, :], 1. + i, rtol=1.e-12) -@skipif_yask def test_loop_bounds_backward(d): """Test the automatic bound detection for backward time loops""" d.data[:] = 5. diff --git a/tests/test_tools.py b/tests/test_tools.py index f0259fa5b5..4891ec9c90 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -1,12 +1,13 @@ import pytest -from conftest import skipif_yask - from sympy.abc import a, b, c, d, e - from devito.tools import toposort +from devito import configuration + +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") -@skipif_yask @pytest.mark.parametrize('elements, expected', [ ([[a, b, c], [c, d, e]], [a, b, c, d, e]), ([[e, d, c], [c, b, a]], [e, d, c, b, a]), diff --git a/tests/test_tti.py b/tests/test_tti.py index 29f3fc602e..bd063d3fac 100644 --- a/tests/test_tti.py +++ b/tests/test_tti.py @@ -1,16 +1,17 @@ import numpy as np import pytest -from conftest import skipif_yask from numpy import linalg - from devito import TimeFunction, configuration from devito.logger import log from examples.seismic import TimeAxis, RickerSource, Receiver, Model, demo_model from examples.seismic.acoustic import AcousticWaveSolver from examples.seismic.tti import AnisotropicWaveSolver +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") + -@skipif_yask @pytest.mark.parametrize('shape', [(120, 140), (120, 140, 150)]) @pytest.mark.parametrize('space_order', [4, 8]) @pytest.mark.parametrize('kernel', ['centered', 'shifted']) @@ -88,7 +89,6 @@ def test_tti(shape, space_order, kernel): assert np.isclose(res, 0.0, atol=1e-4) -@skipif_yask @pytest.mark.parametrize('shape', [(50, 60), (50, 60, 70)]) def test_tti_staggered(shape): spacing = [10. for _ in shape] diff --git a/tests/test_visitors.py b/tests/test_visitors.py index e1a8ae7345..83bb43e66f 100644 --- a/tests/test_visitors.py +++ b/tests/test_visitors.py @@ -1,12 +1,15 @@ import cgen as c import pytest -from conftest import skipif_yask - from devito.ir.equations import DummyEq from devito.ir.iet import (Block, Expression, Callable, FindSections, FindSymbols, IsPerfectIteration, Transformer, Conditional, printAST) from sympy import Mod, Eq +from devito import configuration + +pytestmark = pytest.mark.skipif(configuration['backend'] == 'yask' or + configuration['backend'] == 'ops', + reason="testing is currently restricted") @pytest.fixture(scope="module") @@ -64,7 +67,6 @@ def block4(exprs, iters, dims): return iters[0](Conditional(Eq(Mod(dims['i'], 2), 0), iters[1](exprs[0]))) -@skipif_yask def test_printAST(block1, block2, block3, block4): str1 = printAST(block1) assert str1 in """ @@ -105,7 +107,6 @@ def test_printAST(block1, block2, block3, block4): """ -@skipif_yask def test_create_cgen_tree(block1, block2, block3): assert str(Callable('foo', block1, 'void', ()).ccode) == """\ void foo() @@ -161,7 +162,6 @@ def test_create_cgen_tree(block1, block2, block3): }""" -@skipif_yask def test_find_sections(exprs, block1, block2, block3): finder = FindSections() @@ -182,7 +182,6 @@ def test_find_sections(exprs, block1, block2, block3): assert len(found[2]) == 1 -@skipif_yask def test_is_perfect_iteration(block1, block2, block3, block4): checker = IsPerfectIteration() @@ -204,7 +203,6 @@ def test_is_perfect_iteration(block1, block2, block3, block4): assert checker.visit(block4.nodes[0].then_body) is True -@skipif_yask def test_transformer_wrap(exprs, block1, block2, block3): """Basic transformer test that wraps an expression in comments""" line1 = '// This is the opening comment' @@ -223,7 +221,6 @@ def test_transformer_wrap(exprs, block1, block2, block3): assert "a[i] = a[i] + b[i] + 5.0F;" in newcode -@skipif_yask def test_transformer_replace(exprs, block1, block2, block3): """Basic transformer test that replaces an expression""" line1 = '// Replaced expression' @@ -240,7 +237,6 @@ def test_transformer_replace(exprs, block1, block2, block3): assert "a[i0] = a[i0] + b[i0] + 5.0F;" not in newcode -@skipif_yask def test_transformer_replace_function_body(block1, block2): """Create a Function and replace its body with another.""" args = FindSymbols().visit(block1) @@ -278,7 +274,6 @@ def test_transformer_replace_function_body(block1, block2): }""" -@skipif_yask def test_transformer_add_replace(exprs, block2, block3): """Basic transformer test that adds one expression and replaces another""" line1 = '// Replaced expression' @@ -299,7 +294,6 @@ def test_transformer_add_replace(exprs, block2, block3): assert "a[i0] = a[i0] + b[i0] + 5.0F;" not in newcode -@skipif_yask def test_nested_transformer(exprs, iters, block2): """When created with the kwarg ``nested=True``, a Transformer performs nested replacements. This test simultaneously replace an inner expression diff --git a/tests/test_yask.py b/tests/test_yask.py index 5155bb7342..1b0a09f6af 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -1,5 +1,6 @@ from sympy import cos import numpy as np +from cached_property import cached_property import pytest # noqa @@ -20,7 +21,8 @@ def setup_module(module): """Get rid of any YASK modules generated and JIT-compiled in previous runs. This is not strictly necessary for the tests, but it helps in keeping the - lib directory clean, which may be helpful for offline analysis.""" + lib directory clean, which may be helpful for offline analysis. + """ from devito.yask.wrappers import contexts # noqa contexts.dump() @@ -28,12 +30,12 @@ def setup_module(module): @pytest.fixture(autouse=True) def reset_isa(): """Force back to NO-SIMD after each test, as some tests may optionally - switch on SIMD.""" + switch on SIMD. + """ configuration['develop-mode'] = True class TestOperatorSimple(object): - """ Test execution of "toy" Operators through YASK. """ @@ -347,7 +349,6 @@ def test_repeated_op_calls(self): class TestOperatorAdvanced(object): - """ Test execution of non-trivial Operators through YASK. """ @@ -431,7 +432,6 @@ def test_subsampling(self): class TestIsotropicAcoustic(object): - """ Test the acoustic wave model through YASK. @@ -442,138 +442,139 @@ class TestIsotropicAcoustic(object): def setup_class(cls): clear_cache() - @pytest.fixture + @property def shape(self): return (60, 70, 80) - @pytest.fixture + @cached_property def nbpml(self): return 10 - @pytest.fixture + @cached_property def space_order(self): return 4 - @pytest.fixture + @cached_property def dtype(self): return np.float64 - @pytest.fixture - def model(self, space_order, shape, nbpml, dtype): - return demo_model(spacing=[15., 15., 15.], dtype=dtype, - space_order=space_order, shape=shape, nbpml=nbpml, - preset='layers-isotropic', ratio=3) + @cached_property + def model(self): + return demo_model(spacing=[15., 15., 15.], dtype=self.dtype, + space_order=self.space_order, shape=self.shape, + nbpml=self.nbpml, preset='layers-isotropic', ratio=3) - @pytest.fixture - def time_params(self, model): + @cached_property + def time_params(self): # Derive timestepping from model spacing t0 = 0.0 # Start time tn = 500. # Final time - dt = model.critical_dt + dt = self.model.critical_dt return t0, tn, dt - @pytest.fixture - def m(self, model): - return model.m + @cached_property + def m(self): + return self.model.m - @pytest.fixture - def damp(self, model): - return model.damp + @cached_property + def damp(self): + return self.model.damp - @pytest.fixture + @cached_property def kernel(self): return 'OT2' - @pytest.fixture - def u(self, model, space_order, kernel): - return TimeFunction(name='u', grid=model.grid, - space_order=space_order, time_order=2) + @cached_property + def u(self): + return TimeFunction(name='u', grid=self.model.grid, + space_order=self.space_order, time_order=2) - @pytest.fixture - def eqn(self, m, damp, u, kernel): - t = u.grid.stepping_dim - return iso_stencil(u, m, t.spacing, damp, kernel) + @cached_property + def eqn(self): + t = self.u.grid.stepping_dim + return iso_stencil(self.u, self.m, t.spacing, self.damp, self.kernel) - @pytest.fixture - def src(self, model, dtype): - t0, tn, dt = self.time_params(model) + @cached_property + def src(self): + t0, tn, dt = self.time_params time_range = TimeAxis(start=t0, stop=tn, step=dt) # Discretized time axis # Define source geometry (center of domain, just below surface) - src = RickerSource(name='src', grid=model.grid, f0=0.01, time_range=time_range, - dtype=dtype) - src.coordinates.data[0, :] = np.array(model.domain_size) * .5 + src = RickerSource(name='src', grid=self.model.grid, f0=0.01, + time_range=time_range, dtype=self.dtype) + src.coordinates.data[0, :] = np.array(self.model.domain_size) * .5 src.coordinates.data[0, -1] = 30. return src - @pytest.fixture - def rec(self, model, src, dtype): + @cached_property + def rec(self): nrec = 130 # Number of receivers - t0, tn, dt = self.time_params(model) + t0, tn, dt = self.time_params time_range = TimeAxis(start=t0, stop=tn, step=dt) - rec = Receiver(name='rec', grid=model.grid, + rec = Receiver(name='rec', grid=self.model.grid, time_range=time_range, - npoint=nrec, dtype=dtype) - rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec) - rec.coordinates.data[:, 1:] = src.coordinates.data[0, 1:] + npoint=nrec, dtype=self.dtype) + rec.coordinates.data[:, 0] = np.linspace(0., self.model.domain_size[0], num=nrec) + rec.coordinates.data[:, 1:] = self.src.coordinates.data[0, 1:] return rec - def test_acoustic_wo_src_wo_rec(self, model, eqn, m, damp, u): + def test_acoustic_wo_src_wo_rec(self): """ Test that the acoustic wave equation runs without crashing in absence of sources and receivers. """ - dt = model.critical_dt - u.data[:] = 0.0 - op = Operator(eqn, subs=model.spacing_map) + dt = self.model.critical_dt + self.u.data[:] = 0.0 + op = Operator(self.eqn, subs=self.model.spacing_map) assert 'run_solution' in str(op) - op.apply(u=u, m=m, damp=damp, time=10, dt=dt) + op.apply(u=self.u, m=self.m, damp=self.damp, time=10, dt=dt) - def test_acoustic_w_src_wo_rec(self, model, eqn, m, damp, u, src): + def test_acoustic_w_src_wo_rec(self): """ Test that the acoustic wave equation runs without crashing in absence of receivers. """ - dt = model.critical_dt - u.data[:] = 0.0 - eqns = eqn - eqns += src.inject(field=u.forward, expr=src * dt**2 / m) - op = Operator(eqns, subs=model.spacing_map) + dt = self.model.critical_dt + self.u.data[:] = 0.0 + eqns = self.eqn + eqns += self.src.inject(field=self.u.forward, expr=self.src * dt**2 / self.m) + op = Operator(eqns, subs=self.model.spacing_map) assert 'run_solution' in str(op) - op.apply(u=u, m=m, damp=damp, src=src, dt=dt) + op.apply(u=self.u, m=self.m, damp=self.damp, src=self.src, dt=dt) exp_u = 154.05 - assert np.isclose(np.linalg.norm(u.data[:]), exp_u, atol=exp_u*1.e-2) + assert np.isclose(np.linalg.norm(self.u.data[:]), exp_u, atol=exp_u*1.e-2) - def test_acoustic_w_src_w_rec(self, model, eqn, m, damp, u, src, rec): + def test_acoustic_w_src_w_rec(self): """ Test that the acoustic wave equation forward operator produces the correct results when running a 3D model also used in ``test_adjointA.py``. """ - dt = model.critical_dt - u.data[:] = 0.0 - eqns = eqn - eqns += src.inject(field=u.forward, expr=src * dt**2 / m) - eqns += rec.interpolate(expr=u) - op = Operator(eqns, subs=model.spacing_map) + dt = self.model.critical_dt + self.u.data[:] = 0.0 + eqns = self.eqn + eqns += self.src.inject(field=self.u.forward, expr=self.src * dt**2 / self.m) + eqns += self.rec.interpolate(expr=self.u) + op = Operator(eqns, subs=self.model.spacing_map) assert 'run_solution' in str(op) - op.apply(u=u, m=m, damp=damp, src=src, rec=rec, dt=dt) + op.apply(u=self.u, m=self.m, damp=self.damp, src=self.src, rec=self.rec, dt=dt) # The expected norms have been computed "by hand" looking at the output # of test_adjointA's forward operator w/o using the YASK backend. exp_u = 154.05 exp_rec = 212.15 - assert np.isclose(np.linalg.norm(u.data[:]), exp_u, atol=exp_u*1.e-2) - assert np.isclose(np.linalg.norm(rec.data.reshape(-1)), exp_rec, + assert np.isclose(np.linalg.norm(self.u.data[:]), exp_u, atol=exp_u*1.e-2) + assert np.isclose(np.linalg.norm(self.rec.data.reshape(-1)), exp_rec, atol=exp_rec*1.e-2) - def test_acoustic_adjoint(self, shape, kernel, space_order, nbpml): + def test_acoustic_adjoint(self): """ Full acoustic wave test, forward + adjoint operators """ from test_adjoint import TestAdjoint - TestAdjoint().test_adjoint_F('layers', shape, kernel, space_order, nbpml) + TestAdjoint().test_adjoint_F('layers', self.shape, self.kernel, + self.space_order, self.nbpml)