dsl: Add handling for missing subdomain specification in Eq

Author: EdCaunt

devito/finite_differences/differentiable.py

@@ -318,7 +318,7 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         """
         w = kwargs.get('weights', kwargs.get('w'))
         order = order or self.space_order
-        space_dims = [d.root for d in self.dimensions if d.is_Space]
+        space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
         derivs = tuple('d%s2' % d.name for d in space_dims)
         return Add(*[getattr(self, d)(x0=shift_x0(shift, space_dims[i], None, i),
@@ -345,7 +345,7 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
             Custom weights for the finite difference coefficients.
         """
         w = kwargs.get('weights', kwargs.get('w'))
-        space_dims = [d.root for d in self.dimensions if d.is_Space]
+        space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
         order = order or self.space_order
         return Add(*[getattr(self, 'd%s' % d.name)(x0=shift_x0(shift, d, None, i),
@@ -372,7 +372,7 @@ def grad(self, shift=None, order=None, method='FD', **kwargs):
             Custom weights for the finite
         """
         from devito.types.tensor import VectorFunction, VectorTimeFunction
-        space_dims = [d.root for d in self.dimensions if d.is_Space]
+        space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
         order = order or self.space_order
         w = kwargs.get('weights', kwargs.get('w'))
@@ -388,7 +388,7 @@ def biharmonic(self, weight=1):
         Generates a symbolic expression for the weighted biharmonic operator w.r.t.
         all spatial Dimensions Laplace(weight * Laplace (self))
         """
-        space_dims = [d.root for d in self.dimensions if d.is_Space]
+        space_dims = self.root_dimensions
         derivs = tuple('d%s2' % d.name for d in space_dims)
         return Add(*[getattr(self.laplace * weight, d) for d in derivs])
 

devito/ir/equations/algorithms.py

@@ -113,7 +113,7 @@ def lower_exprs(expressions, subs=None, **kwargs):
 def _lower_exprs(expressions, subs):
     processed = []
     for expr in as_tuple(expressions):
-        dimension_map = make_dimension_map(expr)
+        dimension_map = _make_dimension_map(expr)
 
         # Handle Functions (typical case)
         mapper = {f: _lower_exprs(f.indexify(subs=dimension_map), subs)
@@ -157,7 +157,7 @@ def _lower_exprs(expressions, subs):
         return processed.pop()
 
 
-def make_dimension_map(expr):
+def _make_dimension_map(expr):
     """
     Make the dimension_map for an expression. In the basic case, this is extracted
     directly from the SubDomain attached to the expression.
@@ -233,6 +233,7 @@ def _(expr, mapper, rebuilt, sregistry):
     thicknesses = {i for i in expr.free_symbols if isinstance(i, Thickness)}
     symbols = expr.free_symbols.difference(thicknesses)
 
+    # Iterate over all other symbols before iterating over standalone thicknesses
     for d in tuple(symbols) + tuple(thicknesses):
         _concretize_subdims(d, mapper, rebuilt, sregistry)
 

devito/operations/interpolators.py

@@ -33,7 +33,7 @@ def wrapper(interp, *args, **kwargs):
     return wrapper
 
 
-def extract_subdomain(variables):
+def _extract_subdomain(variables):
     """
     Check if any of the variables provided are defined on a SubDomain
     and extract it if this is the case.
@@ -194,6 +194,10 @@ def __init__(self, sfunction):
     def grid(self):
         return self.sfunction.grid
 
+    @property
+    def r(self):
+        return self.sfunction.r
+
     @memoized_meth
     def _weights(self, subdomain=None):
         raise NotImplementedError
@@ -202,10 +206,6 @@ def _weights(self, subdomain=None):
     def _gdims(self):
         return self.grid.dimensions
 
-    @property
-    def r(self):
-        return self.sfunction.r
-
     @cached_property
     def _cdim(self):
         """Base CustomDimensions used to construct _rdim"""
@@ -366,7 +366,7 @@ def _interpolate(self, expr, increment=False, self_subs={}, implicit_dims=None):
             _expr = expr
 
         variables = list(retrieve_function_carriers(_expr))
-        subdomain = extract_subdomain(variables)
+        subdomain = _extract_subdomain(variables)
 
         # Implicit dimensions
         implicit_dims = self._augment_implicit_dims(implicit_dims, variables)
@@ -416,7 +416,7 @@ def _inject(self, field, expr, implicit_dims=None):
             else:
                 assert len(exprs) == len(fields)
 
-        subdomain = extract_subdomain(fields)
+        subdomain = _extract_subdomain(fields)
 
         # Derivatives must be evaluated before the introduction of indirect accesses
         try:

devito/types/basic.py

@@ -1127,6 +1127,11 @@ def space_dimensions(self):
         """Tuple of Dimensions defining the physical space."""
         return tuple(d for d in self.dimensions if d.is_Space)
 
+    @cached_property
+    def root_dimensions(self):
+        """Tuple of root Dimensions of the physical space Dimensions."""
+        return tuple(d.root for d in self.space_dimensions)
+
     @property
     def base(self):
         return self.indexed
@@ -1461,7 +1466,7 @@ def _offset_subdomain(self):
         """Offset of subdomain indices versus the global index."""
         # If defined on a SubDomain, then need to offset indices accordingly
         if not self._is_on_subdomain:
-            return (0,)*len(self.dimensions)
+            return DimensionTuple(*[0 for _ in self.dimensions], getters=self.dimensions)
         # Symbolic offsets to avoid potential issues with user overrides
         offsets = []
         for d in self.dimensions:

devito/types/equation.py

@@ -1,6 +1,8 @@
 """User API to specify equations."""
 import sympy
 
+from functools import cached_property
+
 from devito.deprecations import deprecations
 from devito.tools import as_tuple, frozendict, Pickable
 from devito.types.lazy import Evaluable
@@ -140,10 +142,26 @@ def _flatten(self):
         else:
             return [self]
 
-    @property
+    @cached_property
     def subdomain(self):
         """The SubDomain in which the Eq is defined."""
-        return self._subdomain
+        if self._subdomain is not None:
+            return self._subdomain
+
+        from devito.symbolics.search import retrieve_functions  # noqa
+
+        funcs = retrieve_functions(self)
+        subdomains = {f.grid for f in funcs} - {None}
+        subdomains = {g for g in subdomains if g.is_SubDomain}
+
+        if len(subdomains) == 0:
+            return None
+        elif len(subdomains) == 1:
+            return subdomains.pop()
+        else:
+            raise ValueError("Multiple `SubDomain`s detected. Provide a `SubDomain`"
+                             " explicitly (i.e., via `Eq(..., subdomain=...)`)"
+                             " to unambiguously define the `Eq`'s iteration domain.")
 
     @property
     def substitutions(self):

devito/types/tensor.py

@@ -159,6 +159,11 @@ def space_dimensions(self):
         """Spatial dimensions."""
         return self._space_dimensions
 
+    @cached_property
+    def root_dimensions(self):
+        """Tuple of root Dimensions of the physical space Dimensions."""
+        return tuple(d.root for d in self.space_dimensions)
+
     @cached_property
     def space_order(self):
         """The space order for all components."""
@@ -209,7 +214,7 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         comps = []
         func = vec_func(self)
         ndim = len(self.space_dimensions)
-        space_dims = [d.root for d in self.space_dimensions]
+        space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (ndim, ndim))
         order = order or self.space_order
         for i in range(len(self.space_dimensions)):
@@ -252,7 +257,7 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         comps = []
         func = vec_func(self)
         order = order or self.space_order
-        space_dims = [d.root for d in self.space_dimensions]
+        space_dims = self.root_dimensions
         ndim = len(self.space_dimensions)
         shift_x0 = make_shift_x0(shift, (ndim, ndim))
         for j in range(ndim):
@@ -345,7 +350,7 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         w = kwargs.get('weights', kwargs.get('w'))
         shift_x0 = make_shift_x0(shift, (len(self.space_dimensions),))
         order = order or self.space_order
-        space_dims = [d.root for d in self.space_dimensions]
+        space_dims = self.root_dimensions
         return sum([getattr(self[i], 'd%s' % d.name)(x0=shift_x0(shift, d, None, i),
                                                      fd_order=order, method=method, w=w)
                     for i, d in enumerate(space_dims)])
@@ -378,7 +383,7 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         func = vec_func(self)
         shift_x0 = make_shift_x0(shift, (len(self.space_dimensions),))
         order = order or self.space_order
-        space_dims = [d.root for d in self.space_dimensions]
+        space_dims = self.root_dimensions
         comps = [sum([getattr(s, 'd%s2' % d.name)(x0=shift_x0(shift, d, None, i),
                                                   fd_order=order, w=w, method=method)
                       for i, d in enumerate(space_dims)])
@@ -406,7 +411,7 @@ def curl(self, shift=None, order=None, method='FD', **kwargs):
             raise AttributeError("Curl only supported for 3D VectorFunction")
         # The curl of a VectorFunction is a VectorFunction
         w = kwargs.get('weights', kwargs.get('w'))
-        dims = [d.root for d in self.space_dimensions]
+        dims = self.root_dimensions
         derivs = ['d%s' % d.name for d in dims]
         shift_x0 = make_shift_x0(shift, (len(dims), len(dims)))
         order = order or self.space_order
@@ -447,7 +452,7 @@ def grad(self, shift=None, order=None, method='FD', **kwargs):
         ndim = len(self.space_dimensions)
         shift_x0 = make_shift_x0(shift, (ndim, ndim))
         order = order or self.space_order
-        space_dims = [d.root for d in self.space_dimensions]
+        space_dims = self.root_dimensions
         comps = [[getattr(f, 'd%s' % d.name)(x0=shift_x0(shift, d, i, j), w=w,
                                              fd_order=order, method=method)
                   for j, d in enumerate(space_dims)]

examples/userapi/03_subdomains.ipynb

@@ -234,7 +234,7 @@
    "source": [
     "The Devito `SubDomain` API features two steps. The first of these is to create a subclass of `SubDomain`, which encapsulates the template by which `SubDomain`s should be constructed. This is not tied to a specific `Grid`. On a basic level, this specification is created by defining a suitable `define` method for the subclass.\n",
     "\n",
-    "Note if one intends to introduce a custom `__init__` when subclassing `SubDomain`, then it is necessary to provide the `kwarg` `grid=None`, and pass this through to `super().__init__()`. Such a custom `__init__` may be useful in applications such as creating many `SubDomain` instances with similar properties from a single template.\n",
+    "Note if one intends to introduce a custom `__init__` when subclassing `SubDomain`, then it is necessary for this `__init__` to take the `kwarg` `grid=None`, passing it through to `super().__init__(grid=grid)`. Such a custom `__init__` may be useful in applications such as creating many `SubDomain` instances with similar properties from a single template.\n",
     "\n",
     "Here is an example of how to define a subdomain that spans an entire 3D computational grid (this is similar to how 'domain' is defined internally):"
    ]
@@ -314,7 +314,7 @@
     "\n",
     "The two other options available are `'left'` and `'right'`. For a statement of the form `d: ('left', N)` the `SubDomain` spans a contiguous region of `N` points starting at `d`\\'s left extreme. A statement of the form `('right', N)` is analogous to the previous case but starting at the dimensions right extreme instead.\n",
     "\n",
-    "We can also attach this subdomain to the grid."
+    "We then attach this subdomain to the grid."
    ]
   },
   {