diff --git a/devito/__init__.py b/devito/__init__.py index be66ba07f5..2eed5c9847 100644 --- a/devito/__init__.py +++ b/devito/__init__.py @@ -1,14 +1,12 @@ from __future__ import absolute_import from devito.base import * # noqa -from devito.finite_difference import * # noqa from devito.dimension import * # noqa +from devito.equation import * # noqa +from devito.finite_difference import * # noqa from devito.grid import * # noqa -from devito.function import Forward, Backward # noqa from devito.logger import error, warning, info, set_log_level # noqa from devito.parameters import * # noqa -from devito.region import * # noqa -from devito.symbolics import * # noqa from devito.tools import * # noqa from devito.compiler import compiler_registry diff --git a/devito/dimension.py b/devito/dimension.py index 6d055e3361..7e8e193472 100644 --- a/devito/dimension.py +++ b/devito/dimension.py @@ -27,13 +27,11 @@ class Dimension(AbstractSymbol): potential iteration space. :param name: Name of the dimension symbol. - :param reverse: Optional, Traverse dimension in reverse order (default False) :param spacing: Optional, symbol for the spacing along this dimension. """ def __new__(cls, name, **kwargs): newobj = sympy.Symbol.__new__(cls, name) - newobj._reverse = kwargs.get('reverse', False) newobj._spacing = kwargs.get('spacing', Scalar(name='h_%s' % name)) return newobj @@ -86,27 +84,16 @@ def start_name(self): def end_name(self): return "%s_e" % self.name - @property - def reverse(self): - return self._reverse - @property def spacing(self): return self._spacing - @reverse.setter - def reverse(self, val): - # TODO: this is an outrageous hack. TimeFunctions are updating this value - # at construction time. - self._reverse = val - @property def base(self): return self def _hashable_content(self): - return super(Dimension, self)._hashable_content() +\ - (self.reverse, self.spacing) + return super(Dimension, self)._hashable_content() + (self.spacing,) def argument_defaults(self, size=None): """ @@ -148,7 +135,6 @@ class SpaceDimension(Dimension): symbols. :param name: Name of the dimension symbol. - :param reverse: Traverse dimension in reverse order (default False) :param spacing: Optional, symbol for the spacing along this dimension. """ @@ -163,7 +149,6 @@ class TimeDimension(Dimension): time dimensions should inherit from :class:`TimeDimension`. :param name: Name of the dimension symbol. - :param reverse: Traverse dimension in reverse order (default False) :param spacing: Optional, symbol for the spacing along this dimension. """ @@ -193,10 +178,6 @@ def __new__(cls, name, parent, **kwargs): def parent(self): return self._parent - @property - def reverse(self): - return self.parent.reverse - @property def spacing(self): return self.parent.spacing diff --git a/devito/dle/blocking_utils.py b/devito/dle/blocking_utils.py index 853ed36c50..b9dc3afabd 100644 --- a/devito/dle/blocking_utils.py +++ b/devito/dle/blocking_utils.py @@ -201,11 +201,9 @@ class IterationFold(Iteration): is_IterationFold = True - def __init__(self, nodes, dimension, limits, index=None, offsets=None, - properties=None, pragmas=None, uindices=None, folds=None): - super(IterationFold, self).__init__(nodes, dimension, limits, index, offsets, - properties, uindices, pragmas) - self.folds = folds + def __init__(self, *args, **kwargs): + self.folds = kwargs.pop('folds', None) + super(IterationFold, self).__init__(*args, **kwargs) def __repr__(self): properties = "" diff --git a/devito/dse/backends/advanced.py b/devito/dse/backends/advanced.py index ea434ae399..fafecd4bb3 100644 --- a/devito/dse/backends/advanced.py +++ b/devito/dse/backends/advanced.py @@ -2,7 +2,7 @@ from collections import OrderedDict -from devito.ir import Cluster, ClusterGroup, groupby +from devito.ir import Cluster, ClusterGroup, IterationSpace, groupby from devito.dse.aliases import collect from devito.dse.backends import BasicRewriter, dse_pass from devito.symbolics import Eq, estimate_cost, xreplace_constrained, iq_timeinvariant @@ -140,9 +140,12 @@ def _eliminate_inter_stencil_redundancies(self, cluster, template, **kwargs): function = make(c) # Build new Cluster expression = Eq(Indexed(function, *indices), origin) - ispace = cluster.ispace.subtract(alias.anti_stencil.boxify().negate()) + intervals, sub_iterators, directions = cluster.ispace.args + # Adjust intervals + intervals = intervals.subtract(alias.anti_stencil.boxify().negate()) if all(time_invariants[i] for i in alias.aliased): - ispace = ispace.drop([i.dim for i in ispace.intervals if i.dim.is_Time]) + intervals = intervals.drop([i for i in intervals.dimensions if i.is_Time]) + ispace = IterationSpace(intervals, sub_iterators, directions) alias_clusters.append(Cluster([expression], ispace)) # Update substitution rules for aliased, distance in alias.with_distance: diff --git a/devito/equation.py b/devito/equation.py new file mode 100644 index 0000000000..b03fd132f6 --- /dev/null +++ b/devito/equation.py @@ -0,0 +1,62 @@ +"""User API to specify equations.""" + +import sympy + +__all__ = ['Eq', 'Inc', 'DOMAIN', 'INTERIOR'] + + +class Eq(sympy.Eq): + + """ + A :class:`sympy.Eq` that accepts the additional keyword parameter ``region``. + + The ``region``, an object of type :class:`Region`, may be used to restrict + the execution of the equation to a sub-domain. + """ + + is_Increment = False + + def __new__(cls, *args, **kwargs): + kwargs['evaluate'] = False + region = kwargs.pop('region', DOMAIN) + obj = sympy.Eq.__new__(cls, *args, **kwargs) + obj._region = region + return obj + + +class Inc(Eq): + + """ + A :class:`Eq` performing a linear increment. + """ + + is_Increment = True + + +class Region(object): + + """ + A region of the computational domain over which a :class:`Function` is + discretized. + """ + + def __init__(self, name): + self._name = name + + def __repr__(self): + return self._name + + def __eq__(self, other): + return isinstance(other, Region) and self._name == other._name + + +DOMAIN = Region('DOMAIN') +""" +Represent the physical domain of the PDE; that is, domain = boundary + interior +""" + +INTERIOR = Region('INTERIOR') +""" +Represent the physical interior domain of the PDE; that is, PDE boundaries are +not included. +""" diff --git a/devito/function.py b/devito/function.py index 4b59063924..fe470bcd94 100644 --- a/devito/function.py +++ b/devito/function.py @@ -6,44 +6,23 @@ import numpy as np from psutil import virtual_memory -from devito.parameters import configuration -from devito.logger import debug, error, warning -from devito.data import Data, first_touch +from devito.arguments import ArgumentMap from devito.cgen_utils import INT, FLOAT +from devito.data import Data, first_touch from devito.dimension import Dimension -from devito.types import SymbolicFunction, AbstractCachedSymbol -from devito.arguments import ArgumentMap +from devito.equation import Eq, Inc from devito.finite_difference import (centered, cross_derivative, first_derivative, left, right, second_derivative, generic_derivative, second_cross_derivative) -from devito.symbolics import Eq, Inc, indexify, retrieve_indexed +from devito.logger import debug, error, warning +from devito.parameters import configuration +from devito.symbolics import indexify, retrieve_indexed +from devito.types import SymbolicFunction, AbstractCachedSymbol from devito.tools import EnrichedTuple __all__ = ['Constant', 'Function', 'TimeFunction', 'SparseFunction', - 'SparseTimeFunction', 'Forward', 'Backward'] - - -class TimeAxis(object): - """Direction in which to advance the time index on - :class:`TimeFunction` objects. - - :param axis: Either 'Forward' or 'Backward' - """ - - def __init__(self, axis): - assert axis in ['Forward', 'Backward'] - self._axis = {'Forward': 1, 'Backward': -1}[axis] - - def __eq__(self, other): - return self._axis == other._axis - - def __repr__(self): - return {-1: 'Backward', 1: 'Forward'}[self._axis] - - -Forward = TimeAxis('Forward') -Backward = TimeAxis('Backward') + 'SparseTimeFunction'] class Constant(AbstractCachedSymbol): @@ -333,6 +312,7 @@ def argument_defaults(self, alias=None): # Collect default dimension arguments from all indices for i, s, o in zip(self.indices, self.shape, self.staggered): args.update(i.argument_defaults(size=s+o)) + return args def argument_values(self, alias=None, **kwargs): @@ -348,6 +328,9 @@ def argument_values(self, alias=None, **kwargs): # Add value override for own data if it is provided if self.name in kwargs: new = kwargs.pop(self.name) + if len(new.shape) != self.ndim: + raise ValueError("Array shape %s does not match" % (new.shape, ) + + "dimensions %s" % (self.indices, )) if isinstance(new, TensorFunction): # Set new values and re-derive defaults values[key] = new._data_buffer @@ -355,11 +338,10 @@ def argument_values(self, alias=None, **kwargs): else: # We've been provided a pure-data replacement (array) values[key] = new - # TODO: Re-derive defaults from shape of array + # Add value overrides for all associated dimensions + for i, s, o in zip(self.indices, new.shape, self.staggered): + values.update(i.argument_defaults(size=s+o)) - # Add value overrides for all associated dimensions - for i in self.indices: - values.update(i.argument_values(**kwargs)) return values @@ -659,6 +641,7 @@ def __init__(self, *args, **kwargs): self.time_dim = kwargs.get('time_dim') self.time_order = kwargs.get('time_order', 1) + self.save = type(kwargs.get('save', None) or None) if not isinstance(self.time_order, int): raise ValueError("'time_order' must be int") @@ -744,6 +727,23 @@ def dt2(self): return self.diff(_t, _t).as_finite_difference(indt) + def argument_values(self, alias=None, **kwargs): + """ + Returns a map of argument values after evaluating user input. + + :param kwargs: Dictionary of user-provided argument overrides. + :param alias: (Optional) name under which to store values. + """ + # Check if data has the right dimension + if self.name in kwargs: + new = kwargs.get(self.name) + if isinstance(new, TimeFunction) and new.save != self.save: + raise TypeError("Incorrect value encountered, save should be %s" % + self.save) + + values = super(TimeFunction, self).argument_values(alias=alias, **kwargs) + return values + class SparseFunction(TensorFunction): """ @@ -816,7 +816,7 @@ def __indices_setup__(cls, **kwargs): if dimensions is not None: return dimensions else: - return (Dimension(name='p'),) + return (Dimension(name='p_%s' % kwargs["name"]),) @classmethod def __shape_setup__(cls, **kwargs): @@ -1022,16 +1022,20 @@ def argument_values(self, alias=None, **kwargs): :param alias: (Optional) name under which to store values. """ # Take a copy of the replacement before super pops it from kwargs + new = kwargs.get(self.name) - values = super(SparseFunction, self).argument_values(alias=alias, **kwargs) + key = alias or self.name + + values = super(SparseFunction, self).argument_values(alias=key, **kwargs) + if new is not None and isinstance(new, SparseFunction): # If we've been replaced with a SparseFunction, # we need to re-derive defaults and values... - values.update(new.argument_defaults(alias=alias).reduce_all()) + values.update(new.argument_defaults(alias=key).reduce_all()) values.update(new.coordinates.argument_defaults(alias=self.coordinates.name)) else: # ..., but if not, we simply need to recurse over children. - values.update(self.coordinates.argument_values(alias=alias, **kwargs)) + values.update(self.coordinates.argument_values(alias=key, **kwargs)) return values @@ -1080,7 +1084,7 @@ def __indices_setup__(cls, **kwargs): if dimensions is not None: return dimensions else: - return (kwargs['grid'].time_dim, Dimension(name='p')) + return (kwargs['grid'].time_dim, Dimension(name='p_%s' % kwargs["name"])) @classmethod def __shape_setup__(cls, **kwargs): diff --git a/devito/ir/clusters/algorithms.py b/devito/ir/clusters/algorithms.py index 1e6ca3ecd8..34578f80a3 100644 --- a/devito/ir/clusters/algorithms.py +++ b/devito/ir/clusters/algorithms.py @@ -1,10 +1,11 @@ from collections import OrderedDict -from devito.ir.support import Scope +from devito.ir.support import (Scope, IterationSpace, detect_flow_directions, + force_directions, group_expressions) from devito.ir.clusters.cluster import PartialCluster, ClusterGroup from devito.symbolics import CondEq, CondNe, IntDiv, xreplace_indices from devito.types import Scalar -from devito.tools import Bunch, flatten, powerset +from devito.tools import flatten, powerset __all__ = ['clusterize', 'groupby'] @@ -34,7 +35,7 @@ def groupby(clusters): anti = scope.d_anti.carried() - scope.d_anti.increment flow = scope.d_flow - (scope.d_flow.inplace() + scope.d_flow.increment) funcs = [i.function for i in anti] - if candidate.ispace == c.ispace and\ + if candidate.ispace.is_compatible(c.ispace) and\ all(is_local(i, candidate, c, clusters) for i in funcs): # /c/ will be fused into /candidate/. All fusion-induced anti # dependences are eliminated through so called "index bumping and @@ -221,33 +222,35 @@ def bump_and_contract(targets, source, sink): def clusterize(exprs): """Group a sequence of :class:`ir.Eq`s into one or more :class:`Cluster`s.""" - # Build a graph capturing the dependencies among the input tensor expressions - mapper = OrderedDict() - for i, e1 in enumerate(exprs): - trace = [e2 for e2 in exprs[:i] if Scope([e2, e1]).has_dep] + [e1] - trace.extend([e2 for e2 in exprs[i+1:] if Scope([e1, e2]).has_dep]) - mapper[e1] = Bunch(trace=trace, ispace=e1.ispace) + # Group expressions based on data dependences + groups = group_expressions(exprs) - # Derive the iteration spaces - queue = list(mapper) - while queue: - target = queue.pop(0) - - ispaces = [mapper[i].ispace for i in mapper[target].trace] - - coerced_ispace = mapper[target].ispace.intersection(*ispaces) - - if coerced_ispace != mapper[target].ispace: - # Something has changed, need to propagate the update - mapper[target].ispace = coerced_ispace - queue.extend([i for i in mapper[target].trace if i not in queue]) - - # Build a PartialCluster for each tensor expression clusters = ClusterGroup() - for k, v in mapper.items(): - if k.is_Tensor: - scalars = [i for i in v.trace[:v.trace.index(k)] if i.is_Scalar] - clusters.append(PartialCluster(scalars + [k], v.ispace)) + for g in groups: + # Coerce iteration space of each expression in each group + mapper = OrderedDict([(e, e.ispace) for e in g]) + flowmap = detect_flow_directions(g) + queue = list(g) + while queue: + v = queue.pop(0) + + intervals, sub_iterators, directions = mapper[v].args + forced, clashes = force_directions(flowmap, lambda i: directions.get(i)) + for e in g: + intervals = intervals.intersection(mapper[e].intervals.drop(clashes)) + directions = {i: forced[i] for i in directions} + coerced_ispace = IterationSpace(intervals, sub_iterators, directions) + + # Need update propagation ? + if coerced_ispace != mapper[v]: + mapper[v] = coerced_ispace + queue.extend([i for i in g if i not in queue]) + + # Wrap each tensor expression in a PartialCluster + for k, v in mapper.items(): + if k.is_Tensor: + scalars = [i for i in g[:g.index(k)] if i.is_Scalar] + clusters.append(PartialCluster(scalars + [k], v)) # Group PartialClusters together where possible clusters = groupby(clusters) diff --git a/devito/ir/clusters/cluster.py b/devito/ir/clusters/cluster.py index f517d656b0..8872bc57ef 100644 --- a/devito/ir/clusters/cluster.py +++ b/devito/ir/clusters/cluster.py @@ -1,9 +1,8 @@ from cached_property import cached_property from frozendict import frozendict -from devito.ir.support import IterationSpace +from devito.ir.equations import ClusterizedEq from devito.ir.clusters.graph import FlowGraph -from devito.symbolics import Eq __all__ = ["Cluster", "ClusterGroup"] @@ -25,7 +24,7 @@ class PartialCluster(object): """ def __init__(self, exprs, ispace, atomics=None, guards=None): - self._exprs = list(exprs) + self._exprs = list(ClusterizedEq(i, ispace) for i in exprs) self._ispace = ispace self._atomics = set(atomics or []) self._guards = guards or {} @@ -74,7 +73,7 @@ def squash(self, other): """Concatenate the expressions in ``other`` to those in ``self``. ``self`` and ``other`` must have same ``ispace``. Duplicate expressions are dropped.""" - assert self.ispace == other.ispace + assert self.ispace.is_compatible(other.ispace) self.exprs.extend([i for i in other.exprs if i not in self.exprs]) @@ -83,7 +82,10 @@ class Cluster(PartialCluster): """A Cluster is an immutable :class:`PartialCluster`.""" def __init__(self, exprs, ispace, atomics=None, guards=None): - self._exprs = tuple(Eq(*i.args, evaluate=False) for i in exprs) + # Keep expressions ordered based on information flow + self._exprs = exprs + self._exprs = tuple(ClusterizedEq(v, ispace) for v in self.trace.values()) + self._ispace = ispace self._atomics = frozenset(atomics or ()) self._guards = frozendict(guards or {}) @@ -119,15 +121,6 @@ class ClusterGroup(list): """An iterable of :class:`PartialCluster`s.""" - @property - def ispace(self): - """ - Return the union of all Clusters' iteration spaces. - """ - if not self: - return IterationSpace([]) - return IterationSpace.generate('intersection', *[i.ispace for i in self]) - def unfreeze(self): """ Return a new ClusterGroup in which all of ``self``'s Clusters have diff --git a/devito/ir/equations/equation.py b/devito/ir/equations/equation.py index 4b34a366ea..01568dd43c 100644 --- a/devito/ir/equations/equation.py +++ b/devito/ir/equations/equation.py @@ -1,16 +1,14 @@ -from collections import OrderedDict - from sympy import Eq from devito.dimension import SubDimension -from devito.ir.support import Interval, DataSpace, IterationSpace, Stencil -from devito.region import DOMAIN, INTERIOR -from devito.symbolics import dimension_sort, indexify +from devito.equation import DOMAIN, INTERIOR +from devito.ir.support import IterationSpace, Any, compute_intervals, compute_directions +from devito.symbolics import FrozenExpr, dimension_sort, indexify -__all__ = ['LoweredEq'] +__all__ = ['LoweredEq', 'ClusterizedEq', 'IREq'] -class EqMixin(object): +class IREq(object): """ A mixin providing operations common to all :mod:`ir` equation types. @@ -24,21 +22,21 @@ def is_Scalar(self): def is_Tensor(self): return self.lhs.is_Indexed + @property + def dimensions(self): + return self.ispace.dimensions + -class LoweredEq(Eq, EqMixin): +class LoweredEq(Eq, IREq): """ LoweredEq(expr, subs=None) - A SymPy equation with associated iteration and data spaces. + A SymPy equation with an associated iteration space. All :class:`Function` objects within ``expr`` get indexified and thus turned into objects of type :class:`types.Indexed`. - The data space is an object of type :class:`DataSpace`. It represents the - data points accessed by the equation along each :class:`Dimension`. The - :class:`Dimension`s are extracted directly from the equation. - The iteration space is an object of type :class:`IterationSpace`. It represents the iteration points along each :class:`Dimension` that the equation may traverse with the guarantee that no out-of-bounds accesses @@ -57,11 +55,10 @@ def __new__(cls, *args, **kwargs): expr = Eq.__new__(cls, *args, evaluate=False) assert isinstance(stamp, Eq) expr.is_Increment = stamp.is_Increment - expr.dspace = stamp.dspace expr.ispace = stamp.ispace return expr else: - raise ValueError("Cannot construct Eq from args=%s " + raise ValueError("Cannot construct LoweredEq from args=%s " "and kwargs=%s" % (str(args), str(kwargs))) # Indexification @@ -83,27 +80,16 @@ def __new__(cls, *args, **kwargs): expr = expr.xreplace(mapper) ordering = [mapper.get(i, i) for i in ordering] - # Get the accessed data points - stencil = Stencil(expr) - - # Split actual Intervals (the data spaces) from the "derived" iterators, - # to build an IterationSpace - iterators = OrderedDict() - for i in ordering: - if i.is_NonlinearDerived: - iterators.setdefault(i.parent, []).append(stencil.entry(i)) - else: - iterators.setdefault(i, []) - intervals = [] - for k, v in iterators.items(): - offs = set.union(set(stencil.get(k)), *[i.ofs for i in v]) - intervals.append(Interval(k, min(offs), max(offs))) + # Compute iteration space + intervals, iterators = compute_intervals(expr) + intervals = sorted(intervals, key=lambda i: ordering.index(i.dim)) + directions, _ = compute_directions(expr, lambda i: Any) + ispace = IterationSpace([i.negate() for i in intervals], iterators, directions) # Finally create the LoweredEq with all metadata attached expr = super(LoweredEq, cls).__new__(cls, expr.lhs, expr.rhs, evaluate=False) expr.is_Increment = getattr(input_expr, 'is_Increment', False) - expr.dspace = DataSpace(intervals) - expr.ispace = IterationSpace([i.negate() for i in intervals], iterators) + expr.ispace = ispace return expr @@ -112,3 +98,38 @@ def xreplace(self, rules): def func(self, *args): return super(LoweredEq, self).func(*args, stamp=self, evaluate=False) + + +class ClusterizedEq(Eq, IREq, FrozenExpr): + + """ + ClusterizedEq(expr, ispace) + + A SymPy equation carrying its own :class:`IterationSpace`. Suitable for + use in a :class:`Cluster`. + + Unlike a :class:`LoweredEq`, a ClusterizedEq is "frozen", meaning that any + call to ``xreplace`` will not trigger re-evaluation (e.g., mathematical + simplification) of the expression. + """ + + def __new__(cls, *args, **kwargs): + # Parse input + if len(args) == 2: + maybe_ispace = args[1] + if isinstance(maybe_ispace, IterationSpace): + input_expr = args[0] + expr = Eq.__new__(cls, *input_expr.args, evaluate=False) + expr.is_Increment = input_expr.is_Increment + expr.ispace = maybe_ispace + else: + expr = Eq.__new__(cls, *args, evaluate=False) + expr.ispace = kwargs['ispace'] + expr.is_Increment = kwargs.get('is_Increment', False) + else: + raise ValueError("Cannot construct ClusterizedEq from args=%s " + "and kwargs=%s" % (str(args), str(kwargs))) + return expr + + def func(self, *args, **kwargs): + return super(ClusterizedEq, self).func(*args, evaluate=False, ispace=self.ispace) diff --git a/devito/ir/iet/nodes.py b/devito/ir/iet/nodes.py index fae965485a..5110709de8 100644 --- a/devito/ir/iet/nodes.py +++ b/devito/ir/iet/nodes.py @@ -14,8 +14,8 @@ from devito.ir.iet import (IterationProperty, SEQUENTIAL, PARALLEL, VECTOR, ELEMENTAL, REMAINDER, WRAPPABLE, tagger, ntags) +from devito.ir.support import Forward from devito.dimension import Dimension -from devito.ir.support import Stencil from devito.symbolics import as_symbol, retrieve_terminals from devito.tools import as_tuple, filter_ordered, filter_sorted, flatten import devito.types as types @@ -277,11 +277,6 @@ def shape(self): """ return () if self.is_scalar else self.expr.lhs.shape - @property - def stencil(self): - """Compute the stencil of the expression.""" - return Stencil(self.expr) - @property def free_symbols(self): """Return all :class:`Symbol` objects used by this :class:`Expression`.""" @@ -297,6 +292,8 @@ class Iteration(Node): tuple of the form (start, finish, stepping). :param index: Symbol to be used as iteration variable. :param offsets: A 2-tuple of start and end offsets to honour in the loop. + :param direction: The :class:`IterationDirection` of the Iteration. Defaults + to ``Forward``. :param properties: A bag of :class:`IterationProperty` objects, decorating the Iteration (sequential, parallel, vectorizable, ...). :param pragmas: A bag of pragmas attached to this Iteration. @@ -310,18 +307,13 @@ class Iteration(Node): _traversable = ['nodes'] def __init__(self, nodes, dimension, limits, index=None, offsets=None, - properties=None, pragmas=None, uindices=None): + direction=None, properties=None, pragmas=None, uindices=None): # Ensure we deal with a list of Expression objects internally - nodes = as_tuple(nodes) - self.nodes = as_tuple([n if isinstance(n, Node) else Expression(n) - for n in nodes]) - assert all(isinstance(i, Node) for i in self.nodes) + self.nodes = as_tuple(nodes) self.dim = dimension self.index = index or self.dim.name - # Store direction, as it might change on the dimension - # before we use it during code generation. - self.reverse = self.dim.reverse + self.direction = direction or Forward # Generate loop limits if isinstance(limits, Iterable): diff --git a/devito/ir/iet/scheduler.py b/devito/ir/iet/scheduler.py index 1d34a64f9a..9aa3eb6869 100644 --- a/devito/ir/iet/scheduler.py +++ b/devito/ir/iet/scheduler.py @@ -62,8 +62,8 @@ def iet_make(clusters, dtype): needed = intervals[index:] # Build Expressions - body = [Expression(v, np.int32 if cluster.trace.is_index(k) else dtype) - for k, v in cluster.trace.items()] + body = [Expression(e, np.int32 if cluster.trace.is_index(e.lhs) else dtype) + for e in cluster.exprs] if not needed: body = List(body=body) @@ -80,18 +80,21 @@ def iet_make(clusters, dtype): value = (i.dim + o) % modulo uindices.append(UnboundedIndex(vname, value, value, j, j + o)) + # Retrieve the iteration direction + direction = cluster.ispace.directions[i.dim] + # Update IET and scheduling if i.dim in cluster.guards: # Must wrap within an if-then scope body = Conditional(cluster.guards[i.dim], body) - iteration = Iteration(body, i.dim, i.dim.limits, - offsets=i.limits, uindices=uindices) + iteration = Iteration(body, i.dim, i.dim.limits, offsets=i.limits, + direction=direction, uindices=uindices) # Adding (None, None) ensures that nested iterations won't # be reused by the next cluster scheduling.extend([(None, None), (i, iteration)]) else: - iteration = Iteration(body, i.dim, i.dim.limits, - offsets=i.limits, uindices=uindices) + iteration = Iteration(body, i.dim, i.dim.limits, offsets=i.limits, + direction=direction, uindices=uindices) scheduling.append((i, iteration)) # Prepare for next dimension diff --git a/devito/ir/iet/visitors.py b/devito/ir/iet/visitors.py index 119aa97bfd..24d14ede0f 100644 --- a/devito/ir/iet/visitors.py +++ b/devito/ir/iet/visitors.py @@ -14,6 +14,7 @@ from devito.cgen_utils import blankline, ccode from devito.exceptions import VisitorException from devito.ir.iet.nodes import Node +from devito.ir.support.space import Backward from devito.tools import as_tuple, filter_sorted, flatten, ctypes_to_C, GenericVisitor @@ -249,8 +250,8 @@ def visit_Iteration(self, o): else: end = o.limits[1] - # For reverse dimensions flip loop bounds - if o.reverse: + # For backward direction flip loop bounds + if o.direction == Backward: loop_init = 'int %s = %s' % (o.index, ccode('%s - 1' % end)) loop_cond = '%s >= %s' % (o.index, ccode(start)) loop_inc = '%s -= %s' % (o.index, o.limits[2]) diff --git a/devito/ir/support/__init__.py b/devito/ir/support/__init__.py index 5fcdaa6afb..abf1681432 100644 --- a/devito/ir/support/__init__.py +++ b/devito/ir/support/__init__.py @@ -1,3 +1,4 @@ from devito.ir.support.basic import * # noqa from devito.ir.support.space import * # noqa from devito.ir.support.stencil import * # noqa +from devito.ir.support.utils import * # noqa diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py index b159a3caa6..76cb7e3afb 100644 --- a/devito/ir/support/basic.py +++ b/devito/ir/support/basic.py @@ -1,12 +1,12 @@ from cached_property import cached_property - -from sympy import Basic, Eq +from sympy import Basic, S from devito.dimension import Dimension +from devito.ir.support.space import Any, Backward from devito.symbolics import retrieve_terminals, q_affine, q_inc from devito.tools import as_tuple, is_integer, filter_sorted -__all__ = ['Scope'] +__all__ = ['Vector', 'IterationInstance', 'Access', 'TimedAccess', 'Scope'] class Vector(tuple): @@ -161,16 +161,21 @@ class IterationInstance(Vector): """ A representation of the iteration and data points accessed by an - :class:`Indexed` object. Two different concepts are distinguished: - - * The index functions; that is, the expressions telling what *iteration* - space point is accessed. - * The ``findices``: that is, the :class:`Dimension`s telling what *data* - space point is accessed. + :class:`Indexed` object. Three different concepts are distinguished: + + * Index functions: the expressions telling what *iteration* space point + is accessed. + * ``aindices``: the :class:`Dimension`s acting as iteration variables. + There is one aindex for each index function. If the index function + is non-affine, then it may not be possible to detect its aindex; + in such a case, None is used as placeholder. + * ``findices``: the :class:`Dimension`s telling what *data* space point + is accessed. """ def __new__(cls, indexed): obj = super(IterationInstance, cls).__new__(cls, *indexed.indices) + # findices obj.findices = tuple(indexed.base.function.indices) if len(obj.findices) != len(set(obj.findices)): raise ValueError("Illegal non-unique `findices`") @@ -209,6 +214,44 @@ def __getitem__(self, index): raise TypeError("IterationInstance indices must be integers, slices, or " "Dimensions, not %s" % type(index)) + @cached_property + def _cached_findices_index(self): + # Avoiding to call self.findices.index repeatedly speeds analysis up + return {fi: i for i, fi in enumerate(self.findices)} + + @cached_property + def index_mode(self): + index_mode = [] + for i, fi in zip(self, self.findices): + if is_integer(i): + index_mode.append('regular') + elif q_affine(i, fi): + index_mode.append('regular') + else: + index_mode.append('irregular') + return tuple(index_mode) + + @cached_property + def aindices(self): + aindices = [] + for i, fi in zip(self, self.findices): + if is_integer(i): + aindices.append(None) + elif q_affine(i, fi): + aindices.append(fi) + else: + dims = {i for i in i.free_symbols if isinstance(i, Dimension)} + aindices.append(dims.pop() if len(dims) == 1 else None) + return tuple(aindices) + + @property + def is_regular(self): + return all(i == 'regular' for i in self.index_mode) + + @property + def is_irregular(self): + return not self.is_regular + def distance(self, other, findex=None): """Compute vector distance from ``self`` to ``other``. If ``findex`` is supplied, compute the vector distance up to and including ``findex``.""" @@ -218,8 +261,8 @@ def distance(self, other, findex=None): raise TypeError("Cannot compute distance due to mismatching `findices`") if findex is not None: try: - limit = self.findices.index(findex) + 1 - except ValueError: + limit = self._cached_findices_index[findex] + 1 + except KeyError: raise TypeError("Cannot compute distance as `findex` not in `findices`") else: limit = self.rank @@ -262,6 +305,10 @@ def __eq__(self, other): isinstance(other, Access) and\ self.function == other.function + @property + def name(self): + return self.function.name + @property def is_read(self): return self.mode in ['R', 'RI'] @@ -284,31 +331,7 @@ def is_increment(self): def __repr__(self): mode = '\033[1;37;31mW\033[0m' if self.is_write else '\033[1;37;32mR\033[0m' - return "%s<%s,[%s]>" % (mode, self.function.name, - ', '.join(str(i) for i in self)) - - -class IterationFunction(object): - - """A representation of a function describing the direction and the step - in which an iteration space is traversed.""" - - _KNOWN = [] - - def __init__(self, name): - assert isinstance(name, str) and name not in IterationFunction._KNOWN - self._name = name - IterationFunction._KNOWN.append(name) - - def __repr__(self): - return self._name - - -INC = IterationFunction('++') -"""Unit-stride increment.""" - -DEC = IterationFunction('--') -"""Unit-stride decrement""" + return "%s<%s,[%s]>" % (mode, self.name, ', '.join(str(i) for i in self)) class TimedAccess(Access): @@ -350,38 +373,25 @@ class TimedAccess(Access): """ - def __new__(cls, indexed, mode, timestamp): + def __new__(cls, indexed, mode, timestamp, directions): assert is_integer(timestamp) obj = super(TimedAccess, cls).__new__(cls, indexed, mode) obj.timestamp = timestamp - obj.direction = [DEC if i.reverse else INC for i in obj.findices] + obj.directions = [directions.get(i, Any) for i in obj.findices] return obj def __eq__(self, other): return super(TimedAccess, self).__eq__(other) and\ isinstance(other, TimedAccess) and\ - self.direction == other.direction + self.directions == other.directions def __lt__(self, other): if not isinstance(other, TimedAccess): raise TypeError("Cannot compare with object of type %s" % type(other)) - if self.direction != other.direction: + if self.directions != other.directions: raise TypeError("Cannot compare due to mismatching `direction`") return super(TimedAccess, self).__lt__(other) - @property - def index_mode(self): - return ['regular' if (is_integer(i) or q_affine(i, fi)) else 'irregular' - for i, fi in zip(self, self.findices)] - - @property - def is_regular(self): - return all(i == 'regular' for i in self.index_mode) - - @property - def is_irregular(self): - return not self.is_regular - def lex_eq(self, other): return self.timestamp == other.timestamp @@ -401,16 +411,23 @@ def lex_lt(self, other): return self.timestamp < other.timestamp def distance(self, other, findex=None): - if (self.direction != other.direction) or (self.rank != other.rank): - raise TypeError("Cannot order due to mismatching `direction` and/or `rank`") - ret = super(TimedAccess, self).distance(other, findex) - if findex is not None: - limit = self.findices.index(findex) + 1 - direction = self.direction[:limit] - else: - direction = self.direction - assert len(direction) == len(ret) - return Vector(*[i if d == INC else (-i) for i, d in zip(ret, direction)]) + if self.rank != other.rank: + raise TypeError("Cannot order due to mismatching `rank`") + if not self.rank: + return Vector() + findex = findex or self.findices[-1] + ret = [] + for i, sd, od in zip(self.findices, self.directions, other.directions): + if sd == od: + ret = list(super(TimedAccess, self).distance(other, i)) + if i == findex: + break + else: + ret.append(S.Infinity) + break + directions = self.directions[:self._cached_findices_index[i] + 1] + assert len(directions) == len(ret) + return Vector(*[(-i) if d == Backward else i for i, d in zip(ret, directions)]) class Dependence(object): @@ -536,12 +553,11 @@ class Scope(object): def __init__(self, exprs): """ - Initialize a Scope, which represents a group of :class:`Access` objects - extracted from some expressions ``exprs``. The expressions are to be - provided as they appear in program order. + A Scope represents a group of :class:`TimedAccess` objects extracted + from some :class:`IREq` ``exprs``. The expressions must be provided + in program order. """ exprs = as_tuple(exprs) - assert all(isinstance(i, Eq) for i in exprs) self.reads = {} self.writes = {} @@ -550,11 +566,11 @@ def __init__(self, exprs): for j in retrieve_terminals(e.rhs): v = self.reads.setdefault(j.base.function, []) mode = 'R' if not q_inc(e) else 'RI' - v.append(TimedAccess(j, mode, i)) + v.append(TimedAccess(j, mode, i, e.ispace.directions)) # write v = self.writes.setdefault(e.lhs.base.function, []) mode = 'W' if not q_inc(e) else 'WI' - v.append(TimedAccess(e.lhs, mode, i)) + v.append(TimedAccess(e.lhs, mode, i, e.ispace.directions)) def getreads(self, function): return as_tuple(self.reads.get(function)) diff --git a/devito/ir/support/space.py b/devito/ir/support/space.py index 0ffe0b2ae6..3840b3b768 100644 --- a/devito/ir/support/space.py +++ b/devito/ir/support/space.py @@ -1,9 +1,10 @@ import abc from collections import OrderedDict -from devito.tools import as_tuple +from devito.tools import as_tuple, filter_ordered -__all__ = ['NullInterval', 'Interval', 'DataSpace', 'IterationSpace'] +__all__ = ['NullInterval', 'Interval', 'IntervalGroup', 'IterationSpace', + 'Forward', 'Backward', 'Any'] class AbstractInterval(object): @@ -75,7 +76,7 @@ def union(self, o): if self.dim == o.dim: return o._rebuild() else: - return Space([self._rebuild(), o._rebuild()]) + return IntervalGroup([self._rebuild(), o._rebuild()]) def overlap(self, o): return False @@ -124,7 +125,7 @@ def union(self, o): elif o.is_Null and self.dim == o.dim: return self._rebuild() else: - return Space([self._rebuild(), o._rebuild()]) + return IntervalGroup([self._rebuild(), o._rebuild()]) def subtract(self, o): if self.dim != o.dim or o.is_Null: @@ -155,115 +156,164 @@ def __hash__(self): return hash((self.dim.name, self.lower, self.upper)) -class Space(object): +class IntervalGroup(tuple): """ - A bag of :class:`Interval`s. - - The intervals input ordering is honored. + A sequence of :class:`Interval`s with set-like operations exposed. """ - def __init__(self, intervals): - self.intervals = as_tuple(intervals) - - def __repr__(self): - return "%s[%s]" % (self.__class__.__name__, - ', '.join([repr(i) for i in self.intervals])) - def __eq__(self, o): - return set(self.intervals) == set(o.intervals) - - def _construct(self, intervals): - return Space(intervals) + return set(self) == set(o) - @property - def size(self): - return len(self.intervals) + def __repr__(self): + return "IntervalGroup[%s]" % (', '.join([repr(i) for i in self])) @property - def empty(self): - return self.size == 0 + def dimensions(self): + return filter_ordered([i.dim for i in self]) @property def is_well_defined(self): """ - Return True if the space :class:`Interval`s are over different - :class:`Dimension`s, False otherwise. + Return True if all :class:`Interval`s are over different :class:`Dimension`s, + False otherwise. """ - dims = [i.dim for i in self.intervals] - return len(dims) == len(set(dims)) + return len(self.dimensions) == len(set(self.dimensions)) @classmethod - def generate(self, op, *spaces): + def generate(self, op, *interval_groups): """ - generate(op, *spaces) + generate(op, *interval_groups) - Create a new :class:`Space` from the iterative application of the - operation ``op`` to the :class:`Space`s in ``spaces``. + Create a new :class:`IntervalGroup` from the iterative application of the + operation ``op`` to the :class:`IntervalGroup`s in ``interval_groups``. :param op: Any legal :class:`Interval` operation, such as ``intersection`` or ``union``. This should be provided as a string. - :param spaces: An iterable of :class:`Space`s. + :param interval_groups: An iterable of :class:`IntervalGroup`s. Example ------- - space0 = Space([Interval(x, 1, -1)]) - space1 = Space([Interval(x, 2, -2), Interval(y, 3, -3)]) - space2 = Space([Interval(y, 2, -2), Interval(z, 1, -1)]) + ig0 = IntervalGroup([Interval(x, 1, -1)]) + ig1 = IntervalGroup([Interval(x, 2, -2), Interval(y, 3, -3)]) + ig2 = IntervalGroup([Interval(y, 2, -2), Interval(z, 1, -1)]) - res = Space.generate('intersection', space0, space1, space2) - res --> Space([Interval(x, 2, -2), Interval(y, 3, -3), Interval(z, 1, -1)]) + ret = IntervalGroup.generate('intersection', ig0, ig1, ig2) + ret -> IntervalGroup([Interval(x, 2, -2), Interval(y, 3, -3), Interval(z, 1, -1)]) """ mapper = OrderedDict() - for i in spaces: - for interval in i.intervals: - mapper.setdefault(interval.dim, []).append(interval) - return Space([Interval._apply_op(v, op) for v in mapper.values()]) - - def intersection(self, *spaces): - mapper = OrderedDict([(i.dim, [i]) for i in self.intervals]) - for i in spaces: - for interval in i.intervals: - mapper.get(interval.dim, []).append(interval) - return self._construct([Interval._apply_op(v, 'intersection') - for v in mapper.values()]) + for ig in interval_groups: + for i in ig: + mapper.setdefault(i.dim, []).append(i) + return IntervalGroup(Interval._apply_op(v, op) for v in mapper.values()) + + def intersection(self, o): + mapper = OrderedDict([(i.dim, i) for i in o]) + intervals = [i.intersection(mapper.get(i.dim, i)) for i in self] + return IntervalGroup(intervals) def subtract(self, o): - mapper = OrderedDict([(i.dim, i) for i in o.intervals]) - intervals = [i.subtract(mapper.get(i.dim, NullInterval(i.dim))) - for i in self.intervals] - return self._construct(intervals) + mapper = OrderedDict([(i.dim, i) for i in o]) + intervals = [i.subtract(mapper.get(i.dim, NullInterval(i.dim))) for i in self] + return IntervalGroup(intervals) def drop(self, d): - return self._construct([i._rebuild() for i in self.intervals - if i.dim not in as_tuple(d)]) + return IntervalGroup([i._rebuild() for i in self if i.dim not in as_tuple(d)]) def negate(self): - return self._construct([i.negate() for i in self.intervals]) + return IntervalGroup([i.negate() for i in self]) - def __getitem__(self, dim): + def __getitem__(self, key): + if isinstance(key, (slice, int)): + return super(IntervalGroup, self).__getitem__(key) if not self.is_well_defined: raise ValueError("Cannot fetch Interval from ill defined Space") - for i in self.intervals: - if i.dim == dim: + for i in self: + if i.dim == key: return i -class DataSpace(Space): - pass +class IterationDirection(object): + + """ + A representation of the direction in which an iteration space is traversed. + """ + + def __init__(self, name): + self._name = name + + def __eq__(self, other): + return self._name == other._name + + def __repr__(self): + return self._name + + def __hash__(self): + return hash(self._name) -class IterationSpace(Space): +Forward = IterationDirection('++') +"""Forward iteration direction ('++').""" + +Backward = IterationDirection('--') +"""Backward iteration direction ('--').""" + +Any = IterationDirection('*') +"""Wildcard direction (both '++' and '--' would be OK).""" + + +class IterationSpace(object): """ - A :class:`Space` associating one or more :class:`Dimension`s to a - :class:`Interval`. The interval is the data space, whereas the dimensions - are the objects used to traverse the data space. + A representation of an iteration space and its traversal through + :class:`Interval`s and :class:`IterationDirection`s. + + :param intervals: An ordered sequence of :class:`Interval`s defining the + iteration space. + :param sub_iterators: (Optional) a mapper from :class:`Dimension`s in + ``intervals`` to iterables of :class:`DerivedDimension`, + which represent sub-dimensions along which the iteration + space is traversed. + :param directions: (Optional) a mapper from :class:`Dimension`s in + ``intervals`` to :class:`IterationDirection`s. """ - def __init__(self, intervals, sub_iterators): - super(IterationSpace, self).__init__(intervals) - self.sub_iterators = sub_iterators + def __init__(self, intervals, sub_iterators=None, directions=None): + self.intervals = IntervalGroup(as_tuple(intervals)) + self.sub_iterators = sub_iterators or {} + self.directions = directions or {} + + def __repr__(self): + ret = ', '.join(["%s%s" % (repr(i), repr(self.directions[i.dim])) + for i in self.intervals]) + return "IterationSpace[%s]" % ret - def _construct(self, intervals): - return IterationSpace(intervals, self.sub_iterators) + def __eq__(self, other): + return self.intervals == other.intervals and self.directions == other.directions + + def is_compatible(self, other): + """A relaxed version of ``__eq__``, in which only non-derived dimensions + are compared for equality.""" + return self.intervals == other.intervals and\ + self.nonderived_directions == other.nonderived_directions + + @property + def args(self): + return (self.intervals, self.sub_iterators, self.directions) + + @property + def size(self): + return len(self.intervals) + + @property + def empty(self): + return self.size == 0 + + @property + def dimensions(self): + sub_dims = [i.dim for v in self.sub_iterators.values() for i in v] + return filter_ordered(self.intervals.dimensions + sub_dims) + + @property + def nonderived_directions(self): + return {k: v for k, v in self.directions.items() if not k.is_Derived} diff --git a/devito/ir/support/stencil.py b/devito/ir/support/stencil.py index 26c9cd41f3..15cae88c4c 100644 --- a/devito/ir/support/stencil.py +++ b/devito/ir/support/stencil.py @@ -1,11 +1,7 @@ from collections import namedtuple -from sympy import Eq - -from devito.dimension import Dimension -from devito.ir.support.space import Interval, Space -from devito.symbolics import retrieve_indexed -from devito.tools import DefaultOrderedDict, flatten +from devito.ir.support.space import Interval, IntervalGroup +from devito.tools import DefaultOrderedDict __all__ = ['Stencil'] @@ -24,58 +20,25 @@ class Stencil(DefaultOrderedDict): Note: Expressions must have been indexified for a Stencil to be computed. """ - def __init__(self, *args): + def __init__(self, entries=None): """ Initialize the Stencil. - :param args: A Stencil may be created in several ways: :: - - * From a SymPy equation, or - * A list of elements of type: :: - * SymPy equation, or - * StencilEntry, or - * 2-tuple (Dimension, set) -- raw initialization + :param entries: An iterable of :class:`StencilEntry` or a 2-tuple + convertible into a :class:`StencilEntry` (i.e., a + :class:`Dimension` and a set). """ processed = [] - for i in args: - if isinstance(i, Eq): - processed.extend(self.extract(i).items()) + for i in (entries or []): + if isinstance(i, StencilEntry): + processed.append((i.dim, i.ofs)) + elif isinstance(i, tuple) and len(i) == 2: + entry = StencilEntry(*i) # Type checking + processed.append((entry.dim, entry.ofs)) else: - for j in i: - if isinstance(j, StencilEntry): - processed.append((j.dim, j.ofs)) - elif isinstance(j, tuple) and len(j) == 2: - entry = StencilEntry(*j) # Type checking - processed.append((entry.dim, entry.ofs)) - else: - raise RuntimeError('Cannot construct a Stencil for %s' % str(j)) + raise TypeError('Cannot construct a Stencil for %s' % str(i)) super(Stencil, self).__init__(set, processed) - @classmethod - def extract(cls, expr): - """ - Compute the stencil of ``expr``. - """ - indexeds = retrieve_indexed(expr, mode='all') - indexeds += flatten([retrieve_indexed(i) for i in e.indices] for e in indexeds) - - stencil = Stencil() - for e in indexeds: - for a in e.indices: - if isinstance(a, Dimension): - stencil[a].update([0]) - d = None - off = [0] - for i in a.args: - if isinstance(i, Dimension): - d = i - elif i.is_integer: - off += [int(i)] - if d is not None: - stencil[d].update(off) - - return stencil - @classmethod def union(cls, *dicts): """ @@ -154,9 +117,10 @@ def copy(self): def boxify(self): """ - Create a :class:`Space` from ``self``, with as many intervals as dimensions. + Create a :class:`IntervalGroup` from ``self``, with as many intervals + as dimensions in ``self``. """ - return Space([Interval(k, min(v), max(v)) for k, v in self.items()]) + return IntervalGroup([Interval(k, min(v), max(v)) for k, v in self.items()]) def __eq__(self, other): return self.entries == other.entries diff --git a/devito/ir/support/utils.py b/devito/ir/support/utils.py new file mode 100644 index 0000000000..b67726fce5 --- /dev/null +++ b/devito/ir/support/utils.py @@ -0,0 +1,193 @@ +from collections import OrderedDict, defaultdict +from itertools import groupby + +from devito.dimension import Dimension +from devito.ir.support.basic import Access, Scope +from devito.ir.support.space import Interval, Backward, Forward, Any +from devito.ir.support.stencil import Stencil +from devito.symbolics import retrieve_indexed +from devito.tools import as_tuple, flatten + +__all__ = ['compute_intervals', 'detect_flow_directions', 'compute_directions', + 'force_directions', 'group_expressions'] + + +def compute_intervals(expr): + """Return an iterable of :class:`Interval`s representing the data items + accessed by the :class:`sympy.Eq` ``expr``.""" + # Deep retrieval of indexed objects in /expr/ + indexeds = retrieve_indexed(expr, mode='all') + indexeds += flatten([retrieve_indexed(i) for i in e.indices] for e in indexeds) + + # Detect the indexeds' offsets along each dimension + stencil = Stencil() + for e in indexeds: + for a in e.indices: + if isinstance(a, Dimension): + stencil[a].update([0]) + d = None + off = [0] + for i in a.args: + if isinstance(i, Dimension): + d = i + elif i.is_integer: + off += [int(i)] + if d is not None: + stencil[d].update(off) + + # Determine intervals and their iterators + iterators = OrderedDict() + for i in stencil.dimensions: + if i.is_NonlinearDerived: + iterators.setdefault(i.parent, []).append(stencil.entry(i)) + else: + iterators.setdefault(i, []) + intervals = [] + for k, v in iterators.items(): + offs = set.union(set(stencil.get(k)), *[i.ofs for i in v]) + intervals.append(Interval(k, min(offs), max(offs))) + + return intervals, iterators + + +def detect_flow_directions(exprs): + """Return a mapper from :class:`Dimension`s to iterables of + :class:`IterationDirection`s representing the theoretically necessary + directions to evaluate ``exprs`` so that the information "naturally + flows" from an iteration to another.""" + exprs = as_tuple(exprs) + + writes = [Access(i.lhs, 'W') for i in exprs] + reads = flatten(retrieve_indexed(i.rhs, mode='all') for i in exprs) + reads = [Access(i, 'R') for i in reads] + + # Determine indexed-wise direction by looking at the vector distance + mapper = defaultdict(set) + for w in writes: + for r in reads: + if r.name != w.name: + continue + dimensions = [d for d in w.aindices if d is not None] + for d in dimensions: + try: + if w.distance(r, d) > 0: + mapper[d].add(Forward) + break + elif w.distance(r, d) < 0: + mapper[d].add(Backward) + break + else: + mapper[d].add(Any) + except TypeError: + # Nothing can be deduced + mapper[d].add(Any) + break + # Remainder + for d in dimensions[dimensions.index(d) + 1:]: + mapper[d].add(Any) + + # Add in any encountered Dimension + mapper.update({d: {Any} for d in flatten(i.aindices for i in reads + writes) + if d is not None and d not in mapper}) + + # Add in stepping dimensions (just in case they haven't been detected yet) + # note: stepping dimensions may force a direction on the parent + assert all(v == {Any} or mapper.get(k.parent, v) in [v, {Any}] + for k, v in mapper.items() if k.is_Stepping) + mapper.update({k.parent: set(v) for k, v in mapper.items() + if k.is_Stepping and mapper.get(k.parent) == {Any}}) + + # Add in derived dimensions parents + mapper.update({k.parent: set(v) for k, v in mapper.items() + if k.is_Derived and k.parent not in mapper}) + + return mapper + + +def compute_directions(exprs, key): + """ + Return a mapper ``M : D -> I`` where D is the set of :class:`Dimension`s + found in the input expressions ``exprs``, while I = {Any, Backward, + Forward} (i.e., the set of possible :class:`IterationDirection`s). + + The iteration direction is chosen so that the information "naturally flows" + from an iteration to another (i.e., to generate "flow" or "read-after-write" + dependencies). + + In the case of a clash (e.g., both Forward and Backward should be used + for a given dimension in order to have a flow dependence), the function + ``key : D -> I`` is used to pick one value. + """ + mapper = detect_flow_directions(exprs) + return force_directions(mapper, key) + + +def force_directions(mapper, key): + """ + Return a mapper ``M : D -> I`` where D is the set of :class:`Dimension`s + found in the input mapper ``M' : D -> {I}``, while I = {Any, Backward, + Forward} (i.e., the set of possible :class:`IterationDirection`s). + + The iteration direction is chosen so that the information "naturally flows" + from an iteration to another (i.e., to generate "flow" or "read-after-write" + dependencies). + + In the case of a clash (e.g., both Forward and Backward should be used + for a given dimension in order to have a flow dependence), the function + ``key : D -> I`` is used to pick one value. + """ + mapper = {k: set(v) for k, v in mapper.items()} + clashes = set(k for k, v in mapper.items() if len(v - {Any}) > 1) + directions = {} + for k, v in mapper.items(): + if len(v) == 1: + directions[k] = v.pop() + elif len(v) == 2: + try: + v.remove(Any) + directions[k] = v.pop() + except KeyError: + assert k in clashes + directions[k] = key(k) + else: + assert k in clashes + directions[k] = key(k) + + return directions, clashes + + +def group_expressions(exprs): + """``{exprs} -> ({exprs'}, {exprs''}, ...)`` where: :: + + * There are data dependences within exprs' and within exprs''; + * There are *no* data dependencies across exprs' and exprs''. + """ + # Partion based on data dependences + mapper = OrderedDict() + ngroups = 0 + for i, e1 in enumerate(exprs): + if e1 in mapper: + continue + found = False + for e2 in exprs[i+1:]: + if Scope([e1, e2]).has_dep: + v = mapper.get(e1, mapper.get(e2)) + if v is None: + ngroups += 1 + v = ngroups + mapper[e1] = mapper[e2] = v + found = True + if not found: + ngroups += 1 + mapper[e1] = ngroups + + # Reorder to match input ordering + groups = [] + data = sorted(mapper, key=lambda i: mapper[i]) + for k, g in groupby(data, key=lambda i: mapper[i]): + groups.append(tuple(sorted(g, key=lambda i: exprs.index(i)))) + + # Sanity check + assert max(mapper.values()) == len(groups) + + return tuple(groups) diff --git a/devito/operator.py b/devito/operator.py index 8e5123c6d6..26e3988ffc 100644 --- a/devito/operator.py +++ b/devito/operator.py @@ -2,7 +2,6 @@ from collections import OrderedDict from operator import attrgetter -from cached_property import cached_property import ctypes import numpy as np @@ -14,7 +13,7 @@ from devito.dle import transform from devito.dse import rewrite from devito.exceptions import InvalidOperator -from devito.function import Forward, Backward, Constant +from devito.function import Constant from devito.logger import bar, info from devito.ir.equations import LoweredEq from devito.ir.clusters import clusterize @@ -43,8 +42,6 @@ class Operator(Callable): * name : Name of the kernel function - defaults to "Kernel". * subs : Dict or list of dicts containing SymPy symbol substitutions for each expression respectively. - * time_axis : :class:`TimeAxis` object to indicate direction in which - to advance time during computation. * dse : Use the Devito Symbolic Engine to optimize the expressions - defaults to ``configuration['dse']``. * dle : Use the Devito Loop Engine to optimize the loops - @@ -59,7 +56,6 @@ def __init__(self, expressions, **kwargs): self.name = kwargs.get("name", "Kernel") subs = kwargs.get("subs", {}) - time_axis = kwargs.get("time_axis", Forward) dse = kwargs.get("dse", configuration['dse']) dle = kwargs.get("dle", configuration['dle']) @@ -81,11 +77,6 @@ def __init__(self, expressions, **kwargs): self.dtype = retrieve_dtype(expressions) self.input, self.output, self.dimensions = retrieve_symbols(expressions) - # Set the direction of time acoording to the given TimeAxis - for time in [d for d in self.dimensions if d.is_Time]: - if not time.is_Derived: - time.reverse = time_axis == Backward - # Group expressions based on their iteration space and data dependences, # and apply the Devito Symbolic Engine (DSE) for flop optimization clusters = clusterize(expressions) @@ -124,31 +115,32 @@ def __init__(self, expressions, **kwargs): # Finish instantiation super(Operator, self).__init__(self.name, nodes, 'int', parameters, ()) - @cached_property - def _argument_defaults(self): + def _argument_defaults(self, arguments): """ Derive all default values from parameters and ensure uniqueness. """ default_args = ArgumentMap() for p in self.input: - default_args.update(p.argument_defaults()) + if p.name not in arguments: + default_args.update(p.argument_defaults()) for p in self.dimensions: - if p.is_Sub: + if p.name not in arguments and p.is_Sub: default_args.update(p.argument_defaults(default_args)) - return {k: default_args.reduce(k) for k in default_args} + return {k: default_args.reduce(k) for k in default_args if k not in arguments} def arguments(self, **kwargs): """ Process runtime arguments passed to ``.apply()` and derive default values for any remaining arguments. """ - # First, derive all default values from parameters - arguments = self._argument_defaults.copy() - - # Next, we insert user-provided overrides + arguments = {} + # First, we insert user-provided override for p in self.input + self.dimensions: arguments.update(p.argument_values(**kwargs)) + # Second, derive all remaining default values from parameters + arguments.update(self._argument_defaults(arguments)) + # Derive additional values for DLE arguments # TODO: This is not pretty, but it works for now. Ideally, the # DLE arguments would be massaged into the IET so as to comply diff --git a/devito/region.py b/devito/region.py deleted file mode 100644 index b1d5dd0d6c..0000000000 --- a/devito/region.py +++ /dev/null @@ -1,29 +0,0 @@ -__all__ = ['DOMAIN', 'INTERIOR'] - - -class Region(object): - """ - A region of the computational domain over which a :class:`Function` is - discretized. - """ - - def __init__(self, name): - self._name = name - - def __repr__(self): - return self._name - - def __eq__(self, other): - return isinstance(other, Region) and self._name == other._name - - -DOMAIN = Region('DOMAIN') -""" -Represent the physical domain of the PDE; that is, domain = boundary + interior -""" - -INTERIOR = Region('INTERIOR') -""" -Represent the physical interior domain of the PDE; that is, PDE boundaries are -not included. -""" diff --git a/devito/symbolics/extended_sympy.py b/devito/symbolics/extended_sympy.py index 639711f6bf..819add281f 100644 --- a/devito/symbolics/extended_sympy.py +++ b/devito/symbolics/extended_sympy.py @@ -7,10 +7,9 @@ from sympy.core.basic import _aresame from sympy.functions.elementary.trigonometric import TrigonometricFunction -from devito.region import DOMAIN from devito.tools import as_tuple -__all__ = ['FrozenExpr', 'Eq', 'CondEq', 'CondNe', 'Inc', 'Mul', 'Add', 'IntDiv', +__all__ = ['FrozenExpr', 'Eq', 'CondEq', 'CondNe', 'Mul', 'Add', 'IntDiv', 'FunctionFromPointer', 'ListInitializer', 'taylor_sin', 'taylor_cos', 'bhaskara_sin', 'bhaskara_cos'] @@ -48,14 +47,8 @@ class Eq(sympy.Eq, FrozenExpr): """A customized version of :class:`sympy.Eq` which suppresses evaluation.""" - is_Increment = False - def __new__(cls, *args, **kwargs): - kwargs['evaluate'] = False - region = kwargs.pop('region', DOMAIN) - obj = sympy.Eq.__new__(cls, *args, **kwargs) - obj._region = region - return obj + return sympy.Eq.__new__(cls, *args, evaluate=False) class CondEq(sympy.Eq, FrozenExpr): @@ -63,8 +56,7 @@ class CondEq(sympy.Eq, FrozenExpr): equality. It suppresses evaluation.""" def __new__(cls, *args, **kwargs): - kwargs['evaluate'] = False - return sympy.Eq.__new__(cls, *args, **kwargs) + return sympy.Eq.__new__(cls, *args, evaluate=False) class CondNe(sympy.Ne, FrozenExpr): @@ -72,17 +64,7 @@ class CondNe(sympy.Ne, FrozenExpr): inequality. It suppresses evaluation.""" def __new__(cls, *args, **kwargs): - kwargs['evaluate'] = False - return sympy.Ne.__new__(cls, *args, **kwargs) - - -class Inc(Eq): - """ - A special :class:`Eq` carrying the information that a linear increment is - performed. - """ - - is_Increment = True + return sympy.Ne.__new__(cls, *args, evaluate=False) class Mul(sympy.Mul, FrozenExpr): diff --git a/devito/symbolics/inspection.py b/devito/symbolics/inspection.py index 629b260541..3d127bf130 100644 --- a/devito/symbolics/inspection.py +++ b/devito/symbolics/inspection.py @@ -138,4 +138,9 @@ def dimension_sort(expr, key=None): dimensions = filter_sorted(dimensions, key=attrgetter('name')) # for determinism ordering.extend([i for i in dimensions if i not in ordering]) + # Add parent dimensions + derived = [i for i in ordering if i.is_Derived] + for i in derived: + ordering.insert(ordering.index(i), i.parent) + return sorted(ordering, key=lambda i: not i.is_Time) diff --git a/devito/yask/operator.py b/devito/yask/operator.py index 7f033f1ed8..06977eed1d 100644 --- a/devito/yask/operator.py +++ b/devito/yask/operator.py @@ -1,6 +1,5 @@ from __future__ import absolute_import -from cached_property import cached_property import cgen as c import numpy as np from sympy import Indexed @@ -116,10 +115,8 @@ def _build_casts(self, nodes): return List(body=casts + [nodes]) - @cached_property - def _argument_defaults(self): - default_args = super(Operator, self)._argument_defaults - + def _argument_defaults(self, arguments): + default_args = super(Operator, self)._argument_defaults(arguments) # Add in solution pointer default_args[namespace['code-soln-name']] = self.yk_soln.rawpointer diff --git a/examples/checkpointing/checkpointing_example.py b/examples/checkpointing/checkpointing_example.py index 3124868eeb..ddf44dace2 100644 --- a/examples/checkpointing/checkpointing_example.py +++ b/examples/checkpointing/checkpointing_example.py @@ -12,7 +12,7 @@ class CheckpointingExample(GradientExample): @cached_property def forward_field(self): return TimeFunction(name="u", grid=self.model.grid, time_order=2, - space_order=self.space_order, save=False) + space_order=self.space_order) @cached_property def forward_operator(self): diff --git a/examples/seismic/acoustic/acoustic_example.py b/examples/seismic/acoustic/acoustic_example.py index a40ee3ebaf..5e5a0d63ac 100644 --- a/examples/seismic/acoustic/acoustic_example.py +++ b/examples/seismic/acoustic/acoustic_example.py @@ -11,7 +11,7 @@ def smooth10(vel, shape): if np.isscalar(vel): return .9 * vel * np.ones(shape, dtype=np.float32) - out = np.ones(shape, dtype=np.float32) + out = np.ones(shape, dtype=vel.dtype) nz = shape[-1] for a in range(5, nz-6): @@ -28,7 +28,7 @@ def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), constant=False, **kwargs): nrec = shape[0] preset = 'constant-isotropic' if constant else 'layers-isotropic' - model = demo_model(preset, shape=shape, + model = demo_model(preset, shape=shape, dtype=kwargs.pop('dtype', np.float32), spacing=spacing, nbpml=nbpml) # Derive timestepping from model spacing @@ -43,7 +43,7 @@ def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), src.coordinates.data[0, -1] = model.origin[-1] + 2 * spacing[-1] # Define receiver geometry (spread across x, just below surface) - rec = Receiver(name='nrec', grid=model.grid, ntime=nt, npoint=nrec) + rec = Receiver(name='rec', grid=model.grid, ntime=nt, npoint=nrec) rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec) rec.coordinates.data[:, 1:] = src.coordinates.data[0, 1:] diff --git a/examples/seismic/acoustic/gradient_example.py b/examples/seismic/acoustic/gradient_example.py index 8333e86265..d9c43f8de5 100644 --- a/examples/seismic/acoustic/gradient_example.py +++ b/examples/seismic/acoustic/gradient_example.py @@ -36,7 +36,7 @@ def _setup_model_and_acquisition(self, shape, spacing, nbpml, tn): # Define receiver geometry (spread across x, just below surface) # We need two receiver fields - one for the true (verification) run - rec_t = Receiver(name='rec_t', grid=model.grid, ntime=self.nt, npoint=nrec) + rec_t = Receiver(name='rec', grid=model.grid, ntime=self.nt, npoint=nrec) rec_t.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec) rec_t.coordinates.data[:, 1:] = src.coordinates.data[0, 1:] @@ -47,7 +47,7 @@ def _setup_model_and_acquisition(self, shape, spacing, nbpml, tn): coordinates=rec_t.coordinates.data) # Receiver for Gradient - self.rec_g = Receiver(name="rec_g", coordinates=self.rec.coordinates.data, + self.rec_g = Receiver(name="rec", coordinates=self.rec.coordinates.data, grid=model.grid, dt=self.dt, ntime=self.nt) # Gradient symbol diff --git a/examples/seismic/acoustic/operators.py b/examples/seismic/acoustic/operators.py index 76ffda4846..1395883d69 100644 --- a/examples/seismic/acoustic/operators.py +++ b/examples/seismic/acoustic/operators.py @@ -1,6 +1,6 @@ from sympy import solve, Symbol -from devito import Eq, Operator, Forward, Backward, Function, TimeFunction +from devito import Eq, Operator, Function, TimeFunction from devito.logger import error from examples.seismic import PointSource, Receiver @@ -86,7 +86,7 @@ def ForwardOperator(model, source, receiver, space_order=4, # Substitute spacing terms to reduce flops return Operator(eqn + src_term + rec_term, subs=model.spacing_map, - time_axis=Forward, name='Forward', **kwargs) + name='Forward', **kwargs) def AdjointOperator(model, source, receiver, space_order=4, @@ -121,7 +121,7 @@ def AdjointOperator(model, source, receiver, space_order=4, # Substitute spacing terms to reduce flops return Operator(eqn + receivers + source_a, subs=model.spacing_map, - time_axis=Backward, name='Adjoint', **kwargs) + name='Adjoint', **kwargs) def GradientOperator(model, source, receiver, space_order=4, save=True, @@ -162,7 +162,7 @@ def GradientOperator(model, source, receiver, space_order=4, save=True, # Substitute spacing terms to reduce flops return Operator(eqn + receivers + [gradient_update], subs=model.spacing_map, - time_axis=Backward, name='Gradient', **kwargs) + name='Gradient', **kwargs) def BornOperator(model, source, receiver, space_order=4, @@ -204,4 +204,4 @@ def BornOperator(model, source, receiver, space_order=4, # Substitute spacing terms to reduce flops return Operator(eqn1 + source + eqn2 + receivers, subs=model.spacing_map, - time_axis=Forward, name='Born', **kwargs) + name='Born', **kwargs) diff --git a/examples/seismic/acoustic/wavesolver.py b/examples/seismic/acoustic/wavesolver.py index 670b49ceab..752303daa0 100644 --- a/examples/seismic/acoustic/wavesolver.py +++ b/examples/seismic/acoustic/wavesolver.py @@ -42,7 +42,7 @@ def __init__(self, model, source, receiver, kernel='OT2', space_order=2, **kwarg self._kwargs = kwargs @memoized_meth - def op_fwd(self, save=False): + def op_fwd(self, save=None): """Cached operator for forward runs with buffered wavefield""" return ForwardOperator(self.model, save=save, source=self.source, receiver=self.receiver, kernel=self.kernel, @@ -51,7 +51,7 @@ def op_fwd(self, save=False): @memoized_meth def op_adj(self): """Cached operator for adjoint runs""" - return AdjointOperator(self.model, save=False, source=self.source, + return AdjointOperator(self.model, save=None, source=self.source, receiver=self.receiver, kernel=self.kernel, space_order=self.space_order, **self._kwargs) @@ -65,11 +65,11 @@ def op_grad(self): @memoized_meth def op_born(self): """Cached operator for born runs""" - return BornOperator(self.model, save=False, source=self.source, + return BornOperator(self.model, save=None, source=self.source, receiver=self.receiver, kernel=self.kernel, space_order=self.space_order, **self._kwargs) - def forward(self, src=None, rec=None, u=None, m=None, save=False, **kwargs): + def forward(self, src=None, rec=None, u=None, m=None, save=None, **kwargs): """ Forward modelling function that creates the necessary data objects for running a forward modelling operator. @@ -128,7 +128,7 @@ def adjoint(self, rec, srca=None, v=None, m=None, **kwargs): # Create the adjoint wavefield if not provided if v is None: - v = TimeFunction(name='v', grid=self.model.grid, save=None, + v = TimeFunction(name='v', grid=self.model.grid, time_order=2, space_order=self.space_order) # Pick m from model unless explicitly provided @@ -160,7 +160,7 @@ def gradient(self, rec, u, v=None, grad=None, m=None, **kwargs): # Create the forward wavefield if v is None: - v = TimeFunction(name='v', grid=self.model.grid, save=None, + v = TimeFunction(name='v', grid=self.model.grid, time_order=2, space_order=self.space_order) # Pick m from model unless explicitly provided @@ -193,10 +193,10 @@ def born(self, dmin, src=None, rec=None, u=None, U=None, m=None, **kwargs): # Create the forward wavefields u and U if not provided if u is None: - u = TimeFunction(name='u', grid=self.model.grid, save=False, + u = TimeFunction(name='u', grid=self.model.grid, time_order=2, space_order=self.space_order) if U is None: - U = TimeFunction(name='U', grid=self.model.grid, save=False, + U = TimeFunction(name='U', grid=self.model.grid, time_order=2, space_order=self.space_order) # Pick m from model unless explicitly provided diff --git a/examples/seismic/model.py b/examples/seismic/model.py index 79cd6fa345..19098204e9 100644 --- a/examples/seismic/model.py +++ b/examples/seismic/model.py @@ -44,9 +44,10 @@ def demo_model(preset, **kwargs): spacing = kwargs.pop('spacing', tuple([10. for _ in shape])) origin = kwargs.pop('origin', tuple([0. for _ in shape])) nbpml = kwargs.pop('nbpml', 10) + dtype = kwargs.pop('dtype', np.float32) vp = kwargs.pop('vp', 1.5) - return Model(vp=vp, origin=origin, shape=shape, + return Model(vp=vp, origin=origin, shape=shape, dtype=dtype, spacing=spacing, nbpml=nbpml, **kwargs) elif preset.lower() in ['constant-tti']: @@ -56,16 +57,17 @@ def demo_model(preset, **kwargs): spacing = kwargs.pop('spacing', tuple([10. for _ in shape])) origin = kwargs.pop('origin', tuple([0. for _ in shape])) nbpml = kwargs.pop('nbpml', 10) - v = np.empty(shape, dtype=np.float32) + dtype = kwargs.pop('dtype', np.float32) + v = np.empty(shape, dtype=dtype) v[:] = 1.5 - epsilon = .3*np.ones(shape) - delta = .2*np.ones(shape) - theta = .7*np.ones(shape) + epsilon = .3*np.ones(shape, dtype=dtype) + delta = .2*np.ones(shape, dtype=dtype) + theta = .7*np.ones(shape, dtype=dtype) phi = None if len(shape) > 2: - phi = .35*np.ones(shape) + phi = .35*np.ones(shape, dtype=dtype) - return Model(vp=v, origin=origin, shape=shape, + return Model(vp=v, origin=origin, shape=shape, dtype=dtype, spacing=spacing, nbpml=nbpml, epsilon=epsilon, delta=delta, theta=theta, phi=phi, **kwargs) @@ -79,17 +81,18 @@ def demo_model(preset, **kwargs): shape = kwargs.pop('shape', (101, 101)) spacing = kwargs.pop('spacing', tuple([10. for _ in shape])) origin = kwargs.pop('origin', tuple([0. for _ in shape])) + dtype = kwargs.pop('dtype', np.float32) nbpml = kwargs.pop('nbpml', 10) ratio = kwargs.pop('ratio', 2) vp_top = kwargs.pop('vp_top', 1.5) vp_bottom = kwargs.pop('vp_bottom', 2.5) # Define a velocity profile in km/s - v = np.empty(shape, dtype=np.float32) + v = np.empty(shape, dtype=dtype) v[:] = vp_top # Top velocity (background) v[..., int(shape[-1] / ratio):] = vp_bottom # Bottom velocity - return Model(vp=v, origin=origin, shape=shape, + return Model(vp=v, origin=origin, shape=shape, dtype=dtype, spacing=spacing, nbpml=nbpml, **kwargs) elif preset.lower() in ['layers-tti', 'twolayer-tti', '2layer-tti']: @@ -100,13 +103,14 @@ def demo_model(preset, **kwargs): shape = kwargs.pop('shape', (101, 101)) spacing = kwargs.pop('spacing', tuple([10. for _ in shape])) origin = kwargs.pop('origin', tuple([0. for _ in shape])) + dtype = kwargs.pop('dtype', np.float32) nbpml = kwargs.pop('nbpml', 10) ratio = kwargs.pop('ratio', 2) vp_top = kwargs.pop('vp_top', 1.5) vp_bottom = kwargs.pop('vp_bottom', 2.5) # Define a velocity profile in km/s - v = np.empty(shape, dtype=np.float32) + v = np.empty(shape, dtype=dtype) v[:] = vp_top # Top velocity (background) v[..., int(shape[-1] / ratio):] = vp_bottom # Bottom velocity @@ -117,7 +121,7 @@ def demo_model(preset, **kwargs): if len(shape) > 2: phi = .1*(v - 1.5) - return Model(vp=v, origin=origin, shape=shape, + return Model(vp=v, origin=origin, shape=shape, dtype=dtype, spacing=spacing, nbpml=nbpml, epsilon=epsilon, delta=delta, theta=theta, phi=phi, **kwargs) @@ -126,6 +130,7 @@ def demo_model(preset, **kwargs): # A simple circle in a 2D domain with a background velocity. # By default, the circle velocity is 2.5 km/s, # and the background veloity is 3.0 km/s. + dtype = kwargs.pop('dtype', np.float32) shape = kwargs.pop('shape', (101, 101)) spacing = kwargs.pop('spacing', tuple([10. for _ in shape])) origin = kwargs.pop('origin', tuple([0. for _ in shape])) @@ -137,15 +142,15 @@ def demo_model(preset, **kwargs): # Only a 2D preset is available currently assert(len(shape) == 2) - v = np.empty(shape, dtype=np.float32) + v = np.empty(shape, dtype=dtype) v[:] = vp_background a, b = shape[0] / 2, shape[1] / 2 y, x = np.ogrid[-a:shape[0]-a, -b:shape[1]-b] v[x*x + y*y <= r*r] = vp - return Model(vp=v, origin=origin, shape=shape, - spacing=spacing, nbpml=nbpml) + return Model(vp=v, origin=origin, shape=shape, dtype=dtype, + spacing=spacing, nbpml=nbpml, **kwargs) elif preset.lower() in ['marmousi-isotropic', 'marmousi2d-isotropic']: shape = (1601, 401) @@ -165,8 +170,8 @@ def demo_model(preset, **kwargs): # Cut the model to make it slightly cheaper v = v[301:-300, :] - return Model(vp=v, origin=origin, shape=v.shape, - spacing=spacing, nbpml=20) + return Model(vp=v, origin=origin, shape=v.shape, dtype=np.float32, + spacing=spacing, nbpml=20, **kwargs) elif preset.lower() in ['marmousi-tti2d', 'marmousi2d-tti']: @@ -205,7 +210,7 @@ def demo_model(preset, **kwargs): theta = np.float32(np.pi / 180 * theta.reshape(shape_full)) theta = theta[101, :, :] - return Model(vp=vp, origin=origin, shape=shape, + return Model(vp=vp, origin=origin, shape=shape, dtype=np.float32, spacing=spacing, nbpml=nbpml, epsilon=epsilon, delta=delta, theta=theta, **kwargs) @@ -246,7 +251,7 @@ def demo_model(preset, **kwargs): dtype='float32', sep="") phi = np.float32(np.pi / 180 * phi.reshape(shape)) - return Model(vp=vp, origin=origin, shape=shape, + return Model(vp=vp, origin=origin, shape=shape, dtype=np.float32, spacing=spacing, nbpml=nbpml, epsilon=epsilon, delta=delta, theta=theta, phi=phi, **kwargs) diff --git a/examples/seismic/tti/operators.py b/examples/seismic/tti/operators.py index 678c28f0b5..6569052d91 100644 --- a/examples/seismic/tti/operators.py +++ b/examples/seismic/tti/operators.py @@ -375,10 +375,10 @@ def ForwardOperator(model, source, receiver, space_order=4, # Create symbols for forward wavefield, source and receivers u = TimeFunction(name='u', grid=model.grid, - save=save, time_dim=source.nt if save else None, + save=source.nt if save else None, time_order=2, space_order=space_order) v = TimeFunction(name='v', grid=model.grid, - save=save, time_dim=source.nt if save else None, + save=source.nt if save else None, time_order=2, space_order=space_order) src = PointSource(name='src', grid=model.grid, ntime=source.nt, npoint=source.npoint) diff --git a/examples/seismic/tti/wavesolver.py b/examples/seismic/tti/wavesolver.py index 7b33d13e0d..d9c2c83c4a 100644 --- a/examples/seismic/tti/wavesolver.py +++ b/examples/seismic/tti/wavesolver.py @@ -68,13 +68,13 @@ def forward(self, src=None, rec=None, u=None, v=None, m=None, # Create the forward wavefield if not provided if u is None: - u = TimeFunction(name='u', grid=self.model.grid, save=save, - time_dim=self.source.nt if save else None, + u = TimeFunction(name='u', grid=self.model.grid, + save=self.source.nt if save else None, time_order=2, space_order=self.space_order) # Create the forward wavefield if not provided if v is None: - v = TimeFunction(name='v', grid=self.model.grid, save=save, - time_dim=self.source.nt if save else None, + v = TimeFunction(name='v', grid=self.model.grid, + save=self.source.nt if save else None, time_order=2, space_order=self.space_order) # Pick m from model unless explicitly provided if m is None: diff --git a/examples/seismic/tutorials/02_rtm.ipynb b/examples/seismic/tutorials/02_rtm.ipynb index 31b3a20245..7a0b03d72f 100644 --- a/examples/seismic/tutorials/02_rtm.ipynb +++ b/examples/seismic/tutorials/02_rtm.ipynb @@ -366,7 +366,7 @@ "outputs": [], "source": [ "# Define gradient operator for imaging\n", - "from devito import TimeFunction, Backward, Operator\n", + "from devito import TimeFunction, Operator\n", "from examples.seismic import PointSource\n", "\n", "from sympy import solve, Eq\n", @@ -394,7 +394,7 @@ " image_update = Eq(image, image - u * v)\n", "\n", " return Operator([stencil] + res_term + [image_update],\n", - " subs=model.spacing_map, time_axis=Backward)" + " subs=model.spacing_map)" ] }, { diff --git a/tests/test_adjointA.py b/tests/test_adjointA.py index 1595e24051..a9426107df 100644 --- a/tests/test_adjointA.py +++ b/tests/test_adjointA.py @@ -33,7 +33,7 @@ def test_acoustic(mkey, shape, kernel, space_order, nbpml): nrec = 130 # Number of receivers # Create model from preset - model = demo_model(spacing=[15. for _ in shape], + model = demo_model(spacing=[15. for _ in shape], dtype=np.float64, shape=shape, nbpml=nbpml, **(presets[mkey])) # Derive timestepping from model spacing @@ -47,7 +47,7 @@ def test_acoustic(mkey, shape, kernel, space_order, nbpml): src.coordinates.data[0, -1] = 30. # Define receiver geometry (same as source, but spread across x) - rec = Receiver(name='nrec', grid=model.grid, ntime=nt, npoint=nrec) + rec = Receiver(name='rec', grid=model.grid, ntime=nt, npoint=nrec) rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec) rec.coordinates.data[:, 1:] = src.coordinates.data[0, 1:] @@ -68,5 +68,5 @@ def test_acoustic(mkey, shape, kernel, space_order, nbpml): term1 = np.dot(srca.data.reshape(-1), solver.source.data) term2 = linalg.norm(rec.data) ** 2 info(': %f, : %f, difference: %12.12f, ratio: %f' - % (term1, term2, term1 - term2, term1 / term2)) - assert np.isclose(term1, term2, rtol=1.e-5) + % (term1, term2, (term1 - term2)/term1, term1 / term2)) + assert np.isclose((term1 - term2)/term1, 0., rtol=1.e-10) diff --git a/tests/test_adjointJ.py b/tests/test_adjointJ.py index 7b358190d4..bf6cd17188 100644 --- a/tests/test_adjointJ.py +++ b/tests/test_adjointJ.py @@ -20,7 +20,7 @@ def test_acousticJ(shape, space_order): # Create two-layer "true" model from preset with a fault 1/3 way down model = demo_model('layers-isotropic', ratio=3, vp_top=1.5, vp_bottom=2.5, - spacing=spacing, shape=shape, nbpml=nbpml) + spacing=spacing, shape=shape, nbpml=nbpml, dtype=np.float64) # Derive timestepping from model spacing dt = model.critical_dt @@ -43,7 +43,7 @@ def test_acousticJ(shape, space_order): # Create initial model (m0) with a constant velocity throughout model0 = demo_model('layers-isotropic', ratio=3, vp_top=1.5, vp_bottom=1.5, - spacing=spacing, shape=shape, nbpml=nbpml) + spacing=spacing, shape=shape, nbpml=nbpml, dtype=np.float64) # Compute the full wavefield u0 _, u0, _ = solver.forward(save=True, m=model0.m) diff --git a/tests/test_checkpointing.py b/tests/test_checkpointing.py index 152a604cf5..07b2358054 100644 --- a/tests/test_checkpointing.py +++ b/tests/test_checkpointing.py @@ -7,7 +7,7 @@ import pytest from functools import reduce -from devito import Grid, TimeFunction, Operator, Backward, Function, Eq, silencio +from devito import Grid, TimeFunction, Operator, Function, Eq, silencio @silencio(log_level='WARNING') @@ -208,7 +208,7 @@ def test_index_alignment(const): and hence grad = 0*0 + 1*1 + 2*2 + 3*3 = sum(n^2) where n -> [0, nt] If the test fails, the resulting number can tell you how the fields are misaligned """ - nt = 10 + nt = 9 grid = Grid(shape=(3, 5)) order_of_eqn = 1 modulo_factor = order_of_eqn + 1 @@ -228,7 +228,7 @@ def test_index_alignment(const): v.data[last_time_step_v, :, :] = u.data[last_time_step_u, :, :] # Decrement one in the reverse pass 3 -> 2 -> 1 -> 0 adj_eqn = Eq(v.backward, v - 1.*const) - adj_op = Operator(adj_eqn, time_axis=Backward) + adj_op = Operator(adj_eqn) # Invocation 2 adj_op(t=nt, constant=1) @@ -241,7 +241,7 @@ def test_index_alignment(const): # Multiply u and v and add them # = 3*3 + 2*2 + 1*1 + 0*0 prod_eqn = Eq(prod, prod + u * v) - comb_op = Operator([adj_eqn, prod_eqn], time_axis=Backward) + comb_op = Operator([adj_eqn, prod_eqn]) # Invocation 3 comb_op(time=nt, constant=1) @@ -262,7 +262,7 @@ def test_index_alignment(const): time_order=order_of_eqn) prod_eqn_2 = Eq(prod, prod + u_nosave * v) - comb_op_2 = Operator([adj_eqn, prod_eqn_2], time_axis=Backward) + comb_op_2 = Operator([adj_eqn, prod_eqn_2]) wrap_rev = CheckpointOperator(comb_op_2, constant=1, time_order=order_of_eqn) wrp = Revolver(cp, wrap_fw, wrap_rev, None, nt-order_of_eqn) diff --git a/tests/test_dimension.py b/tests/test_dimension.py index ebc9f0c23f..741d837759 100644 --- a/tests/test_dimension.py +++ b/tests/test_dimension.py @@ -83,6 +83,7 @@ def test_conditional_shifted(): Eq(usave.subs(time_subsampled, time_subsampled - t_sub_shift), u)] op = Operator(eqns) + # Starting at time_s=10, so time_subsampled - t_sub_shift is in range op.apply(time_s=10, time_e=nt, t_sub_shift=3) assert np.all(np.allclose(u.data[0], 8)) assert np.all([np.allclose(u2.data[i], i - 10) for i in range(10, nt)]) diff --git a/tests/test_dle.py b/tests/test_dle.py index 2c74cdd0b8..09865e2ad1 100644 --- a/tests/test_dle.py +++ b/tests/test_dle.py @@ -12,6 +12,7 @@ from devito.dle import transform from devito.dle.backends import DevitoRewriter as Rewriter from devito import Grid, Function, TimeFunction, Eq, Operator +from devito.ir.equations import LoweredEq from devito.ir.iet import (ELEMENTAL, Expression, Callable, Iteration, List, tagger, Transformer, FindNodes, iet_analyze, retrieve_iteration_tree) @@ -391,7 +392,7 @@ def test_cache_blocking_edge_cases_highorder(shape, blockshape): ]) def test_loops_ompized(fa, fb, fc, fd, t0, t1, t2, t3, exprs, expected, iters): scope = [fa, fb, fc, fd, t0, t1, t2, t3] - node_exprs = [Expression(EVAL(i, *scope)) for i in exprs] + node_exprs = [Expression(LoweredEq(EVAL(i, *scope))) for i in exprs] ast = iters[6](iters[7](node_exprs)) ast = iet_analyze(ast) diff --git a/tests/test_gradient.py b/tests/test_gradient.py index cfcd9f2587..64ba8f32bf 100644 --- a/tests/test_gradient.py +++ b/tests/test_gradient.py @@ -34,7 +34,7 @@ def test_gradientFWI(shape, kernel, space_order): :return: assertion that the Taylor properties are satisfied """ spacing = tuple(15. for _ in shape) - wave = setup(shape=shape, spacing=spacing, + wave = setup(shape=shape, spacing=spacing, dtype=np.float64, kernel=kernel, space_order=space_order, nbpml=10+space_order/2) m0 = smooth10(wave.model.m.data, wave.model.shape_domain) @@ -97,11 +97,11 @@ def test_gradientJ(shape, kernel, space_order): :return: assertion that the Taylor properties are satisfied """ spacing = tuple(15. for _ in shape) - wave = setup(shape=shape, spacing=spacing, + wave = setup(shape=shape, spacing=spacing, dtype=np.float64, kernel=kernel, space_order=space_order, tn=1000., nbpml=10+space_order/2) m0 = smooth10(wave.model.m.data, wave.model.shape_domain) - dm = np.float32(wave.model.m.data - m0) + dm = np.float64(wave.model.m.data - m0) linrec = Receiver(name='rec', grid=wave.model.grid, ntime=wave.receiver.nt, coordinates=wave.receiver.coordinates.data) # Compute receiver data and full wavefield for the smooth velocity @@ -133,4 +133,4 @@ def test_gradientJ(shape, kernel, space_order): if __name__ == "__main__": - test_gradientJ(shape=(60, 70), kernel='OT2', space_order=4) + test_gradientFWI(shape=(60, 70), kernel='OT2', space_order=4) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index c96e9adbe5..cc8191ed37 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -208,7 +208,7 @@ def test_position(shape): src.coordinates.data[0, -1] = 30. # Define receiver geometry (same as source, but spread across x) - rec = Receiver(name='nrec', grid=model.grid, ntime=nt, npoint=nrec) + rec = Receiver(name='rec', grid=model.grid, ntime=nt, npoint=nrec) rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec) rec.coordinates.data[:, 1:] = src.coordinates.data[0, 1:] @@ -225,7 +225,7 @@ def test_position(shape): src.coordinates.data[0, -1] = 130. # Define receiver geometry (same as source, but spread across x) - rec2 = Receiver(name='rec2', grid=model.grid, ntime=nt, npoint=nrec) + rec2 = Receiver(name='rec', grid=model.grid, ntime=nt, npoint=nrec) rec2.coordinates.data[:, 0] = np.linspace(100., 100. + model.domain_size[0], num=nrec) rec2.coordinates.data[:, 1:] = src.coordinates.data[0, 1:] diff --git a/tests/test_ir.py b/tests/test_ir.py index dcc2194744..c73a758ad0 100644 --- a/tests/test_ir.py +++ b/tests/test_ir.py @@ -3,9 +3,11 @@ from conftest import EVAL, time, x, y, z, skipif_yask # noqa from devito import Eq # noqa +from devito.ir.equations import LoweredEq from devito.ir.iet.nodes import Conditional, Expression from devito.ir.support.basic import IterationInstance, TimedAccess, Scope -from devito.ir.support.space import NullInterval, Interval, Space +from devito.ir.support.space import (NullInterval, Interval, IntervalGroup, + Forward, Backward) @pytest.fixture(scope="session") @@ -27,14 +29,14 @@ def ii_literal(fa, fc): @pytest.fixture(scope="session") def ta_literal(fc): - tcxy_w0 = TimedAccess(fc[x, y], 'W', 0) - tcxy_r0 = TimedAccess(fc[x, y], 'R', 0) - tcx1y1_r1 = TimedAccess(fc[x + 1, y + 1], 'R', 1) - tcx1y_r1 = TimedAccess(fc[x + 1, y], 'R', 1) - x.reverse = True - rev_tcxy_w0 = TimedAccess(fc[x, y], 'W', 0) - rev_tcx1y1_r1 = TimedAccess(fc[x + 1, y + 1], 'R', 1) - x.reverse = False + fwd_directions = {x: Forward, y: Forward} + mixed_directions = {x: Backward, y: Forward} + tcxy_w0 = TimedAccess(fc[x, y], 'W', 0, fwd_directions) + tcxy_r0 = TimedAccess(fc[x, y], 'R', 0, fwd_directions) + tcx1y1_r1 = TimedAccess(fc[x + 1, y + 1], 'R', 1, fwd_directions) + tcx1y_r1 = TimedAccess(fc[x + 1, y], 'R', 1, fwd_directions) + rev_tcxy_w0 = TimedAccess(fc[x, y], 'W', 0, mixed_directions) + rev_tcx1y1_r1 = TimedAccess(fc[x + 1, y + 1], 'R', 1, mixed_directions) return tcxy_w0, tcxy_r0, tcx1y1_r1, tcx1y_r1, rev_tcxy_w0, rev_tcx1y1_r1 @@ -203,7 +205,11 @@ def test_dependences_eq(expr, expected, ti0, ti1, fa): * the dimension causing the dependence * whether it's direct or indirect (i.e., through A[B[i]]) """ - expr = EVAL(expr, ti0.base, ti1.base, fa) + expr = LoweredEq(EVAL(expr, ti0.base, ti1.base, fa)) + + # Force innatural flow, only to stress the compiler to see if it was + # capable of detecting anti-dependences + expr.ispace.directions = {i: Forward for i in expr.ispace.directions} scope = Scope(expr) deps = scope.d_all @@ -300,9 +306,14 @@ def test_dependences_scope(exprs, expected, ti0, ti1, ti3, fa): * if it's a flow, anti, or output dependence * the dimension causing the dependence """ - exprs = EVAL(exprs, ti0.base, ti1.base, ti3.base, fa) + exprs = [LoweredEq(i) for i in EVAL(exprs, ti0.base, ti1.base, ti3.base, fa)] expected = [tuple(i.split(',')) for i in expected] + # Force innatural flow, only to stress the compiler to see if it was + # capable of detecting anti-dependences + for i in exprs: + i.ispace.directions = {i: Forward for i in i.ispace.directions} + scope = Scope(exprs) assert len(scope.d_all) == len(expected) @@ -379,13 +390,13 @@ def test_intervals_union(): nully = NullInterval(y) iy = Interval(y, -2, 2) - # Mixed disjoint (note: Space input order is irrelevant) - assert ix.union(ix4) == Space([ix4, ix]) + # Mixed disjoint (note: IntervalGroup input order is irrelevant) + assert ix.union(ix4) == IntervalGroup([ix4, ix]) assert ix.union(ix5) == Interval(x, -3, 2) - assert ix6.union(ix) == Space([ix, ix6]) - assert ix.union(nully) == Space([ix, nully]) - assert ix.union(iy) == Space([iy, ix]) - assert iy.union(ix) == Space([iy, ix]) + assert ix6.union(ix) == IntervalGroup([ix, ix6]) + assert ix.union(nully) == IntervalGroup([ix, nully]) + assert ix.union(iy) == IntervalGroup([iy, ix]) + assert iy.union(ix) == IntervalGroup([iy, ix]) @skipif_yask diff --git a/tests/test_operator.py b/tests/test_operator.py index 1e707ed001..67d09196fe 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -7,12 +7,13 @@ import numpy as np import pytest -from devito import (clear_cache, Grid, Eq, Operator, Constant, Function, Backward, - Forward, TimeFunction, SparseTimeFunction, Dimension, configuration, +from devito import (clear_cache, Grid, Eq, Operator, Constant, Function, + TimeFunction, SparseTimeFunction, Dimension, configuration, error, INTERIOR) from devito.foreign import Operator as OperatorForeign from devito.ir.iet import (Expression, Iteration, FindNodes, IsPerfectIteration, retrieve_iteration_tree) +from devito.ir.support import Any, Backward, Forward def dimify(dimensions): @@ -332,8 +333,8 @@ def test_default_composite_functions(self): expected = { 's': s.data, 's_coords': s.coordinates.data, # Default dimensions of the sparse data - 'p_size': 3, 'p_s': 0, 'p_e': 3, - 'd_size': 3, 'p_s': 0, 'p_e': 3, + 'p_s_size': 3, 'p_s_s': 0, 'p_s_e': 3, + 'd_size': 3, 'd_s': 0, 'd_e': 3, 'time_size': 4, 'time_s': 0, 'time_e': 4, } self.verify_arguments(op.arguments(), expected) @@ -471,13 +472,13 @@ def test_override_timefunction_data(self): assert (a.data[:] == 4.).all() # Override with symbol (different name) - a1 = TimeFunction(name='a1', grid=grid) + a1 = TimeFunction(name='a1', grid=grid, save=2) a1.data[:] = 2. op(t=2, a=a1) assert (a1.data[:] == 5.).all() # Override with symbol (same name as original) - a2 = TimeFunction(name='a', grid=grid) + a2 = TimeFunction(name='a', grid=grid, save=2) a2.data[:] = 3. op(t=2, a=a2) assert (a2.data[:] == 6.).all() @@ -774,7 +775,7 @@ def test_consistency_coupled_wo_ofs(self, tu, tv, ti0, t0, t1): ('Eq(ti0[x,y,z+2], ti0[x,y,z-1] + ti1[x,y,z+1])', 'Eq(ti1[x,y,z+3], ti3[x,y,z+1])', 'Eq(ti3[x,y,z+2], ti0[x,y,z+1]*ti3[x,y,z-1])'), - ('Eq(ti0[x,y,z], ti0[x-2,y-1,z-1] + ti1[x+2,y+3,z+1])', + ('Eq(ti0[x,y,z], ti0[x-2,y-1,z-1] + ti1[x-3,y+3,z+1])', 'Eq(ti1[x+4,y+5,z+3], ti3[x+1,y-4,z+1])', 'Eq(ti3[x+7,y,z+2], ti3[x+5,y,z-1] - ti0[x-3,y-2,z-4])') ]) @@ -799,75 +800,81 @@ def test_consistency_coupled_w_ofs(self, exprs, ti0, ti1, ti3): exprs = FindNodes(Expression).visit(tree[-1]) assert len(exprs) == 3 - @pytest.mark.parametrize('exprs,axis,expected,visit', [ + @pytest.mark.parametrize('exprs,directions,expected,visit', [ # WAR 2->3; expected=2 (('Eq(ti0[x,y,z], ti0[x,y,z] + ti1[x,y,z])', 'Eq(ti1[x,y,z], ti3[x,y,z])', 'Eq(ti3[x,y,z], ti1[x,y,z+1] + 1.)'), - Forward, ['xyz'], 'xyz'), + '**-', ['xyz'], 'xyz'), # WAR 1->2, 2->3; one may think it should be expected=3, but these are all # Arrays, so ti0 gets optimized through index bumping and array contraction, - # which results in expected=2 + # which results in expected=3 (('Eq(ti0[x,y,z], ti0[x,y,z] + ti1[x,y,z])', 'Eq(ti1[x,y,z], ti0[x,y,z+1])', - 'Eq(ti3[x,y,z], ti1[x,y,z+2] + 1.)'), - Forward, ['xyz', 'xyz'], 'xyzz'), + 'Eq(ti3[x,y,z], ti1[x,y,z-2] + 1.)'), + '*****', ['xyz', 'xyz', 'xyz'], 'xyzzz'), # WAR 1->3; expected=1 (('Eq(ti0[x,y,z], ti0[x,y,z] + ti1[x,y,z])', 'Eq(ti1[x,y,z], ti3[x,y,z])', 'Eq(ti3[x,y,z], ti0[x,y,z+1] + 1.)'), - Forward, ['xyz'], 'xyz'), + '**-', ['xyz'], 'xyz'), # WAR 1->2, 2->3; WAW 1->3; expected=2 # ti0 is an Array, so the observation made above still holds (expected=2 # rather than 3) (('Eq(ti0[x,y,z], ti0[x,y,z] + ti1[x,y,z])', 'Eq(ti1[x,y,z], 3*ti0[x,y,z+2])', 'Eq(ti0[x,y,0], ti0[x,y,0] + 1.)'), - Forward, ['xyz', 'xy'], 'xyz'), - # WAR 1->2; WAW 1->3; expected=3 + '**-', ['xyz', 'xy'], 'xyz'), + # WAR 1->2; WAW 1->3; expected=2 # Now tu, tv, tw are not Arrays, so they must end up in separate loops (('Eq(tu[t,x,y,z], tu[t,x,y,z] + tv[t,x,y,z])', 'Eq(tv[t,x,y,z], tu[t,x,y,z+2])', 'Eq(tu[t,x,y,0], tu[t,x,y,0] + 1.)'), - Forward, ['txyz', 'txyz', 'txy'], 'txyzz'), - # WAR 1->2; WAW 2->3; expected=3 + '***-', ['txyz', 'txy'], 'txyz'), + # WAR 1->2; RAW 2->3; expected=3 (('Eq(tu[t,x,y,z], tu[t,x,y,z] + tv[t,x,y,z])', 'Eq(tv[t,x,y,z], tu[t,x,y,z+2])', - 'Eq(tw[t,x,y,z], tv[t,x,y,z+3] + 1.)'), - Forward, ['txyz', 'txyz', 'txyz'], 'txyzzz'), - # WAR 1->2; WAW 1->3; expected=3 + 'Eq(tw[t,x,y,z], tv[t,x,y,z-3] + 1.)'), + '******', ['txyz', 'txyz', 'txyz'], 'txyzzz'), + # WAR 1->2; WAW 1->3; expected=2 (('Eq(tu[t,x,y,z], tu[t,x,y,z] + tv[t,x,y,z])', 'Eq(tv[t,x,y,z], tu[t,x+2,y,z])', 'Eq(tu[t,3,y,0], tu[t,3,y,0] + 1.)'), - Forward, ['txyz', 'txyz', 'ty'], 'txyzxyzy'), - # WAR 1->2, WAW 2->3; expected=3 + '*-***', ['txyz', 'ty'], 'txyzy'), + # RAW 1->2, WAR 2->3; expected=1 (('Eq(tu[t,x,y,z], tu[t,x,y,z] + tv[t,x,y,z])', - 'Eq(tv[t,x,y,z], tu[t,x,y,z+2])', + 'Eq(tv[t,x,y,z], tu[t,x,y,z-2])', 'Eq(tw[t,x,y,z], tv[t,x,y+1,z] + 1.)'), - Forward, ['txyz', 'txyz', 'txyz'], 'txyzzyz'), - # WAR 1->2; WAW 1->3; expected=3 - # Time Forward, anti dependence in time, end up in different loop nests + '**-+', ['txyz'], 'txyz'), + # WAR 1->2; WAW 1->3; expected=2 (('Eq(tu[t-1,x,y,z], tu[t,x,y,z] + tv[t,x,y,z])', 'Eq(tv[t,x,y,z], tu[t,x,y,z+2])', 'Eq(tu[t-2,x,y,0], tu[t,x,y,0] + 1.)'), - Forward, ['txyz', 'txyz', 'txy'], 'txyztxyztxy'), - # Time Backward, so flow dependences in time + '-***', ['txyz', 'txy'], 'txyz'), + # WAR 1->2; expected=1 + (('Eq(tu[t-1,x,y,z], tu[t,x,y,z] + tv[t,x,y,z])', + 'Eq(tv[t,x,y,z], tu[t,x,y,z+2] + tu[t,x,y,z-2])', + 'Eq(tw[t,x,y,z], tv[t,x,y,z] + 2)'), + '-***', ['txyz'], 'txyz'), + # Time goes backward so that information flows in time (('Eq(tu[t-1,x,y,z], tu[t,x+3,y,z] + tv[t,x,y,z])', 'Eq(tv[t-1,x,y,z], tu[t,x,y,z+2])', 'Eq(tw[t-1,x,y,z], tu[t,x,y+1,z] + tv[t,x,y-1,z])'), - Backward, ['txyz'], 'txyz'), - # Time Backward, so flow dependences in time, interleaved with independent Eq + '-***', ['txyz'], 'txyz'), + # Time goes backward so that information flows in time, interleaved + # with independent Eq (('Eq(tu[t-1,x,y,z], tu[t,x+3,y,z] + tv[t,x,y,z])', 'Eq(ti0[x,y,z], ti1[x,y,z+2])', 'Eq(tw[t-1,x,y,z], tu[t,x,y+1,z] + tv[t,x,y-1,z])'), - Backward, ['txyz', 'xyz'], 'txyzxyz'), - # Time Backward, so flow dependences in time, interleaved with dependent Eq + '-******', ['txyz', 'xyz'], 'txyzxyz'), + # Time goes backward so that information flows in time, interleaved + # with independent Eq (('Eq(ti0[x,y,z], ti1[x,y,z+2])', 'Eq(tu[t-1,x,y,z], tu[t,x+3,y,z] + tv[t,x,y,z])', 'Eq(tw[t-1,x,y,z], tu[t,x,y+1,z] + ti0[x,y-1,z])'), - Backward, ['xyz', 'txyz'], 'xyztxyz') + '*+*-*+*', ['xyz', 'txyz'], 'xyztxyz'), ]) - def test_consistency_anti_dependences(self, exprs, axis, expected, visit, + def test_consistency_anti_dependences(self, exprs, directions, expected, visit, ti0, ti1, ti3, tu, tv, tw): """ Test that anti dependences end up generating multi loop nests, rather @@ -875,15 +882,19 @@ def test_consistency_anti_dependences(self, exprs, axis, expected, visit, """ eq1, eq2, eq3 = EVAL(exprs, ti0.base, ti1.base, ti3.base, tu.base, tv.base, tw.base) - op = Operator([eq1, eq2, eq3], dse='noop', dle='noop', time_axis=axis) + op = Operator([eq1, eq2, eq3], dse='noop', dle='noop') trees = retrieve_iteration_tree(op) iters = FindNodes(Iteration).visit(op) assert len(trees) == len(expected) + assert len(iters) == len(directions) # mapper just makes it quicker to write out the test parametrization mapper = {'time': 't'} assert ["".join(mapper.get(i.dim.name, i.dim.name) for i in j) for j in trees] == expected assert "".join(mapper.get(i.dim.name, i.dim.name) for i in iters) == visit + # mapper just makes it quicker to write out the test parametrization + mapper = {'+': Forward, '-': Backward, '*': Any} + assert all(i.direction == mapper[j] for i, j in zip(iters, directions)) def test_expressions_imperfect_loops(self, ti0, ti1, ti2, t0): """ diff --git a/tests/test_save.py b/tests/test_save.py index 930f99143c..f67789b942 100644 --- a/tests/test_save.py +++ b/tests/test_save.py @@ -2,7 +2,7 @@ from sympy import solve from conftest import skipif_yask -from devito import Grid, Eq, Operator, TimeFunction, Forward +from devito import Grid, Eq, Operator, TimeFunction def initial(dx=0.01, dy=0.01): @@ -34,7 +34,7 @@ def run_simulation(save=False, dx=0.01, dy=0.01, a=0.5, timesteps=100): eqn = Eq(u.dt, a * (u.dx2 + u.dy2)) stencil = solve(eqn, u.forward)[0] - op = Operator(Eq(u.forward, stencil), time_axis=Forward) + op = Operator(Eq(u.forward, stencil)) op.apply(time=timesteps, dt=dt) if save: diff --git a/tests/test_timestepping.py b/tests/test_timestepping.py index 5ff795c0dd..e0e75f6bb6 100644 --- a/tests/test_timestepping.py +++ b/tests/test_timestepping.py @@ -3,7 +3,7 @@ import pytest from conftest import skipif_yask -from devito import Grid, Eq, Operator, Forward, Backward, TimeFunction +from devito import Grid, Eq, Operator, TimeFunction @pytest.fixture @@ -26,7 +26,7 @@ def b(grid): @pytest.fixture def c(grid): """Forward time data object, buffered (save=False)""" - return TimeFunction(name='c', grid=grid, time_order=1, save=None, time_axis=Forward) + return TimeFunction(name='c', grid=grid, time_order=1, save=None) @pytest.fixture @@ -46,7 +46,7 @@ def test_forward(a): @skipif_yask def test_backward(b): b.data[-1, :] = 7. - Operator(Eq(b.backward, b - 1.), time_axis=Backward)() + Operator(Eq(b.backward, b - 1.))() for i in range(b.shape[0]): assert np.allclose(b.data[i, :], 2. + i, rtol=1.e-12) @@ -69,10 +69,10 @@ def test_forward_backward(a, b, nt=5): a.data[0, :] = 1. b.data[0, :] = 1. eqn_a = Eq(a.forward, a + 1.) - Operator(eqn_a, time_axis=Forward)(time=nt) + Operator(eqn_a)(time=nt) eqn_b = Eq(b, a + 1.) - Operator(eqn_b, time_axis=Backward)(time=nt) + Operator(eqn_b)(time=nt) for i in range(nt): assert np.allclose(b.data[i, :], 2. + i, rtol=1.e-12) @@ -85,8 +85,8 @@ def test_forward_backward_overlapping(a, b, nt=5): """ a.data[0, :] = 1. b.data[0, :] = 1. - op_fwd = Operator(Eq(a.forward, a + 1.), time_axis=Forward) - op_bwd = Operator(Eq(b, a + 1.), time_axis=Backward) + op_fwd = Operator(Eq(a.forward, a + 1.)) + op_bwd = Operator(Eq(b, a + 1.)) op_fwd(time=nt) op_bwd(time=nt) @@ -99,7 +99,7 @@ def test_loop_bounds_forward(d): """Test the automatic bound detection for forward time loops""" d.data[:] = 1. eqn = Eq(d, 2. + d.dt2) - Operator(eqn, dle=None, dse=None, time_axis=Forward)(dt=1.) + Operator(eqn, dle=None, dse=None)(dt=1.) assert np.allclose(d.data[0, :], 1., rtol=1.e-12) assert np.allclose(d.data[-1, :], 1., rtol=1.e-12) for i in range(1, d.data.shape[0]-1): @@ -109,10 +109,11 @@ def test_loop_bounds_forward(d): @skipif_yask def test_loop_bounds_backward(d): """Test the automatic bound detection for backward time loops""" - d.data[:] = 1. - eqn = Eq(d, 2. + d.dt2) - Operator(eqn, dle=None, dse=None, time_axis=Backward)(dt=1.) - assert np.allclose(d.data[0, :], 1., rtol=1.e-12) - assert np.allclose(d.data[-1, :], 1., rtol=1.e-12) + d.data[:] = 5. + eqn = Eq(d.backward, d - 1) + op = Operator(eqn, dle=None, dse=None) + op(dt=1.) + assert np.allclose(d.data[0, :], 0., rtol=1.e-12) + assert np.allclose(d.data[-1, :], 5., rtol=1.e-12) for i in range(1, d.data.shape[0]-1): - assert np.allclose(d.data[i, :], d.data.shape[0] - i, rtol=1.e-12) + assert np.allclose(d.data[i, :], i, rtol=1.e-12) diff --git a/tests/test_tti.py b/tests/test_tti.py index bef9a40151..1e196bd556 100644 --- a/tests/test_tti.py +++ b/tests/test_tti.py @@ -82,14 +82,14 @@ def ricker_source(t, f0): time_order=2, space_order=space_order) # Create new wavefield object restart forward computation - u = TimeFunction(name='u', grid=model.grid, save=False, + u = TimeFunction(name='u', grid=model.grid, time_order=2, space_order=space_order, dtype=model.dtype) u.data[0:3, :] = u1.data[indlast, :] rec, _, _ = acoustic.forward(save=False, u=u) - utti = TimeFunction(name='u', grid=model.grid, save=False, + utti = TimeFunction(name='u', grid=model.grid, time_order=2, space_order=space_order, dtype=model.dtype) - vtti = TimeFunction(name='v', grid=model.grid, save=False, + vtti = TimeFunction(name='v', grid=model.grid, time_order=2, space_order=space_order, dtype=model.dtype) utti.data[0:3, :] = u1.data[indlast, :] vtti.data[0:3, :] = u1.data[indlast, :] diff --git a/tests/test_visitors.py b/tests/test_visitors.py index 9595839d45..6ce1faf5fe 100644 --- a/tests/test_visitors.py +++ b/tests/test_visitors.py @@ -154,20 +154,14 @@ def test_find_sections(exprs, block1, block2, block3): assert len(sections) == 2 found = list(sections.values()) assert len(found[0]) == 1 - assert found[0][0].stencil == exprs[0].stencil assert len(found[1]) == 1 - assert found[1][0].stencil == exprs[1].stencil sections = finder.visit(block3) assert len(sections) == 3 found = list(sections.values()) assert len(found[0]) == 1 - assert found[0][0].stencil == exprs[0].stencil assert len(found[1]) == 2 - assert found[1][0].stencil == exprs[1].stencil - assert found[1][1].stencil == exprs[2].stencil assert len(found[2]) == 1 - assert found[2][0].stencil == exprs[3].stencil @skipif_yask diff --git a/tests/test_yask.py b/tests/test_yask.py index 48eeb0cb41..7d2e6f2a6d 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -6,7 +6,7 @@ pexpect = pytest.importorskip('yask') # Run only if YASK is available from devito import (Eq, Grid, Operator, Constant, Function, TimeFunction, - SparseTimeFunction, Backward, configuration, clear_cache) # noqa + SparseTimeFunction, configuration, clear_cache) # noqa from devito.ir.iet import retrieve_iteration_tree # noqa from devito.yask.wrappers import contexts # noqa @@ -236,7 +236,7 @@ def test_reverse_time_loop(self): u = TimeFunction(name='yu4D', grid=grid, space_order=0, time_order=2) u.data[:] = 2. eq = Eq(u.backward, u - 1.) - op = Operator(eq, time_axis=Backward) + op = Operator(eq) op(yu4D=u, t=3) assert 'run_solution' in str(op) assert np.all(u.data[2] == 2.)