diff --git a/devito/ir/clusters/cluster.py b/devito/ir/clusters/cluster.py index de4cd8f29f2..0dc3200b4f8 100644 --- a/devito/ir/clusters/cluster.py +++ b/devito/ir/clusters/cluster.py @@ -334,7 +334,6 @@ def dspace(self): # Construct the `intervals` of the DataSpace, that is a global, # Dimension-centric view of the data space - intervals = IntervalGroup.generate('union', *parts.values()) # E.g., `db0 -> time`, but `xi NOT-> x` intervals = intervals.promote(lambda d: not d.is_Sub) diff --git a/devito/ir/equations/algorithms.py b/devito/ir/equations/algorithms.py index adc462d0592..2606d2dd2dc 100644 --- a/devito/ir/equations/algorithms.py +++ b/devito/ir/equations/algorithms.py @@ -56,9 +56,7 @@ def handle_indexed(indexed): # such as A[3]) indexeds = retrieve_indexed(expr, deep=True) for i in indexeds: - expl_dims = {d for (d, e) in zip(i.function.dimensions, i.indices) - if e.is_integer} - extra.update(expl_dims) + extra.update({d for d in i.function.dimensions if i.indices[d].is_integer}) # Enforce determinism extra = filter_sorted(extra) @@ -75,8 +73,9 @@ def handle_indexed(indexed): # of `x`, besides `(x, xi)`, we also have to add `(time, x)` so that we # obtain the desired ordering `(time, x, xi)`. W/o `(time, x)`, the ordering # `(x, time, xi)` might be returned instead, which would be non-sense - implicit_relations.update({tuple(filter_ordered(d.root for d in i)) - for i in relations}) + for i in relations: + dims = [di for d in i for di in (d.index, d)] + implicit_relations.update({tuple(filter_ordered(dims))}) ordering = PartialOrderTuple(extra, relations=(relations | implicit_relations)) diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py index 2f7a47c0e95..98ba9da51a8 100644 --- a/devito/ir/support/basic.py +++ b/devito/ir/support/basic.py @@ -660,8 +660,8 @@ def is_const(self, dim): """ True if a constant dependence, that is no Dimensions involved, False otherwise. """ - return (self.source.aindices.get(dim, None) is None and - self.sink.aindices.get(dim, None) is None and + return (self.source.aindices.get(dim) is None and + self.sink.aindices.get(dim) is None and self.distance_mapper.get(dim, 0) == 0) @memoized_meth diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index 4550b5b143e..b480c9f1139 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -141,7 +141,7 @@ def _weights(self): raise NotImplementedError @property - def _gdim(self): + def _gdims(self): return self.grid.dimensions @property @@ -153,12 +153,15 @@ def _rdim(self): parent = self.sfunction.dimensions[-1] dims = [CustomDimension("r%s%s" % (self.sfunction.name, d.name), -self.r+1, self.r, 2*self.r, parent) - for d in self._gdim] + for d in self._gdims] - return DimensionTuple(*dims, getters=self._gdim) + return DimensionTuple(*dims, getters=self._gdims) def _augment_implicit_dims(self, implicit_dims): - return as_tuple(implicit_dims) + self.sfunction.dimensions + if self.sfunction._sparse_position == -1: + return self.sfunction.dimensions + as_tuple(implicit_dims) + else: + return as_tuple(implicit_dims) + self.sfunction.dimensions def _coeff_temps(self, implicit_dims): return [] @@ -178,7 +181,7 @@ def _interp_idx(self, variables, implicit_dims=None): # Coefficient symbol expression temps.extend(self._coeff_temps(implicit_dims)) - for ((di, d), rd, p) in zip(enumerate(self._gdim), self._rdim, pos): + for ((di, d), rd, p) in zip(enumerate(self._gdims), self._rdim, pos): # Add conditional to avoid OOB lb = sympy.And(rd + p >= d.symbolic_min - self.r, evaluate=False) ub = sympy.And(rd + p <= d.symbolic_max + self.r, evaluate=False) @@ -330,7 +333,7 @@ class LinearInterpolator(WeightedInterpolator): @property def _weights(self): c = [(1 - p) * (1 - r) + p * r - for (p, d, r) in zip(self._point_symbols, self._gdim, self._rdim)] + for (p, d, r) in zip(self._point_symbols, self._gdims, self._rdim)] return Mul(*c) @cached_property @@ -345,7 +348,7 @@ def _coeff_temps(self, implicit_dims): pmap = self.sfunction._position_map poseq = [Eq(self._point_symbols[d], pos - floor(pos), implicit_dims=implicit_dims) - for (d, pos) in zip(self._gdim, pmap.keys())] + for (d, pos) in zip(self._gdims, pmap.keys())] return poseq diff --git a/devito/passes/clusters/aliases.py b/devito/passes/clusters/aliases.py index 099c3a929eb..cfe9645b1b6 100644 --- a/devito/passes/clusters/aliases.py +++ b/devito/passes/clusters/aliases.py @@ -837,7 +837,7 @@ def lower_schedule(schedule, meta, sregistry, ftemps): # This prevents cases such as `floor(a*b)` with `a` and `b` floats # that would creat a temporary `int r = b` leading to erronous numerical results # Such cases happen with the positions for sparse functions for example. - dtype = sympy_dtype(pivot, meta.dtype) or meta.dtype + dtype = sympy_dtype(pivot, meta.dtype) if writeto: # The Dimensions defining the shape of Array diff --git a/devito/passes/iet/mpi.py b/devito/passes/iet/mpi.py index 1343a33b8a5..de87bb13a82 100644 --- a/devito/passes/iet/mpi.py +++ b/devito/passes/iet/mpi.py @@ -47,9 +47,11 @@ def _drop_halospots(iet): # If a HaloSpot is outside any iteration it is not needed for iters, halo_spots in MapNodes(Iteration, HaloSpot, 'groupby').visit(iet).items(): + if iters: + continue for hs in halo_spots: for f, v in hs.fmapper.items(): - if not iters and v.loc_indices: + if v.loc_indices: mapper[hs].add(f) # Transform the IET introducing the "reduced" HaloSpots diff --git a/devito/passes/iet/parpragma.py b/devito/passes/iet/parpragma.py index c0b6f2d1554..016726c0017 100644 --- a/devito/passes/iet/parpragma.py +++ b/devito/passes/iet/parpragma.py @@ -426,7 +426,12 @@ def _make_guard(self, parregion): def _make_nested_partree(self, partree): # Apply heuristic - if self.nhyperthreads <= self.nested or partree.root.is_ParallelAtomic: + if self.nhyperthreads <= self.nested: + return partree + + # Loop nest with atomic reductions are more likely to have less latency + # keep outer loop parallel + if partree.root.is_ParallelAtomic: return partree # Note: there might be multiple sub-trees amenable to nested parallelism, diff --git a/devito/symbolics/inspection.py b/devito/symbolics/inspection.py index 235a708fa4e..c1659c8cf67 100644 --- a/devito/symbolics/inspection.py +++ b/devito/symbolics/inspection.py @@ -281,4 +281,5 @@ def sympy_dtype(expr, default): return default else: # Infer expression dtype from its arguments - return infer_dtype([sympy_dtype(a, default) for a in expr.args]) + dtype = infer_dtype([sympy_dtype(a, default) for a in expr.args]) + return dtype or default diff --git a/devito/tools/data_structures.py b/devito/tools/data_structures.py index be1bd4edc91..539f75d5934 100644 --- a/devito/tools/data_structures.py +++ b/devito/tools/data_structures.py @@ -66,7 +66,7 @@ def __getnewargs_ex__(self): # objects with varying number of attributes return (tuple(self), dict(self.__dict__)) - def get(self, key, val): + def get(self, key, val=None): return self._getters.get(key, val) @@ -605,6 +605,7 @@ class UnboundTuple(object): """ A simple data structure that returns the last element forever once reached """ + def __init__(self, items): self.items = as_tuple(items) self.last = len(self.items) diff --git a/devito/tools/dtypes_lowering.py b/devito/tools/dtypes_lowering.py index c8fe8a3fa55..0b3cd53ebfb 100644 --- a/devito/tools/dtypes_lowering.py +++ b/devito/tools/dtypes_lowering.py @@ -130,6 +130,9 @@ def dtype_to_mpitype(dtype): def dtype_to_mpidtype(dtype): + """ + Map numpy type to MPI internal types for communication + """ from devito.mpi import MPI return MPI._typedict[np.dtype(dtype).char] diff --git a/devito/types/dense.py b/devito/types/dense.py index a6ac169f3cd..94fa7d8b778 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -1460,9 +1460,7 @@ def parent(self): @property def origin(self): - """ - SubFunction have zero origin - """ + # SubFunction have zero origin return DimensionTuple(*(0 for _ in range(self.ndim)), getters=self.dimensions) diff --git a/devito/types/dimension.py b/devito/types/dimension.py index b65cb76309d..20772441f1a 100644 --- a/devito/types/dimension.py +++ b/devito/types/dimension.py @@ -107,8 +107,6 @@ class Dimension(ArgProvider): is_Incr = False is_Block = False - indirect = False - # Prioritize self's __add__ and __sub__ to construct AffineIndexAccessFunction _op_priority = sympy.Expr._op_priority + 1. @@ -183,6 +181,14 @@ def min_name(self): def max_name(self): return "%s_M" % self.name + @property + def indirect(self): + return False + + @property + def index(self): + return self if self.indirect is True else getattr(self, 'parent', self) + @property def is_const(self): return False @@ -456,7 +462,6 @@ class DerivedDimension(BasicDimension): """ is_Derived = True - indirect = False __rargs__ = Dimension.__rargs__ + ('parent',) __rkwargs__ = () @@ -819,10 +824,6 @@ def condition(self): def indirect(self): return self._indirect - @property - def index(self): - return self if self.indirect is True else self.parent - @cached_property def free_symbols(self): retval = set(super().free_symbols) @@ -1216,7 +1217,7 @@ def __init_finalize__(self, name, symbolic_min=None, symbolic_max=None, self._symbolic_min = symbolic_min self._symbolic_max = symbolic_max self._symbolic_size = symbolic_size - self._parent = parent + self._parent = parent or BOTTOM super().__init_finalize__(name) @property diff --git a/devito/types/sparse.py b/devito/types/sparse.py index a2f0fd63c4e..162d7d08976 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -136,19 +136,6 @@ def __subfunc_setup__(self, key, suffix, dtype=None): return sf - @property - def npoint(self): - return self.shape[self._sparse_position] - - @property - def space_order(self): - """The space order.""" - return self._space_order - - @property - def r(self): - return self._radius - @property def _sparse_dim(self): return self.dimensions[self._sparse_position] @@ -179,6 +166,37 @@ def _coords_indices(self): np.floor((self.coordinates_data - self.grid.origin) / self.grid.spacing) ).astype(int) + @property + def _support(self): + """ + The grid points surrounding each sparse point within the radius of self's + injection/interpolation operators. + """ + max_shape = np.array(self.grid.shape).reshape(1, self.grid.dim) + minmax = lambda arr: np.minimum(max_shape, np.maximum(0, arr)) + return np.stack([minmax(self._coords_indices + s) for s in self._point_support], + axis=2) + + @property + def _dist_datamap(self): + """ + Mapper ``M : MPI rank -> required sparse data``. + """ + return self.grid.distributor.glb_to_rank(self._support) or {} + + @property + def npoint(self): + return self.shape[self._sparse_position] + + @property + def space_order(self): + """The space order.""" + return self._space_order + + @property + def r(self): + return self._radius + @property def gridpoints(self): try: @@ -207,28 +225,6 @@ def coordinates_data(self): except AttributeError: return None - @property - def _support(self): - """ - The grid points surrounding each sparse point within the radius of self's - injection/interpolation operators. - """ - max_shape = np.array(self.grid.shape).reshape(1, self.grid.dim) - minmax = lambda arr: np.minimum(max_shape, np.maximum(0, arr)) - return np.stack([minmax(self._coords_indices + s) for s in self._point_support], - axis=2) - - @property - def _dist_datamap(self): - """ - Mapper ``M : MPI rank -> required sparse data``. - """ - return self.grid.distributor.glb_to_rank(self._support) or {} - - @cached_property - def dist_origin(self): - return self._dist_origin - @cached_property def _pos_symbols(self): return [Symbol(name='pos%s' % d, dtype=np.int32) @@ -265,6 +261,10 @@ def _dist_reorder_mask(self): if d is not self._sparse_dim) return ret + @cached_property + def dist_origin(self): + return self._dist_origin + def interpolate(self, *args, **kwargs): """ Implement an interpolation operation from the grid onto the given sparse points diff --git a/tests/test_dimension.py b/tests/test_dimension.py index b9fb7ea2f52..7f867a88e47 100644 --- a/tests/test_dimension.py +++ b/tests/test_dimension.py @@ -1450,7 +1450,6 @@ def test_sparse_time_function(self): # Note the endpoint of the range is 12 because we inject at p.forward for i in range(1, 12): assert p.data[i].sum() == i - 1 - print(p.data[i, 10, 10, 10]) assert p.data[i, 10, 10, 10] == i - 1 for i in range(12, 20): assert np.all(p.data[i] == 0) diff --git a/tests/test_pickle.py b/tests/test_pickle.py index 58842b75800..f95deb6048e 100644 --- a/tests/test_pickle.py +++ b/tests/test_pickle.py @@ -135,6 +135,24 @@ def test_precomputed_sparse_function(self, mode, pickle): assert sf.dtype == new_sf.dtype assert sf.npoint == new_sf.npoint == 3 + def test_alias_sparse_function(self, pickle): + grid = Grid(shape=(3,)) + sf = SparseFunction(name='sf', grid=grid, npoint=3, space_order=2, + coordinates=[(0.,), (1.,), (2.,)]) + sf.data[0] = 1. + + # Create alias + f0 = sf._rebuild(name='f0', alias=True) + pkl_f0 = pickle.dumps(f0) + new_f0 = pickle.loads(pkl_f0) + + assert f0.data is None and new_f0.data is None + assert f0.coordinates.data is None and new_f0.coordinates.data is None + + assert sf.space_order == f0.space_order == new_f0.space_order + assert sf.dtype == f0.dtype == new_f0.dtype + assert sf.npoint == f0.npoint == new_f0.npoint + def test_internal_symbols(self, pickle): s = dSymbol(name='s', dtype=np.float32) pkl_s = pickle.dumps(s)