From e1e6c839ce268b0a0dede59cf165693f8ea32c55 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Sat, 12 Aug 2023 14:32:33 +0000 Subject: [PATCH 01/12] compiler: Workaround subprocess-induced leak at jit_compile --- devito/arch/compiler.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/devito/arch/compiler.py b/devito/arch/compiler.py index 56f0c7a6e5..4b7542eb70 100644 --- a/devito/arch/compiler.py +++ b/devito/arch/compiler.py @@ -10,7 +10,7 @@ import numpy.ctypeslib as npct from codepy.jit import compile_from_string -from codepy.toolchain import GCCToolchain, call_capture_output +from codepy.toolchain import GCCToolchain, call_capture_output as _call_capture_output from devito.arch import (AMDGPUX, Cpu64, M1, NVIDIAX, POWER8, POWER9, GRAVITON, INTELGPUX, get_nvidia_cc, check_cuda_runtime, @@ -19,7 +19,7 @@ from devito.logger import debug, warning, error from devito.parameters import configuration from devito.tools import (as_list, change_directory, filter_ordered, - memoized_func, memoized_meth, make_tempdir) + memoized_func, make_tempdir) __all__ = ['sniff_mpi_distro', 'compiler_registry'] @@ -123,6 +123,15 @@ def sniff_mpi_flags(mpicc='mpicc'): return compile_flags.split(), link_flags.split() +@memoized_func +def call_capture_output(cmd): + """ + Memoize calls to codepy's `call_capture_output` to avoid leaking memory due + to some prefork/subprocess voodoo. + """ + return _call_capture_output(cmd) + + class Compiler(GCCToolchain): """ Base class for all compiler classes. @@ -220,12 +229,16 @@ def __new_with__(self, **kwargs): def name(self): return self.__class__.__name__ - @memoized_meth + def get_version(self): + result, stdout, stderr = call_capture_output((self.cc, "--version")) + if result != 0: + raise RuntimeError(f"version query failed: {stderr}") + return stdout + def get_jit_dir(self): """A deterministic temporary directory for jit-compiled objects.""" return make_tempdir('jitcache') - @memoized_meth def get_codepy_dir(self): """A deterministic temporary directory for the codepy cache.""" return make_tempdir('codepy') @@ -729,9 +742,9 @@ def __init__(self, *args, **kwargs): def get_version(self): if configuration['mpi']: - cmd = [self.cc, "-cc=%s" % self.CC, "--version"] + cmd = (self.cc, "-cc=%s" % self.CC, "--version") else: - cmd = [self.cc, "--version"] + cmd = (self.cc, "--version") result, stdout, stderr = call_capture_output(cmd) if result != 0: raise RuntimeError(f"version query failed: {stderr}") From a3d3f4aded44a7abc418ab14cc5230b009aae76b Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Fri, 11 Aug 2023 09:41:26 +0000 Subject: [PATCH 02/12] compiler: Overhaul creation of AbstractFunctions --- devito/finite_differences/derivative.py | 14 ++- devito/operations/solve.py | 4 +- devito/symbolics/printer.py | 11 +- devito/types/array.py | 11 +- devito/types/basic.py | 132 ++++++++++++++---------- devito/types/dense.py | 36 ++++--- devito/types/sparse.py | 30 ++++-- devito/types/utils.py | 7 +- tests/test_caching.py | 63 ++++++----- 9 files changed, 195 insertions(+), 113 deletions(-) diff --git a/devito/finite_differences/derivative.py b/devito/finite_differences/derivative.py index e8db6e882e..6e57ac285a 100644 --- a/devito/finite_differences/derivative.py +++ b/devito/finite_differences/derivative.py @@ -106,8 +106,18 @@ def __new__(cls, expr, *dims, **kwargs): obj._deriv_order = orders if skip else DimensionTuple(*orders, getters=obj._dims) obj._side = kwargs.get("side") obj._transpose = kwargs.get("transpose", direct) - obj._ppsubs = as_tuple(frozendict(i) for i in - kwargs.get("subs", kwargs.get("_ppsubs", []))) + + ppsubs = kwargs.get("subs", kwargs.get("_ppsubs", [])) + processed = [] + if ppsubs: + for i in ppsubs: + try: + processed.append(frozendict(i)) + except AttributeError: + # E.g. `i` is a Transform object + processed.append(i) + obj._ppsubs = tuple(processed) + obj._x0 = frozendict(kwargs.get('x0', {})) return obj diff --git a/devito/operations/solve.py b/devito/operations/solve.py index 41ec329b83..12b19697c7 100644 --- a/devito/operations/solve.py +++ b/devito/operations/solve.py @@ -102,7 +102,7 @@ def _(expr): @singledispatch def factorize_target(expr, target): - return 1 if expr is target else 0 + return 1 if expr == target else 0 @factorize_target.register(Add) @@ -114,6 +114,7 @@ def _(expr, target): for a in expr.args: c += factorize_target(a, target) + return c @@ -125,4 +126,5 @@ def _(expr, target): c = 1 for a in expr.args: c *= a if not a.has(target) else factorize_target(a, target) + return c diff --git a/devito/symbolics/printer.py b/devito/symbolics/printer.py index c5f4e19e4c..c47ef95bfc 100644 --- a/devito/symbolics/printer.py +++ b/devito/symbolics/printer.py @@ -13,6 +13,7 @@ from devito.arch.compiler import AOMPCompiler from devito.symbolics.inspection import has_integer_args +from devito.types.basic import AbstractFunction __all__ = ['ccode'] @@ -44,10 +45,12 @@ def parenthesize(self, item, level, strict=False): return super().parenthesize(item, level, strict=strict) def _print_Function(self, expr): - # There exist no unknown Functions - if expr.func.__name__ not in self.known_functions: - self.known_functions[expr.func.__name__] = expr.func.__name__ - return super()._print_Function(expr) + if isinstance(expr, AbstractFunction): + return str(expr) + else: + if expr.func.__name__ not in self.known_functions: + self.known_functions[expr.func.__name__] = expr.func.__name__ + return super()._print_Function(expr) def _print_CondEq(self, expr): return "%s == %s" % (self._print(expr.lhs), self._print(expr.rhs)) diff --git a/devito/types/array.py b/devito/types/array.py index 8059471b4e..f12922f330 100644 --- a/devito/types/array.py +++ b/devito/types/array.py @@ -19,8 +19,15 @@ class ArrayBasic(AbstractFunction): is_ArrayBasic = True @classmethod - def __indices_setup__(cls, **kwargs): - return as_tuple(kwargs['dimensions']), as_tuple(kwargs['dimensions']) + def __indices_setup__(cls, *args, **kwargs): + dimensions = kwargs['dimensions'] + + if args: + indices = args + else: + indices = dimensions + + return as_tuple(dimensions), as_tuple(indices) @property def _C_name(self): diff --git a/devito/types/basic.py b/devito/types/basic.py index c450dd0335..10912259e6 100644 --- a/devito/types/basic.py +++ b/devito/types/basic.py @@ -727,7 +727,7 @@ def __subfunc_setup__(cls, *args, **kwargs): return [] -class AbstractFunction(sympy.Function, Basic, Cached, Pickable, Evaluable): +class AbstractFunction(sympy.Function, Basic, Pickable, Evaluable): """ Base class for tensor symbols, cached by both SymPy and Devito. It inherits @@ -797,54 +797,39 @@ class AbstractFunction(sympy.Function, Basic, Cached, Pickable, Evaluable): True if data is allocated as a single, contiguous chunk of memory. """ - __rkwargs__ = ('name', 'dtype', 'grid', 'halo', 'padding', 'alias') - - @classmethod - def _cache_key(cls, *args, **kwargs): - return cls, args - - @classmethod - def _cache_get(cls, *args, **kwargs): - # Is the object already in cache (e.g., f(x), f(x+1)) ? - key = cls._cache_key(*args, **kwargs) - obj = super()._cache_get(key) - if obj is not None: - return obj - - # Does the base object exist at least (e.g. f(x))? - obj = super()._cache_get(cls) - if obj is not None: - options = kwargs.get('options', {'evaluate': False}) - with sympy_mutex: - newobj = sympy.Function.__new__(cls, *args, **options) - newobj.__init_cached__(obj) - Cached.__init__(newobj, key) - return newobj - - # Not in cache - return None + __rkwargs__ = ('name', 'dtype', 'grid', 'halo', 'padding', 'alias', 'function') def __new__(cls, *args, **kwargs): - # Only perform a construction if the object isn't in cache - obj = cls._cache_get(*args, **kwargs) - if obj is not None: - return obj - - options = kwargs.get('options', {'evaluate': False}) - # Preprocess arguments args, kwargs = cls.__args_setup__(*args, **kwargs) - # Not in cache. Create a new Function via sympy.Function - name = kwargs.get('name') - dimensions, indices = cls.__indices_setup__(**kwargs) + # Extract the `indices`, as perhaps they're explicitly provided + dimensions, indices = cls.__indices_setup__(*args, **kwargs) - # Create new, unique type instance from cls and the symbol name - newcls = type(name, (cls,), dict(cls.__dict__)) + # If it's an alias or simply has a different name, ignore `function`. + # These cases imply the construction of a new AbstractFunction off + # an existing one! This relieves the pressure on the caller by not + # requiring `function=None` explicitly at rebuild + name = kwargs.get('name') + alias = kwargs.get('alias') + function = kwargs.get('function') + if alias or (function and function.name != name): + function = kwargs['function'] = None + + # If same name/indices and `function` isn't None, then it's + # definitely a reconstruction + if function is not None and \ + function.name == name and \ + function.indices == indices: + # Special case: a syntactically identical alias of `function`, so + # let's just return `function` itself + return function - # Create the new Function object and invoke __init__ with sympy_mutex: - newobj = sympy.Function.__new__(newcls, *indices, **options) + # Go straight through Basic, thus bypassing caching and machinery + # in sympy.Application/Function that isn't really necessary + # AbstractFunctions are unique by construction! + newobj = sympy.Basic.__new__(cls, *indices) # Initialization. The following attributes must be available # when executing __init_finalize__ @@ -854,14 +839,10 @@ def __new__(cls, *args, **kwargs): newobj._dtype = cls.__dtype_setup__(**kwargs) newobj.__init_finalize__(*args, **kwargs) - # All objects cached on the AbstractFunction `newobj` keep a reference - # to `newobj` through the `function` field. Thus, all indexified - # object will point to `newobj`, the "actual Function". - newobj.function = newobj - - # Store new instance in symbol cache - key = (newcls, indices) - Cached.__init__(newobj, key, newcls) + # All objects created off an existing AbstractFunction `f` (e.g., + # via .func, or .subs, such as `f(x + 1)`) keep a reference to `f` + # through the `function` field + newobj.function = function or newobj return newobj @@ -869,6 +850,36 @@ def __init__(self, *args, **kwargs): # no-op, the true init is performed by __init_finalize__ pass + def __str__(self): + return "%s(%s)" % (self.name, ', '.join(str(i) for i in self.indices)) + + __repr__ = __str__ + + def _sympystr(self, printer): + return str(self) + + def __eq__(self, other): + try: + return (self.function is other.function and + self.indices == other.indices) + except AttributeError: + # `other` not even an AbstractFunction + return False + + __hash__ = sympy.Function.__hash__ + + def _hashable_content(self): + return super()._hashable_content() + (id(self.function), self.indices) + + @sympy.cacheit + def sort_key(self, order=None): + # Ensure that `f(x)` appears before `g(x)` + # With the legacy caching framework this wasn't necessary because + # the function name was already encoded in the class_key + class_key, args, exp, coeff = super().sort_key(order=order) + args = (len(args[1]) + 1, (self.name,) + args[1]) + return class_key, args, exp, coeff + def __init_finalize__(self, *args, **kwargs): # Setup halo and padding regions self._is_halo_dirty = False @@ -888,8 +899,6 @@ def __init_finalize__(self, *args, **kwargs): # routines is applied to several actual DiscreteFunctions self._alias = kwargs.get('alias', False) - __hash__ = Cached.__hash__ - @classmethod def __args_setup__(cls, *args, **kwargs): """ @@ -902,7 +911,7 @@ def __args_setup__(cls, *args, **kwargs): return args, kwargs @classmethod - def __indices_setup__(cls, **kwargs): + def __indices_setup__(cls, *args, **kwargs): """Extract the object indices from ``kwargs``.""" return (), () @@ -1240,13 +1249,24 @@ def __getitem__(self, index): """Shortcut for ``self.indexed[index]``.""" return self.indexed[index] + # Reconstruction support + func = Pickable._rebuild + # Pickling support __reduce_ex__ = Pickable.__reduce_ex__ - # Reconstruction support - @property - def _rcls(self): - return self.__class__.__base__ + def __getnewargs_ex__(self): + args, kwargs = super().__getnewargs_ex__() + + # Since f(x) stashes a reference to itself f(x), we must drop it here + # or we'll end up with infinite recursion while attempting to serialize + # the object. Upon unpickling, we'll then get `function=None`, which is + # perfectly fine based on how `__new__`, and in particular the + # initialization of the `.function` attribute, works + if self is kwargs.get('function'): + kwargs.pop('function') + + return args, kwargs # Extended SymPy hierarchy follows, for essentially two reasons: diff --git a/devito/types/dense.py b/devito/types/dense.py index 5d50e22d9d..6bb5d1a34f 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -59,7 +59,7 @@ class DiscreteFunction(AbstractFunction, ArgProvider, Differentiable): __rkwargs__ = AbstractFunction.__rkwargs__ + ('staggered', 'initializer') - def __init_finalize__(self, *args, **kwargs): + def __init_finalize__(self, *args, function=None, **kwargs): # Staggering metadata self._staggered = self.__staggered_setup__(**kwargs) @@ -72,12 +72,21 @@ def __init_finalize__(self, *args, **kwargs): if self._coefficients not in ('standard', 'symbolic'): raise ValueError("coefficients must be `standard` or `symbolic`") - # Data-related properties and data initialization + # Data-related properties self._data = None self._first_touch = kwargs.get('first_touch', configuration['first-touch']) self._allocator = kwargs.get('allocator') or default_allocator() + + # Data initialization initializer = kwargs.get('initializer') - if initializer is None or callable(initializer) or self.alias: + if self.alias: + self._initializer = None + elif function is not None: + # An object derived from user-level AbstractFunction (e.g., + # `f(x+1)`), so we just copy the reference to the original data + self._initializer = None + self._data = function._data + elif initializer is None or callable(initializer) or self.alias: # Initialization postponed until the first access to .data self._initializer = initializer elif isinstance(initializer, (np.ndarray, list, tuple)): @@ -96,14 +105,6 @@ def __init_finalize__(self, *args, **kwargs): raise ValueError("`initializer` must be callable or buffer, not %s" % type(initializer)) - def __eq__(self, other): - # The only possibility for two DiscreteFunctions to be considered equal - # is that they are indeed the same exact object - return self is other - - def __hash__(self): - return id(self) - _subs = Differentiable._subs def _allocate_memory(func): @@ -1018,7 +1019,7 @@ def _eval_at(self, func): return self @classmethod - def __indices_setup__(cls, **kwargs): + def __indices_setup__(cls, *args, **kwargs): grid = kwargs.get('grid') dimensions = kwargs.get('dimensions') if grid is None: @@ -1027,6 +1028,9 @@ def __indices_setup__(cls, **kwargs): elif dimensions is None: dimensions = grid.dimensions + if args: + return tuple(dimensions), tuple(args) + # Staggered indices staggered = kwargs.get("staggered", None) if staggered in [CELL, NODE]: @@ -1324,9 +1328,10 @@ def __fd_setup__(self): to=self.time_order) @classmethod - def __indices_setup__(cls, **kwargs): + def __indices_setup__(cls, *args, **kwargs): dimensions = kwargs.get('dimensions') staggered = kwargs.get('staggered') + if dimensions is None: save = kwargs.get('save') grid = kwargs.get('grid') @@ -1338,7 +1343,10 @@ def __indices_setup__(cls, **kwargs): raise TypeError("`time_dim` must be a time dimension") dimensions = list(Function.__indices_setup__(**kwargs)[0]) dimensions.insert(cls._time_position, time_dim) - return Function.__indices_setup__(dimensions=dimensions, staggered=staggered) + + return Function.__indices_setup__( + *args, dimensions=dimensions, staggered=staggered + ) @classmethod def __shape_setup__(cls, **kwargs): diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 8eda411e25..3b5678a665 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -57,11 +57,17 @@ def __fd_setup__(self): return generate_fd_shortcuts(self.dimensions, self.space_order) @classmethod - def __indices_setup__(cls, **kwargs): + def __indices_setup__(cls, *args, **kwargs): dimensions = as_tuple(kwargs.get('dimensions')) if not dimensions: dimensions = (Dimension(name='p_%s' % kwargs["name"]),) - return dimensions, dimensions + + if args: + indices = args + else: + indices = dimensions + + return dimensions, indices @classmethod def __shape_setup__(cls, **kwargs): @@ -319,12 +325,18 @@ def time_dim(self): return self._time_dim @classmethod - def __indices_setup__(cls, **kwargs): + def __indices_setup__(cls, *args, **kwargs): dimensions = as_tuple(kwargs.get('dimensions')) if not dimensions: dimensions = (kwargs['grid'].time_dim, Dimension(name='p_%s' % kwargs["name"])) - return dimensions, dimensions + + if args: + indices = args + else: + indices = dimensions + + return dimensions, indices @classmethod def __shape_setup__(cls, **kwargs): @@ -1582,7 +1594,7 @@ def inject(self, field, expr, offset=0, u_t=None, p_t=None): return out @classmethod - def __indices_setup__(cls, **kwargs): + def __indices_setup__(cls, *args, **kwargs): """ Return the default dimension indices for a given data shape. """ @@ -1590,7 +1602,13 @@ def __indices_setup__(cls, **kwargs): if dimensions is None: dimensions = (kwargs['grid'].time_dim, Dimension( name='p_%s' % kwargs["name"])) - return dimensions, dimensions + + if args: + indices = args + else: + indices = dimensions + + return dimensions, indices @classmethod def __shape_setup__(cls, **kwargs): diff --git a/devito/types/utils.py b/devito/types/utils.py index 5d1d79f9c6..81daf1eb5c 100644 --- a/devito/types/utils.py +++ b/devito/types/utils.py @@ -42,9 +42,12 @@ class CtypesFactory(object): @classmethod def generate(cls, pname, pfields): - dtype = POINTER(type(pname, (Structure,), {'_fields_': pfields})) key = (pname, tuple(pfields)) - return cls.cache.setdefault(key, dtype) + try: + return cls.cache[key] + except KeyError: + dtype = POINTER(type(pname, (Structure,), {'_fields_': pfields})) + return cls.cache.setdefault(key, dtype) class HierarchyLayer(object): diff --git a/tests/test_caching.py b/tests/test_caching.py index 776cb9d2fe..408a49453e 100644 --- a/tests/test_caching.py +++ b/tests/test_caching.py @@ -263,12 +263,16 @@ def test_function_duplicates(self, FunctionType): # Make sure ub is u0 assert ub is u0 assert hash(ub) == hash(u0) - # Three new cache entries: u, u(t,x,y), u(t, x+h_x, y) - ncreated = 3 + # With the legacy caching model we would have created three + # entries: u, u(t,x,y), u(t, x+h_x, y) + # But now Devito doesn't cache AbstractFunctions anymore! + ncreated = 0 assert len(_SymbolCache) == _cache_size + ncreated - # shift again, no new entry should be created + # With the legacy caching model identical shifts such as two + # `u(x + h_x, y)` would have been the same physical object. Now they + # are two distinct objects, both uncached uf2 = ub.subs({x: x + x.spacing}) - assert uf is uf2 + assert uf is not uf2 assert len(_SymbolCache) == _cache_size + ncreated def test_symbols(self): @@ -546,8 +550,11 @@ def test_clear_cache(self, operate_on_empty_cache, nx=1000, ny=1000): assert(len(_SymbolCache) == cache_size) Function(name='u', grid=grid, space_order=2) - # Both u and u(inds) added to cache - assert(len(_SymbolCache) == cache_size + 2) + # With the legacy caching model we would have created two new + # entries in the cache, u and u(inds) + # But now Devito doesn't cache AbstractFunctions anymore! + ncreated = 0 + assert(len(_SymbolCache) == cache_size + ncreated) clear_cache() @@ -556,13 +563,16 @@ def test_clear_cache_with_Csymbol(self, operate_on_empty_cache, nx=1000, ny=1000 cache_size = len(_SymbolCache) u = Function(name='u', grid=grid, space_order=2) - # Both u and u(inds) added to cache - assert(len(_SymbolCache) == cache_size + 2) + # With the legacy caching model we would have created two new + # entries in the cache, u and u(inds) + # But now Devito doesn't cache AbstractFunctions anymore! + ncreated = 0 + assert(len(_SymbolCache) == cache_size + ncreated) u._C_symbol # Cache size won't change since _C_symbol isn't cached by devito to # avoid circular references in the cache - assert(len(_SymbolCache) == cache_size + 2) + assert(len(_SymbolCache) == cache_size + ncreated) def test_clear_cache_with_alive_symbols(self, operate_on_empty_cache, nx=1000, ny=1000): @@ -617,8 +627,10 @@ def test_sparse_function(self, operate_on_empty_cache): u = SparseFunction(name='u', grid=grid, npoint=1, nt=10) - # created: u, u(inds), p_u, h_p_u, u_coords, u_coords(inds), d, h_d - ncreated = 8 + # created: p_u, h_p_u, d, h_d + # With the legacy caching model also u, u(inds), u_coords, and + # u_coords(inds) would have been added to the cache; not anymore! + ncreated = 4 assert len(_SymbolCache) == cur_cache_size + ncreated cur_cache_size = len(_SymbolCache) @@ -642,11 +654,12 @@ def test_sparse_function(self, operate_on_empty_cache): del u del i clear_cache() - # At this point, not all children objects have been cleared. In particular, the - # ii_u_* Symbols are still alive, as well as p_u and h_p_u. This is because - # in the first clear_cache they were still referenced by their "parent" objects - # (e.g., ii_u_* by ConditionalDimensions, through `condition`) - assert len(_SymbolCache) == init_cache_size + 8 + # At this point, not all children objects have been cleared. In + # particular, the ii_u_* Symbols are still alive, as well as p_u and + # h_p_u. This is because in the first clear_cache they were still + # referenced by their "parent" objects (e.g., ii_u_* by + # ConditionalDimensions, through `condition`) + assert len(_SymbolCache) == init_cache_size + 6 clear_cache() # Now we should be back to the original state assert len(_SymbolCache) == init_cache_size @@ -795,17 +808,15 @@ def test_solve(self, operate_on_empty_cache): del eqn del grid - # We only deleted `u`, however we also cache shifted version created by the - # finite difference (u.dt, u.dx2). In this case we have three extra references - # to u(t + dt), u(x - h_x) and u(x + h_x) that have to be cleared. - # Then `u` points to the various Dimensions, the Dimensions point to the various - # spacing symbols, hence, we need four sweeps to clear up the cache. - assert len(_SymbolCache) == 14 + # We deleted `u`. + # With the legacy caching model, we would also have cache shifted versions + # created by the finite difference (u.dt, u.dx2). We would have had + # three extra references to u(t + dt), u(x - h_x) and u(x + h_x). + # But this is not the case anymore! + assert len(_SymbolCache) == 11 clear_cache() - assert len(_SymbolCache) == 9 + assert len(_SymbolCache) == 8 clear_cache() - assert len(_SymbolCache) == 3 - clear_cache() - assert len(_SymbolCache) == 1 + assert len(_SymbolCache) == 2 clear_cache() assert len(_SymbolCache) == 0 From 853ff2676267ee65a83099483855a88c51499aaf Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 16 Aug 2023 15:15:42 +0000 Subject: [PATCH 03/12] sympy: Monkey-patch together(), temporarily --- devito/__init__.py | 3 + devito/mpatches/.__init__.py.swp | Bin 0 -> 12288 bytes devito/mpatches/.rationaltools.py.swp | Bin 0 -> 12288 bytes devito/mpatches/__init__.py | 1 + devito/mpatches/rationaltools.py | 95 ++++++++++++++++++++++++++ 5 files changed, 99 insertions(+) create mode 100644 devito/mpatches/.__init__.py.swp create mode 100644 devito/mpatches/.rationaltools.py.swp create mode 100644 devito/mpatches/__init__.py create mode 100644 devito/mpatches/rationaltools.py diff --git a/devito/__init__.py b/devito/__init__.py index 1ba848d621..535b93ecbb 100644 --- a/devito/__init__.py +++ b/devito/__init__.py @@ -38,6 +38,9 @@ from devito.mpi.routines import mpi_registry from devito.operator import profiler_registry, operator_registry +# Apply monkey-patching while we wait for our patches to be upstreamed and released +from devito.mpatches import * # noqa + from ._version import get_versions # noqa __version__ = get_versions()['version'] diff --git a/devito/mpatches/.__init__.py.swp b/devito/mpatches/.__init__.py.swp new file mode 100644 index 0000000000000000000000000000000000000000..48f848b3e8678b2eef767d23fc7e5bf397b8c404 GIT binary patch literal 12288 zcmeI%K}*9h6u|K}@1p1zhz^7a>zWDTSMcC*9xd70VW4el8wOtVgYCL+;m-R8yW4i| zScIK6_&*4b@RCRJyM;`{*XI}Y6pW>gOOe%P9ls6kWP2epPxYg23r(>*k}4T4bQtRH zC2i$v*Kiak3+vrJWrgvx_oNy{kz^+5I_I{1Q0*~fm0RmMNhu#!=L?Kpx2R4^?sZi0R#|0009ILKmY**5I_I{1Vl_M zvaJ95f9lcv|L@=b?HnKk5I_I{1Q0*~0R#|0009ILICg<~=`t0R##@)0)H|0}%4UTt Ny}DPbuX6WcNmVs>m1OB*{vEPB;fM0;qz`H;SBtRGVX**-5fTO^Tz?auB_A2lS@F;LS z@Ws`Py$QSlJPF(ce03FLF9NH;-N0U8CvXSw&y|dQ2Yd^h2Hpam2d2OTxEELgJm4bW zn=2Um8h96Y4bZ?*;0Q1PZUJrve!ZNr&w$d1SsGh-~iAC z{=O8k0v`ZpfLDO$fED0Y;0JL0J@7H`5pW849e53R9ykb49-g_Fv4??&fSZB8>wlcN z0PE|2?6NMCtn`eUNHdi=-^c4)W=FOd_A(Wh+FSS1iY~XFw6(-pYB&10xr)u7BC;x+ zhQH@8>0Zat)h9}_I2xaWGB*esai?+<%~xuzO(#Tm(fi6NQBKJBXcoA9M=`) zhCEfqa@bwyraIG6s*Om4s?z%WfY4FkWI{RKK&%_sv13P-iO#GCuRAE`U3}2B_q+G@ zX5FQkb#njy{qsK887z_L9*mOklAB7>^`h#`8fyMa$$AwF1!TY~DZ8YV#0kxR4Q34E8m6BQ%4+R=5C%;~T?oJ7;`f6Z}zlEQJ00<0RZ z7zM61^I4yleSTtYb|F3nvO-;TWq8sH zXmP|tdY;J~T|40-p2$%yqQI(|s^fA^=6#qHW-+f0YDnrdS1EQJa5GYI6z~M602S&! z*UIY9iVue(6S?Ws^IW4C4!NGHA`a_bg96Pte4C1NB14ZKQbt-G2TFuE1{1X8`lpWJ`5y-GxF+etGuFTgY7n&zUFs1996jzMeB-B}K66|qe@n2{wdq4mH literal 0 HcmV?d00001 diff --git a/devito/mpatches/__init__.py b/devito/mpatches/__init__.py new file mode 100644 index 0000000000..834d6078d0 --- /dev/null +++ b/devito/mpatches/__init__.py @@ -0,0 +1 @@ +from .rationaltools import * # noqa diff --git a/devito/mpatches/rationaltools.py b/devito/mpatches/rationaltools.py new file mode 100644 index 0000000000..4e4cea1539 --- /dev/null +++ b/devito/mpatches/rationaltools.py @@ -0,0 +1,95 @@ +"""Tools for manipulation of rational expressions. """ + +import importlib + +import sympy +from sympy.core import Basic, Add, sympify +from sympy.core.exprtools import gcd_terms +from sympy.utilities import public +from sympy.utilities.iterables import iterable + +__all__ = [] + + +@public +def together(expr, deep=False, fraction=True): + """ + Denest and combine rational expressions using symbolic methods. + + This function takes an expression or a container of expressions + and puts it (them) together by denesting and combining rational + subexpressions. No heroic measures are taken to minimize degree + of the resulting numerator and denominator. To obtain completely + reduced expression use :func:`~.cancel`. However, :func:`~.together` + can preserve as much as possible of the structure of the input + expression in the output (no expansion is performed). + + A wide variety of objects can be put together including lists, + tuples, sets, relational objects, integrals and others. It is + also possible to transform interior of function applications, + by setting ``deep`` flag to ``True``. + + By definition, :func:`~.together` is a complement to :func:`~.apart`, + so ``apart(together(expr))`` should return expr unchanged. Note + however, that :func:`~.together` uses only symbolic methods, so + it might be necessary to use :func:`~.cancel` to perform algebraic + simplification and minimize degree of the numerator and denominator. + + Examples + ======== + + >>> from sympy import together, exp + >>> from sympy.abc import x, y, z + + >>> together(1/x + 1/y) + (x + y)/(x*y) + >>> together(1/x + 1/y + 1/z) + (x*y + x*z + y*z)/(x*y*z) + + >>> together(1/(x*y) + 1/y**2) + (x + y)/(x*y**2) + + >>> together(1/(1 + 1/x) + 1/(1 + 1/y)) + (x*(y + 1) + y*(x + 1))/((x + 1)*(y + 1)) + + >>> together(exp(1/x + 1/y)) + exp(1/y + 1/x) + >>> together(exp(1/x + 1/y), deep=True) + exp((x + y)/(x*y)) + + >>> together(1/exp(x) + 1/(x*exp(x))) + (x + 1)*exp(-x)/x + + >>> together(1/exp(2*x) + 1/(x*exp(3*x))) + (x*exp(x) + 1)*exp(-3*x)/x + + """ + def _together(expr): + if isinstance(expr, Basic): + if expr.is_Atom or (expr.is_Function and not deep): + return expr + elif expr.is_Add: + return gcd_terms(list(map(_together, Add.make_args(expr))), + fraction=fraction) + elif expr.is_Pow: + base = _together(expr.base) + + if deep: + exp = _together(expr.exp) + else: + exp = expr.exp + + return expr.func(base, exp) + else: + return expr.func(*[_together(arg) for arg in expr.args]) + elif iterable(expr): + return expr.__class__([_together(ex) for ex in expr]) + + return expr + + return _together(sympify(expr)) + + +# Apply the monkey patch +simplify = importlib.import_module(sympy.simplify.__module__) +simplify.together = together From d405d6d84b4835551d8bccd80ca4c49eb918999b Mon Sep 17 00:00:00 2001 From: mloubout Date: Thu, 17 Aug 2023 08:36:56 -0400 Subject: [PATCH 04/12] api: revamp solve to avoid xreplace and fix Derivative subs --- devito/finite_differences/derivative.py | 47 +++++++++++++++++-------- devito/operations/solve.py | 30 ++++++++++------ devito/symbolics/queries.py | 5 +++ devito/symbolics/search.py | 10 ++++-- tests/test_derivatives.py | 6 +++- tests/test_symbolics.py | 9 ++--- 6 files changed, 75 insertions(+), 32 deletions(-) diff --git a/devito/finite_differences/derivative.py b/devito/finite_differences/derivative.py index 6e57ac285a..0d559aba36 100644 --- a/devito/finite_differences/derivative.py +++ b/devito/finite_differences/derivative.py @@ -217,27 +217,43 @@ def _new_from_self(self, **kwargs): def func(self, expr, *args, **kwargs): return self._new_from_self(expr=expr, **kwargs) - def subs(self, *args, **kwargs): - """ - Bypass sympy.Subs as Devito has its own lazy evaluation mechanism. - """ - # Check if we are calling subs(self, old, new, **hint) in which case - # return the standard substitution. Need to check `==` rather than `is` - # because a new derivative could be created i.e `f.dx.subs(f.dx, y)` - if len(args) == 2 and args[0] == self: - return args[1] - try: - rules = dict(*args) - except TypeError: - rules = dict((args,)) - kwargs.pop('simultaneous', None) - return self.xreplace(rules, **kwargs) + def _subs(self, old, new, **hints): + # Basic case + if old == self: + return new + # Is it in expr? + if self.expr.has(old): + newexpr = self.expr._subs(old, new, **hints) + try: + return self._new_from_self(expr=newexpr) + except ValueError: + # Expr replacement leads to non-differentiable expression + # e.g `f.dx.subs(f: 1) = 1.dx = 0` + # returning zero + return sympy.S.Zero + + # In case `x0` was passed as a substitution instead of `(x0=` + if str(old) == 'x0': + return self._new_from_self(x0={self.dims[0]: new}) + + # Trying to substitute by another derivative with different metadata + # Only need to check if is a Derivative since one for the cases above would + # have found it + if isinstance(old, Derivative): + return self + + # Fall back if we didn't catch any special case + return self.xreplace({old: new}, **hints) def _xreplace(self, subs): """ This is a helper method used internally by SymPy. We exploit it to postpone substitutions until evaluation. """ + # Return if no subs + if not subs: + return self, False + # Check if trying to replace the whole expression if self in subs: new = subs.pop(self) @@ -245,6 +261,7 @@ def _xreplace(self, subs): return new._xreplace(subs) except AttributeError: return new, True + subs = self._ppsubs + (subs,) # Postponed substitutions return self._new_from_self(subs=subs), True diff --git a/devito/operations/solve.py b/devito/operations/solve.py index 12b19697c7..0203dbe26d 100644 --- a/devito/operations/solve.py +++ b/devito/operations/solve.py @@ -73,9 +73,9 @@ def linsolve(expr, target, **kwargs): The symbol w.r.t. which the equation is rearranged. May be a `Function` or any other symbolic object. """ - c = factorize_target(expr, target) + c, expr = factorize_target(expr, target) if c != 0: - return -expr.xreplace({target: 0})/c + return -expr/c raise SolveError("No linear solution found") @@ -102,7 +102,7 @@ def _(expr): @singledispatch def factorize_target(expr, target): - return 1 if expr == target else 0 + return (1, 0) if expr == target else (0, expr) @factorize_target.register(Add) @@ -110,21 +110,31 @@ def factorize_target(expr, target): def _(expr, target): c = 0 if not expr.has(target): - return c + return c, expr + args = [] for a in expr.args: - c += factorize_target(a, target) + c1, a1 = factorize_target(a, target) + c += c1 + args.append(a1) - return c + return c, expr.func(*args, evaluate=False) @factorize_target.register(Mul) def _(expr, target): if not expr.has(target): - return 0 + return 0, expr c = 1 + args = [] for a in expr.args: - c *= a if not a.has(target) else factorize_target(a, target) - - return c + if not a.has(target): + c *= a + args.append(a) + else: + c1, a1 = factorize_target(a, target) + c *= c1 + args.append(a1) + + return c, expr.func(*args, evaluate=False) diff --git a/devito/symbolics/queries.py b/devito/symbolics/queries.py index 1d8d6aff1c..ec86ae7809 100644 --- a/devito/symbolics/queries.py +++ b/devito/symbolics/queries.py @@ -43,6 +43,11 @@ def q_function(expr): return isinstance(expr, DiscreteFunction) +def q_derivative(expr): + from devito.finite_differences.derivative import Derivative + return isinstance(expr, Derivative) + + def q_terminal(expr): return (expr.is_Symbol or expr.is_Indexed or diff --git a/devito/symbolics/search.py b/devito/symbolics/search.py index feb4d622aa..94e533a692 100644 --- a/devito/symbolics/search.py +++ b/devito/symbolics/search.py @@ -1,9 +1,10 @@ from devito.symbolics.queries import (q_indexed, q_function, q_terminal, q_leaf, - q_symbol, q_dimension) + q_symbol, q_dimension, q_derivative) from devito.tools import as_tuple __all__ = ['retrieve_indexed', 'retrieve_functions', 'retrieve_function_carriers', - 'retrieve_terminals', 'retrieve_symbols', 'retrieve_dimensions', 'search'] + 'retrieve_terminals', 'retrieve_symbols', 'retrieve_dimensions', + 'retrieve_derivatives', 'search'] class Search(object): @@ -170,3 +171,8 @@ def retrieve_terminals(exprs, mode='all', deep=False): def retrieve_dimensions(exprs, mode='all', deep=False): """Shorthand to retrieve the dimensions in ``exprs``.""" return search(exprs, q_dimension, mode, 'dfs', deep) + + +def retrieve_derivatives(exprs, mode='all', deep=False): + """Shorthand to retrieve the Derivatives in ``exprs``.""" + return search(exprs, q_derivative, mode, 'dfs', deep) diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index a303e929ae..d8c2ca68b7 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -521,8 +521,12 @@ def test_substitution(self): grid = Grid((11, 11)) f = Function(name="f", grid=grid, space_order=4) expr = f.dx + f + 1 + assert simplify(expr.subs(f.dx, 1) - (f + 2)) == 0 - assert simplify(expr.subs(f, -1) - f.dx) == 0 + # f.dx.subs(f, -1) = 0 and f.subs(f, -1) = -1 so + # expr.subs(f, -1) = 0 + assert simplify(expr.subs(f, -1)) == 0 + # f.dx -> 1, f -> -1 assert simplify(expr.subs({f.dx: 1, f: -1}) - 1) == 0 expr2 = expr.subs({'x0': 2}) diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index cb2602507b..046a0b01d3 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -10,7 +10,8 @@ from devito.ir import Expression, FindNodes from devito.symbolics import (retrieve_functions, retrieve_indexed, evalrel, # noqa CallFromPointer, Cast, DefFunction, FieldFromPointer, - INT, FieldFromComposite, IntDiv, ccode, uxreplace) + INT, FieldFromComposite, IntDiv, ccode, uxreplace, + retrieve_derivatives) from devito.tools import as_tuple from devito.types import (Array, Bundle, FIndexed, LocalObject, Object, Symbol as dSymbol) @@ -360,9 +361,9 @@ def test_solve_time(): eq = m * u.dt2 + u.dx sol = solve(eq, u.forward) # Check u.dx is not evaluated. Need to simplify because the solution - # contains some Dummy in the Derivatibe subs that make equality break. - # TODO: replace by retreive_derivatives after Fabio's PR - assert sympy.simplify(u.dx - sol.args[2].args[0].args[1]) == 0 + # contains some Dummy in the Derivative subs that make equality break. + assert len(retrieve_derivatives(sol)) == 1 + assert sympy.simplify(u.dx - retrieve_derivatives(sol)[0]) == 0 assert sympy.simplify(sol - (-dt**2*u.dx/m + 2.0*u - u.backward)) == 0 From ff0a65e5022f2cf84e7ff951f47584d60165014d Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Sat, 19 Aug 2023 15:01:40 +0000 Subject: [PATCH 05/12] compiler: Patch DiscreteFunction reconstruction --- devito/builtins/arithmetic.py | 4 +++- devito/types/dense.py | 11 +++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/devito/builtins/arithmetic.py b/devito/builtins/arithmetic.py index f5f8002cb1..48df46e5ae 100644 --- a/devito/builtins/arithmetic.py +++ b/devito/builtins/arithmetic.py @@ -83,7 +83,9 @@ def sum(f, dims=None): elif f.is_SparseTimeFunction: if f.time_dim in dims: # Sum over time -> SparseFunction - new_coords = f.coordinates._rebuild(name="%ssum_coords" % f.name) + new_coords = f.coordinates._rebuild( + name="%ssum_coords" % f.name, initializer=f.coordinates.initializer + ) out = dv.SparseFunction(name="%ssum" % f.name, grid=f.grid, dimensions=new_dims, npoint=f.shape[1], coordinates=new_coords) diff --git a/devito/types/dense.py b/devito/types/dense.py index 6bb5d1a34f..b7ca53a0c0 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -57,7 +57,7 @@ class DiscreteFunction(AbstractFunction, ArgProvider, Differentiable): The type of the underlying data object. """ - __rkwargs__ = AbstractFunction.__rkwargs__ + ('staggered', 'initializer') + __rkwargs__ = AbstractFunction.__rkwargs__ + ('staggered',) def __init_finalize__(self, *args, function=None, **kwargs): # Staggering metadata @@ -638,7 +638,7 @@ def space_dimensions(self): @property def initializer(self): - if self._data is not None: + if isinstance(self._data, np.ndarray): return self.data_with_halo.view(np.ndarray) else: return self._initializer @@ -868,6 +868,13 @@ def _arg_finalize(self, args, alias=None): key = alias or self return {key.name: self._C_make_dataobj(args[key.name])} + # Pickling support + + @property + def _pickle_rkwargs(self): + # Picklying carries data over, if available + return tuple(self.__rkwargs__) + ('initializer',) + class Function(DiscreteFunction): From e981c407b50e3f5f06853e8bf3c6945d48ef4b5c Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 21 Aug 2023 07:49:26 +0000 Subject: [PATCH 06/12] compiler: Patch Weights __eq__ and hashing --- devito/finite_differences/differentiable.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index 25e258de81..1d09488e95 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -613,15 +613,15 @@ def __init_finalize__(self, *args, **kwargs): def __eq__(self, other): return (isinstance(other, Weights) and - self.dimension is other.dimension and self.name == other.name and + self.dimension == other.dimension and self.indices == other.indices and self.weights == other.weights) __hash__ = sympy.Basic.__hash__ def _hashable_content(self): - return super()._hashable_content() + (self.name,) + tuple(self.weights) + return (self.name, self.dimension, hash(tuple(self.weights))) @property def dimension(self): From f1d1567f0a7d78709d85844421b54cff2fb5585d Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 21 Aug 2023 10:57:19 +0000 Subject: [PATCH 07/12] compiler: Fix Differentiable._subs --- devito/finite_differences/differentiable.py | 4 ++-- tests/test_derivatives.py | 7 +++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index 1d09488e95..1ca7e29506 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -121,9 +121,9 @@ def _eval_at(self, func): return self.func(*[getattr(a, '_eval_at', lambda x: a)(func) for a in self.args]) def _subs(self, old, new, **hints): - if old is self: + if old == self: return new - if old is new: + if old == new: return self args = list(self.args) for i, arg in enumerate(args): diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index d8c2ca68b7..123975d63d 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -539,6 +539,13 @@ def test_substitution(self): expr3 = expr.subs({f.dx: f.dx2, 'x0': 2}) assert simplify(expr3 - (f.dx2(x0=2) + f + 1)) == 0 + # Test substitution with reconstructed objects + x, y = grid.dimensions + f1 = f.func(x + 1, y) + f2 = f.func(x + 1, y) + assert f1 is not f2 + assert f1.subs(f2, -1) == -1 + class TestTwoStageEvaluation(object): From 1a61319f1afc1176e61e4ff75107e1bcba11e030 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 21 Aug 2023 11:01:46 +0000 Subject: [PATCH 08/12] tests: Update expected output --- tests/test_gpu_common.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_gpu_common.py b/tests/test_gpu_common.py index 9e7af8c5fa..a93d280fc7 100644 --- a/tests/test_gpu_common.py +++ b/tests/test_gpu_common.py @@ -220,7 +220,7 @@ def test_tasking_unfused_two_locks(self): # Check generated code assert len(retrieve_iteration_tree(op)) == 3 - assert len([i for i in FindSymbols().visit(op) if isinstance(i, Lock)]) == 1 + 2 + assert len([i for i in FindSymbols().visit(op) if isinstance(i, Lock)]) == 2 sections = FindNodes(Section).visit(op) assert len(sections) == 4 assert (str(sections[1].body[0].body[0].body[0].body[0]) == @@ -412,7 +412,7 @@ def test_streaming_basic(self, opt, ntmps): assert np.all(u.data[1] == 36) @pytest.mark.parametrize('opt,ntmps,nfuncs', [ - (('buffering', 'streaming', 'orchestrate'), 9, 5), + (('buffering', 'streaming', 'orchestrate'), 8, 5), (('buffering', 'streaming', 'fuse', 'orchestrate', {'fuse-tasks': True}), 6, 5), ]) def test_streaming_two_buffers(self, opt, ntmps, nfuncs): @@ -724,7 +724,7 @@ def test_composite_buffering_tasking_multi_output(self): assert len(retrieve_iteration_tree(op1)) == 8 assert len(retrieve_iteration_tree(op2)) == 5 symbols = FindSymbols().visit(op1) - assert len([i for i in symbols if isinstance(i, Lock)]) == 1 + 2 + assert len([i for i in symbols if isinstance(i, Lock)]) == 2 threads = [i for i in symbols if isinstance(i, PThreadArray)] assert len(threads) == 2 assert threads[0].size.size == async_degree From ed6031ea3d250db8aa492551542e30b71a44a23e Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 21 Aug 2023 14:25:00 +0000 Subject: [PATCH 09/12] compiler: Patch reconstruction of MPI-distributed SparseFunctions --- devito/types/dense.py | 4 ++++ devito/types/sparse.py | 28 +++++++++++++++++++++------- examples/seismic/source.py | 2 +- 3 files changed, 26 insertions(+), 8 deletions(-) diff --git a/devito/types/dense.py b/devito/types/dense.py index b7ca53a0c0..19f6aa4769 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -287,6 +287,10 @@ def shape_global(self): ----- In an MPI context, this is the *global* domain region shape, which is therefore identical on all MPI ranks. + + Issues + ------ + * https://github.com/devitocodes/devito/issues/1498 """ if self.grid is None: return self.shape diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 3b5678a665..7ce74e0586 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -40,11 +40,11 @@ class AbstractSparseFunction(DiscreteFunction): _sub_functions = () """SubFunctions encapsulated within this AbstractSparseFunction.""" - __rkwargs__ = DiscreteFunction.__rkwargs__ + ('npoint', 'space_order') + __rkwargs__ = DiscreteFunction.__rkwargs__ + ('npoint_global', 'space_order') def __init_finalize__(self, *args, **kwargs): super(AbstractSparseFunction, self).__init_finalize__(*args, **kwargs) - self._npoint = kwargs['npoint'] + self._npoint = kwargs.get('npoint', kwargs.get('npoint_global')) self._space_order = kwargs.get('space_order', 0) # Dynamically add derivative short-cuts @@ -76,7 +76,7 @@ def __shape_setup__(cls, **kwargs): if grid is None: raise TypeError('Need `grid` argument') shape = kwargs.get('shape') - npoint = kwargs['npoint'] + npoint = kwargs.get('npoint', kwargs.get('npoint_global')) if shape is None: glb_npoint = SparseDistributor.decompose(npoint, grid.distributor) shape = (glb_npoint[grid.distributor.myrank],) @@ -90,6 +90,17 @@ def _halo_exchange(self): def npoint(self): return self.shape[self._sparse_position] + @property + def npoint_global(self): + """ + Global `npoint`s. This only differs from `self.npoint` in an MPI context. + + Issues + ------ + * https://github.com/devitocodes/devito/issues/1498 + """ + return self._npoint + @property def space_order(self): """The space order.""" @@ -490,8 +501,11 @@ def __distributor_setup__(self, **kwargs): A `SparseDistributor` handles the SparseFunction decomposition based on physical ownership, and allows to convert between global and local indices. """ - return SparseDistributor(kwargs['npoint'], self._sparse_dim, - kwargs['grid'].distributor) + return SparseDistributor( + kwargs.get('npoint', kwargs.get('npoint_global')), + self._sparse_dim, + kwargs['grid'].distributor + ) @property def coordinates(self): @@ -831,8 +845,8 @@ class SparseTimeFunction(AbstractSparseTimeFunction, SparseFunction): is_SparseTimeFunction = True - __rkwargs__ = (AbstractSparseTimeFunction.__rkwargs__ + - SparseFunction.__rkwargs__) + __rkwargs__ = tuple(filter_ordered(AbstractSparseTimeFunction.__rkwargs__ + + SparseFunction.__rkwargs__)) def interpolate(self, expr, offset=0, u_t=None, p_t=None, increment=False): """ diff --git a/examples/seismic/source.py b/examples/seismic/source.py index 692746e648..cdaf1ef4ae 100644 --- a/examples/seismic/source.py +++ b/examples/seismic/source.py @@ -110,7 +110,7 @@ def __args_setup__(cls, *args, **kwargs): kwargs['nt'] = kwargs['time_range'].num # Either `npoint` or `coordinates` must be provided - npoint = kwargs.get('npoint') + npoint = kwargs.get('npoint', kwargs.get('npoint_global')) if npoint is None: coordinates = kwargs.get('coordinates', kwargs.get('coordinates_data')) if coordinates is None: From e0d79c3b3171b612793b5c234a4985fb258ed35e Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 21 Aug 2023 20:53:05 -0400 Subject: [PATCH 10/12] ci: fix solve test fircing expansion --- tests/test_symbolics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py index 046a0b01d3..27d40c93c2 100644 --- a/tests/test_symbolics.py +++ b/tests/test_symbolics.py @@ -364,7 +364,7 @@ def test_solve_time(): # contains some Dummy in the Derivative subs that make equality break. assert len(retrieve_derivatives(sol)) == 1 assert sympy.simplify(u.dx - retrieve_derivatives(sol)[0]) == 0 - assert sympy.simplify(sol - (-dt**2*u.dx/m + 2.0*u - u.backward)) == 0 + assert sympy.simplify(sympy.expand(sol - (-dt**2*u.dx/m + 2.0*u - u.backward))) == 0 @pytest.mark.parametrize('expr,subs,expected', [ From 9b63f9131f69e95dba6848c085960dc75d1015da Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 21 Aug 2023 22:06:40 -0400 Subject: [PATCH 11/12] api: add _latex print for function for notebooks --- devito/types/basic.py | 2 ++ examples/seismic/tutorials/07_DRP_schemes.ipynb | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/devito/types/basic.py b/devito/types/basic.py index 10912259e6..51dbba181e 100644 --- a/devito/types/basic.py +++ b/devito/types/basic.py @@ -858,6 +858,8 @@ def __str__(self): def _sympystr(self, printer): return str(self) + _latex = _sympystr + def __eq__(self, other): try: return (self.function is other.function and diff --git a/examples/seismic/tutorials/07_DRP_schemes.ipynb b/examples/seismic/tutorials/07_DRP_schemes.ipynb index 80734fb4d6..3e6b16e601 100644 --- a/examples/seismic/tutorials/07_DRP_schemes.ipynb +++ b/examples/seismic/tutorials/07_DRP_schemes.ipynb @@ -99,7 +99,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "W(x, 1, u(t, x, y), x)*u(t, x, y) + W(x - h_x, 1, u(t, x, y), x)*u(t, x - h_x, y) + W(x + h_x, 1, u(t, x, y), x)*u(t, x + h_x, y)\n" + "u(t, x, y)*W(x, 1, u(t, x, y), x) + u(t, x - h_x, y)*W(x - h_x, 1, u(t, x, y), x) + u(t, x + h_x, y)*W(x + h_x, 1, u(t, x, y), x)\n" ] } ], From 25c856ad944e1253f0cf959b3560e18cead6ca1a Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 21 Aug 2023 22:12:38 -0400 Subject: [PATCH 12/12] api: fix symbolic coefficient conservation at rebuild --- devito/mpatches/.__init__.py.swp | Bin 12288 -> 0 bytes devito/mpatches/.rationaltools.py.swp | Bin 12288 -> 0 bytes devito/types/basic.py | 6 +++++- devito/types/dense.py | 2 +- 4 files changed, 6 insertions(+), 2 deletions(-) delete mode 100644 devito/mpatches/.__init__.py.swp delete mode 100644 devito/mpatches/.rationaltools.py.swp diff --git a/devito/mpatches/.__init__.py.swp b/devito/mpatches/.__init__.py.swp deleted file mode 100644 index 48f848b3e8678b2eef767d23fc7e5bf397b8c404..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12288 zcmeI%K}*9h6u|K}@1p1zhz^7a>zWDTSMcC*9xd70VW4el8wOtVgYCL+;m-R8yW4i| zScIK6_&*4b@RCRJyM;`{*XI}Y6pW>gOOe%P9ls6kWP2epPxYg23r(>*k}4T4bQtRH zC2i$v*Kiak3+vrJWrgvx_oNy{kz^+5I_I{1Q0*~fm0RmMNhu#!=L?Kpx2R4^?sZi0R#|0009ILKmY**5I_I{1Vl_M zvaJ95f9lcv|L@=b?HnKk5I_I{1Q0*~0R#|0009ILICg<~=`t0R##@)0)H|0}%4UTt Ny}DPbuX6WcNmVs>m1OB*{vEPB;fM0;qz`H;SBtRGVX**-5fTO^Tz?auB_A2lS@F;LS z@Ws`Py$QSlJPF(ce03FLF9NH;-N0U8CvXSw&y|dQ2Yd^h2Hpam2d2OTxEELgJm4bW zn=2Um8h96Y4bZ?*;0Q1PZUJrve!ZNr&w$d1SsGh-~iAC z{=O8k0v`ZpfLDO$fED0Y;0JL0J@7H`5pW849e53R9ykb49-g_Fv4??&fSZB8>wlcN z0PE|2?6NMCtn`eUNHdi=-^c4)W=FOd_A(Wh+FSS1iY~XFw6(-pYB&10xr)u7BC;x+ zhQH@8>0Zat)h9}_I2xaWGB*esai?+<%~xuzO(#Tm(fi6NQBKJBXcoA9M=`) zhCEfqa@bwyraIG6s*Om4s?z%WfY4FkWI{RKK&%_sv13P-iO#GCuRAE`U3}2B_q+G@ zX5FQkb#njy{qsK887z_L9*mOklAB7>^`h#`8fyMa$$AwF1!TY~DZ8YV#0kxR4Q34E8m6BQ%4+R=5C%;~T?oJ7;`f6Z}zlEQJ00<0RZ z7zM61^I4yleSTtYb|F3nvO-;TWq8sH zXmP|tdY;J~T|40-p2$%yqQI(|s^fA^=6#qHW-+f0YDnrdS1EQJa5GYI6z~M602S&! z*UIY9iVue(6S?Ws^IW4C4!NGHA`a_bg96Pte4C1NB14ZKQbt-G2TFuE1{1X8`lpWJ`5y-GxF+etGuFTgY7n&zUFs1996jzMeB-B}K66|qe@n2{wdq4mH diff --git a/devito/types/basic.py b/devito/types/basic.py index 51dbba181e..d7a422b39e 100644 --- a/devito/types/basic.py +++ b/devito/types/basic.py @@ -6,6 +6,7 @@ import numpy as np import sympy + from sympy.core.assumptions import _assume_rules from cached_property import cached_property @@ -855,11 +856,14 @@ def __str__(self): __repr__ = __str__ - def _sympystr(self, printer): + def _sympystr(self, printer, **kwargs): return str(self) _latex = _sympystr + def _pretty(self, printer, **kwargs): + return printer._print_Function(self, func_name=self.name) + def __eq__(self, other): try: return (self.function is other.function and diff --git a/devito/types/dense.py b/devito/types/dense.py index 19f6aa4769..aec69bc1eb 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -57,7 +57,7 @@ class DiscreteFunction(AbstractFunction, ArgProvider, Differentiable): The type of the underlying data object. """ - __rkwargs__ = AbstractFunction.__rkwargs__ + ('staggered',) + __rkwargs__ = AbstractFunction.__rkwargs__ + ('staggered', 'coefficients') def __init_finalize__(self, *args, function=None, **kwargs): # Staggering metadata