From 76c4bd6ba9a65d38724f688913cda35b32812df7 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 15:39:15 -0300
Subject: [PATCH 1/9] misc: Reduction in size of types/dimension

---
 devito/types/dimension.py | 183 +++++++++++---------------------------
 1 file changed, 50 insertions(+), 133 deletions(-)

diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index fec9fb78a3..945b4f9a6b 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -18,28 +18,24 @@
 from devito.types.relational import relational_max, relational_min
 
 __all__ = [
-    'BlockDimension',
-    'ConditionalDimension',
-    'CustomDimension',
-    'DefaultDimension',
-    'Dimension',
-    'IncrDimension',
-    'ModuloDimension',
-    'MultiSubDimension',
-    'SpaceDimension',
-    'Spacing',
-    'StencilDimension',
-    'SteppingDimension',
-    'SubDimension',
-    'TimeDimension',
-    'VirtualDimension',
-    'dimensions',
+    'BlockDimension', 'ConditionalDimension', 'CustomDimension', 'DefaultDimension',
+    'Dimension', 'IncrDimension', 'ModuloDimension', 'MultiSubDimension',
+    'SpaceDimension', 'Spacing', 'StencilDimension', 'SteppingDimension',
+    'SubDimension', 'TimeDimension', 'VirtualDimension', 'dimensions',
 ]
 
 
 SubDimensionThickness = namedtuple('SubDimensionThickness', 'left right')
 
 
+def _as_number(v):
+    """Return `v` as a SymPy Number if possible, otherwise `v` unchanged."""
+    try:
+        return sympy.Number(v)
+    except (TypeError, ValueError):
+        return v
+
+
 class Dimension(ArgProvider):
 
     """
@@ -359,14 +355,12 @@ def _arg_check(self, args, size, interval):
 
         # Allow the specific case of max=min-1, which disables the loop
         if args[self.max_name] < args[self.min_name]-1:
-            raise InvalidArgument(
-                f'Illegal {self.max_name}={args[self.max_name]} < '
-                f'{self.min_name}={args[self.min_name]}'
-            )
+            raise InvalidArgument(f'Illegal {self.max_name}={args[self.max_name]} < '
+                                  f'{self.min_name}={args[self.min_name]}')
         elif args[self.max_name] == args[self.min_name]-1:
             debug("%s=%d and %s=%d might cause no iterations along Dimension %s",
-                  self.min_name, args[self.min_name],
-                  self.max_name, args[self.max_name], self.name)
+                  self.min_name, args[self.min_name], self.max_name, args[self.max_name],
+                  self.name)
 
     # Pickling support
     __reduce_ex__ = Pickable.__reduce_ex__
@@ -546,7 +540,6 @@ def _arg_names(self):
 
     def _arg_check(self, *args, **kwargs):
         """A DerivedDimension performs no runtime checks."""
-        return
 
 
 # ***
@@ -627,7 +620,6 @@ class AbstractSubDimension(DerivedDimension):
     is_AbstractSub = True
 
     __rargs__ = DerivedDimension.__rargs__ + ('thickness',)
-    __rkwargs__ = ()
 
     _thickness_type = Thickness
 
@@ -737,10 +729,7 @@ class SubDimension(AbstractSubDimension):
 
     __rargs__ = AbstractSubDimension.__rargs__ + ('local',)
 
-    _thickness_type = Thickness
-
-    def __init_finalize__(self, name, parent, thickness, local,
-                          **kwargs):
+    def __init_finalize__(self, name, parent, thickness, local, **kwargs):
         self._local = local
         super().__init_finalize__(name, parent, thickness)
 
@@ -763,10 +752,9 @@ def _symbolic_thickness(self, thickness=None):
 
         names = [f"{self.parent.name}_{s}tkn" for s in ('l', 'r')]
         sides = [LEFT, RIGHT]
-        return SubDimensionThickness(*[
-            Thickness(name=n, side=s, value=t, **kwargs)
-            for n, s, t in zip(names, sides, thickness, strict=True)
-        ])
+        return SubDimensionThickness(*[Thickness(name=n, side=s, value=t, **kwargs)
+                                       for n, s, t in zip(names, sides, thickness,
+                                                          strict=True)])
 
     @cached_property
     def _interval(self):
@@ -1003,10 +991,7 @@ def uses_symbolic_factor(self):
     def factor_data(self):
         if isinstance(self.factor, Constant):
             return self.factor.data
-        elif self.factor is not None:
-            return self.factor
-        else:
-            return 1
+        return 1 if self.factor is None else self.factor
 
     @property
     def spacing(self):
@@ -1023,9 +1008,7 @@ def symbolic_factor(self):
         elif isinstance(self.factor, Constant):
             return self.factor
         else:
-            return SubsamplingFactor(
-                name=f'{self.name}f', dtype=np.int32, is_const=True
-            )
+            return SubsamplingFactor(name=f'{self.name}f', dtype=np.int32, is_const=True)
 
     @property
     def condition(self):
@@ -1178,14 +1161,9 @@ def origin(self):
 
     @cached_property
     def symbolic_size(self):
-        try:
+        with suppress(TypeError, ValueError):
             return sympy.Number(self.modulo)
-        except (TypeError, ValueError):
-            pass
-        try:
-            return sympy.Number(self.incr)
-        except (TypeError, ValueError):
-            return self.incr
+        return _as_number(self.incr)
 
     @cached_property
     def symbolic_min(self):
@@ -1193,10 +1171,7 @@ def symbolic_min(self):
             return self.offset % self.modulo
         # Make sure we return a symbolic object as this point `offset` may well
         # be a pure Python number
-        try:
-            return sympy.Number(self.offset)
-        except (TypeError, ValueError):
-            return self.offset
+        return _as_number(self.offset)
 
     @cached_property
     def symbolic_incr(self):
@@ -1205,10 +1180,7 @@ def symbolic_incr(self):
             incr = incr % self.modulo
         # Make sure we return a symbolic object as this point `incr` may well
         # be a pure Python number
-        try:
-            return sympy.Number(incr)
-        except (TypeError, ValueError):
-            return incr
+        return _as_number(incr)
 
     @cached_property
     def bound_symbols(self):
@@ -1225,21 +1197,17 @@ def _arg_values(self, *args, **kwargs):
     def __add__(self, other):
         # Exploit compatibility with addition:
         # `a1 ≡ b1 (mod n) and a2 ≡ b2 (mod n)` => `a1 + a2 ≡ b1 + b2 (mod n)`
-        try:
+        with suppress(AttributeError, TypeError, sympy.SympifyError):
             if self.modulo == other.modulo:
                 return self.origin + other.origin
-        except (AttributeError, TypeError, sympy.SympifyError):
-            pass
         return super().__add__(other)
 
     def __sub__(self, other):
         # Exploit compatibility with subtraction:
         # `a1 ≡ b1 (mod n) and a2 ≡ b2 (mod n)` => `a1 – a2 ≡ b1 – b2 (mod n)`
-        try:
+        with suppress(AttributeError, TypeError, sympy.SympifyError):
             if self.modulo == other.modulo:
                 return self.origin - other.origin
-        except (AttributeError, TypeError, sympy.SympifyError):
-            pass
         return super().__sub__(other)
 
 
@@ -1307,10 +1275,7 @@ def symbolic_size(self):
         if self.size is not None:
             # Make sure we return a symbolic object as the provided size might
             # be for example a pure int
-            try:
-                return sympy.Number(self.size)
-            except (TypeError, ValueError):
-                return self._size
+            return _as_number(self.size)
         else:
             # The size must be given as a function of the parent's symbols
             return self.symbolic_max - self.symbolic_min + 1
@@ -1319,26 +1284,17 @@ def symbolic_size(self):
     def symbolic_min(self):
         # Make sure we return a symbolic object as the provided min might
         # be for example a pure int
-        try:
-            return sympy.Number(self._min)
-        except (TypeError, ValueError):
-            return self._min
+        return _as_number(self._min)
 
     @cached_property
     def symbolic_max(self):
         # Make sure we return a symbolic object as the provided max might
         # be for example a pure int
-        try:
-            return sympy.Number(self._max)
-        except (TypeError, ValueError):
-            return self._max
+        return _as_number(self._max)
 
     @cached_property
     def symbolic_incr(self):
-        try:
-            return sympy.Number(self.step)
-        except (TypeError, ValueError):
-            return self.step
+        return _as_number(self.step)
 
     @cached_property
     def bound_symbols(self):
@@ -1379,11 +1335,7 @@ def _rebuild_hierarchy(self, callback=None, step=None):
 
         name0 = pp.name
 
-        if callback is None:
-            name1 = p.name
-        else:
-            base = callback(name0)
-            name1 = callback(f'{base}_blk')
+        name1 = p.name if callback is None else callback(f'{callback(name0)}_blk')
 
         bd = p._rebuild(name1, pp, step=step or p.step)
 
@@ -1391,10 +1343,8 @@ def _rebuild_hierarchy(self, callback=None, step=None):
         if step is not None:
             subs[p.step] = step
 
-        d = self._rebuild(
-            name0, bd,
-            self._min.subs(subs), self._max.subs(subs), size=self.size.subs(subs)
-        )
+        d = self._rebuild(name0, bd, self._min.subs(subs), self._max.subs(subs),
+                          size=self.size.subs(subs))
 
         return {self: d, p: bd}
 
@@ -1455,9 +1405,8 @@ def _arg_check(self, args, *_args):
                 )
         else:
             if value < 0:
-                raise InvalidArgument(
-                    f'Illegal block size `{name}={value}`: it should be > 0'
-                )
+                raise InvalidArgument(f'Illegal block size `{name}={value}`: '
+                                      'it should be > 0')
             if value > args[self.root.max_name] - args[self.root.min_name] + 1:
                 # Avoid OOB
                 raise InvalidArgument(
@@ -1482,8 +1431,7 @@ class CustomDimension(BasicDimension):
 
     is_Custom = True
 
-    __rkwargs__ = ('symbolic_min', 'symbolic_max', 'symbolic_size', 'parent',
-                   'local')
+    __rkwargs__ = ('symbolic_min', 'symbolic_max', 'symbolic_size', 'parent', 'local')
 
     def __init_finalize__(self, name, symbolic_min=None, symbolic_max=None,
                           symbolic_size=None, parent=None, local=True, **kwargs):
@@ -1512,17 +1460,11 @@ def index(self):
 
     @property
     def root(self):
-        if self.is_Derived:
-            return self.parent.root
-        else:
-            return self
+        return self.parent.root if self.is_Derived else self
 
     @property
     def spacing(self):
-        if self.is_Derived:
-            return self.parent.spacing
-        else:
-            return self._spacing
+        return self.parent.spacing if self.is_Derived else self._spacing
 
     @property
     def local(self):
@@ -1548,40 +1490,22 @@ def _defines(self):
 
     @cached_property
     def symbolic_min(self):
-        try:
-            return sympy.Number(self._symbolic_min)
-        except (TypeError, ValueError):
-            pass
         if self._symbolic_min is None:
             return super().symbolic_min
-        else:
-            return self._symbolic_min
+        return _as_number(self._symbolic_min)
 
     @cached_property
     def symbolic_max(self):
-        try:
-            return sympy.Number(self._symbolic_max)
-        except (TypeError, ValueError):
-            pass
         if self._symbolic_max is None:
             return super().symbolic_max
-        else:
-            return self._symbolic_max
+        return _as_number(self._symbolic_max)
 
     @cached_property
     def symbolic_size(self):
-        try:
-            return sympy.Number(self._symbolic_size)
-        except (TypeError, ValueError):
-            pass
         if self._symbolic_size is None:
             v = self.symbolic_max - self.symbolic_min + 1
-            if v.is_Number:
-                return v
-            else:
-                return super().symbolic_size
-        else:
-            return self._symbolic_size
+            return v if v.is_Number else super().symbolic_size
+        return _as_number(self._symbolic_size)
 
     def _arg_defaults(self, **kwargs):
         return {}
@@ -1591,7 +1515,6 @@ def _arg_values(self, *args, **kwargs):
 
     def _arg_check(self, *args, **kwargs):
         """A CustomDimension performs no runtime checks."""
-        return
 
 
 class DynamicDimensionMixin:
@@ -1647,18 +1570,13 @@ class StencilDimension(BasicDimension):
     __rargs__ = BasicDimension.__rargs__ + ('_min', '_max')
     __rkwargs__ = BasicDimension.__rkwargs__ + ('step',)
 
-    def __init_finalize__(self, name, _min, _max, spacing=1, step=1,
-                          **kwargs):
+    def __init_finalize__(self, name, _min, _max, spacing=1, step=1, **kwargs):
         self._spacing = sympy.sympify(spacing)
 
-        if not is_integer(_min):
-            raise ValueError(f"Expected integer `min` (got {_min})")
-        if not is_integer(_max):
-            raise ValueError(f"Expected integer `max` (got {_max})")
-        if not is_integer(self._spacing):
-            raise ValueError(f"Expected integer `spacing` (got {self._spacing})")
-        if not is_integer(step):
-            raise ValueError(f"Expected integer `step` (got {step})")
+        for k, v in (('min', _min), ('max', _max), ('spacing', self._spacing),
+                     ('step', step)):
+            if not is_integer(v):
+                raise ValueError(f"Expected integer `{k}` (got {v})")
 
         self._min = int(_min)
         self._max = int(_max)
@@ -1740,8 +1658,7 @@ class VirtualDimension(CustomDimension):
     __rkwargs__ = ('parent',)
 
     def __init_finalize__(self, name, parent=None):
-        super().__init_finalize__(name, parent=parent,
-                                  symbolic_min=sympy.S.Zero,
+        super().__init_finalize__(name, parent=parent, symbolic_min=sympy.S.Zero,
                                   symbolic_max=sympy.S.Zero)
 
 

From e03a32f2c80d4df3399827209283df252cec62f7 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 15:52:43 -0300
Subject: [PATCH 2/9] Further reduction of types

---
 devito/types/basic.py  | 122 ++++----------
 devito/types/sparse.py | 372 +++++++++++++----------------------------
 2 files changed, 151 insertions(+), 343 deletions(-)

diff --git a/devito/types/basic.py b/devito/types/basic.py
index 0effda1c86..699b644cca 100644
--- a/devito/types/basic.py
+++ b/devito/types/basic.py
@@ -24,14 +24,8 @@
 from devito.types.lazy import Evaluable
 from devito.types.utils import DimensionTuple, Offset, Size
 
-__all__ = [
-    'DeviceMap',
-    'Indexed',
-    'IndexedData',
-    'IrregularFunctionInterface',
-    'Scalar',
-    'Symbol',
-]
+__all__ = ['DeviceMap', 'Indexed', 'IndexedData', 'IrregularFunctionInterface',
+           'Scalar', 'Symbol']
 
 
 class CodeSymbol:
@@ -72,7 +66,6 @@ def dtype(self):
             * `numpy.dtype`: basic data types. For example, `np.float64 -> double`.
             * `ctypes`: composite objects (e.g., structs), foreign types.
         """
-        return
 
     @property
     @abc.abstractmethod
@@ -84,7 +77,6 @@ def _C_name(self):
         -------
         str
         """
-        return
 
     @property
     def _C_typedata(self):
@@ -102,13 +94,10 @@ def _C_typedata(self):
         if _type is c_char_p:
             _type = c_char
 
-        try:
-            # We have internal types such as c_complex that are
-            # Structure too but should be treated as plain c_type
-            _ = _type._base_dtype
-        except AttributeError:
-            if issubclass(_type, Structure):
-                _type = f'struct {_type.__name__}'
+        # We have internal types such as c_complex that are
+        # Structure too but should be treated as plain c_type
+        if not hasattr(_type, '_base_dtype') and issubclass(_type, Structure):
+            _type = f'struct {_type.__name__}'
 
         return _type
 
@@ -118,7 +107,6 @@ def _C_ctype(self):
         """
         The type of the object in the generated code as a `ctypes` class.
         """
-        return
 
     @property
     def _C_symbol(self):
@@ -380,15 +368,13 @@ class AbstractSymbol(sympy.Symbol, Basic, Pickable, Evaluable):
     @classmethod
     def _filter_assumptions(cls, **kwargs):
         """Extract sympy.Symbol-specific kwargs."""
-        assumptions = {}
         # Pop predefined assumptions
         for key in ('real', 'imaginary', 'commutative'):
             kwargs.pop(key, None)
 
         # Extract sympy.Symbol-specific kwargs
-        for i in list(kwargs):
-            if i in _assume_rules.defined_facts:
-                assumptions[i] = kwargs.pop(i)
+        assumptions = {i: kwargs.pop(i) for i in list(kwargs)
+                       if i in _assume_rules.defined_facts}
 
         return assumptions, kwargs
 
@@ -509,11 +495,9 @@ def _subs(self, old, new, **hints):
         involving devito Scalars.  Ordinarily the comparisons between
         devito subclasses of sympy types are quite strict.
         """
-        try:
+        with suppress(AttributeError):
             if old.is_Symbol and old.name == self.name:
                 return new
-        except AttributeError:
-            pass
 
         return self
 
@@ -635,14 +619,12 @@ def _arg_defaults(self, **kwargs):
             # through a wrapper object (e.g., a Dimension spacing `h_x` gets its
             # value via a Grid object)
             return {}
-        else:
-            return {self.name: self.default_value}
+        return {self.name: self.default_value}
 
     def _arg_values(self, **kwargs):
         if self.name in kwargs:
             return {self.name: kwargs.pop(self.name)}
-        else:
-            return self._arg_defaults(**kwargs)
+        return self._arg_defaults(**kwargs)
 
 
 class AbstractFunction(sympy.Function, Basic, Pickable, Evaluable):
@@ -744,9 +726,8 @@ def __new__(cls, *args, **kwargs):
 
         # If same name/indices and `function` isn't None, then it's
         # definitely a reconstruction
-        if function is not None and \
-           function.name == name and \
-           function.indices == indices:
+        if (function is not None and function.name == name and
+                function.indices == indices):
             # Special case: a syntactically identical alias of `function`, so
             # let's just return `function` itself
             return function
@@ -863,10 +844,8 @@ def __init_finalize__(self, *args, **kwargs):
         # Averaging mode for off the grid evaluation
         self._avg_mode = kwargs.get('avg_mode', 'arithmetic')
         if self._avg_mode not in ['arithmetic', 'harmonic', 'safe_harmonic']:
-            raise ValueError(
-                f"Invalid averaging mode_mode {self._avg_mode}, accepted values are"
-                " arithmetic or harmonic"
-            )
+            raise ValueError(f"Invalid averaging mode_mode {self._avg_mode}, accepted "
+                             "values are arithmetic or harmonic")
 
     @classmethod
     def __args_setup__(cls, *args, **kwargs):
@@ -972,14 +951,9 @@ def origin(self):
         f(x) : origin = 0
         f(x + hx/2) : origin = hx/2
         """
-        return DimensionTuple(*(
-            r - d + o
-            for d, r, o in zip(
-                self.dimensions,
-                self.indices_ref,
-                self._offset_subdomain, strict=True
-            )
-        ), getters=self.dimensions)
+        return DimensionTuple(*(r - d + o for d, r, o in zip(
+            self.dimensions, self.indices_ref, self._offset_subdomain, strict=True
+        )), getters=self.dimensions)
 
     @property
     def dimensions(self):
@@ -993,10 +967,7 @@ def _padded_dimensions(self):
         except IndexError:
             return ()
 
-        if d is self.dimensions[-1]:
-            return (d,)
-        else:
-            return ()
+        return (d,) if d is self.dimensions[-1] else ()
 
     @cached_property
     def space_dimensions(self):
@@ -1192,12 +1163,9 @@ def dmap(self):
         return None.
         """
         if self._mem_mapped:
-            return DeviceMap(f'd_{self.name}', shape=self._shape,
-                             function=self.function)
+            return DeviceMap(f'd_{self.name}', shape=self._shape, function=self.function)
         elif self._mem_local:
             return self.indexed
-        else:
-            return None
 
     @property
     def size(self):
@@ -1300,23 +1268,23 @@ def _dist_dimensions(self):
         """The Dimensions decomposed for distributed-parallelism."""
         if self._distributor is None:
             return ()
-        else:
-            return tuple(d for d in self.dimensions if d in self._distributor.dimensions)
+        return tuple(d for d in self.dimensions if d in self._distributor.dimensions)
 
     @cached_property
     def _size_domain(self):
         """Number of points in the domain region."""
         return DimensionTuple(*self.shape, getters=self.dimensions)
 
+    def _make_sizes(self, pairs):
+        left = tuple(zip(*pairs, strict=True))[0]
+        right = tuple(zip(*pairs, strict=True))[1]
+        sizes = tuple(Size(i, j) for i, j in pairs)
+        return DimensionTuple(*sizes, getters=self.dimensions, left=left, right=right)
+
     @cached_property
     def _size_halo(self):
         """Number of points in the halo region."""
-        left = tuple(zip(*self._halo, strict=True))[0]
-        right = tuple(zip(*self._halo, strict=True))[1]
-
-        sizes = tuple(Size(i, j) for i, j in self._halo)
-
-        return DimensionTuple(*sizes, getters=self.dimensions, left=left, right=right)
+        return self._make_sizes(self._halo)
 
     @cached_property
     def _size_owned(self):
@@ -1331,20 +1299,13 @@ def _size_owned(self):
     @cached_property
     def _size_padding(self):
         """Number of points in the padding region."""
-        left = tuple(zip(*self._padding, strict=True))[0]
-        right = tuple(zip(*self._padding, strict=True))[1]
-
-        sizes = tuple(Size(i, j) for i, j in self._padding)
-
-        return DimensionTuple(*sizes, getters=self.dimensions, left=left, right=right)
+        return self._make_sizes(self._padding)
 
     @cached_property
     def _size_nopad(self):
         """Number of points in the domain+halo region."""
-        sizes = tuple(
-            i+sum(j)
-            for i, j in zip(self._size_domain, self._size_halo, strict=True)
-        )
+        sizes = tuple(i+sum(j) for i, j in zip(self._size_domain, self._size_halo,
+                                               strict=True))
         return DimensionTuple(*sizes, getters=self.dimensions)
 
     @cached_property
@@ -1756,10 +1717,8 @@ def _eval_matrix_mul(self, other):
                 row, col = i // other.cols, i % other.cols
                 row_indices = range(self_cols*row, self_cols*(row+1))
                 col_indices = range(col, other_len, other.cols)
-                vec = [
-                    mat[a]*other_mat[b]
-                    for a, b in zip(row_indices, col_indices, strict=True)
-                ]
+                vec = [mat[a]*other_mat[b]
+                       for a, b in zip(row_indices, col_indices, strict=True)]
                 new_mat[i] = sum(vec)
 
         # Get new class and return product
@@ -1801,10 +1760,7 @@ def sort_key(self, order=None):
     def __getitem__(self, indices, **kwargs):
         """Produce a types.Indexed, rather than a sympy.Indexed."""
         # Is there a specific Indexed class to use?
-        try:
-            cls = self.function._indexed_cls
-        except AttributeError:
-            cls = Indexed
+        cls = getattr(self.function, '_indexed_cls', Indexed)
 
         return cls(self, *as_tuple(indices))
 
@@ -1977,14 +1933,12 @@ def _subs(self, old, new, **hints):
         # Wrap in a try to make sure no substitution happens when
         # old is an Indexed as only checkink `old is new` would lead to
         # incorrect substitution of `old.base` by `new`
-        try:
+        with suppress(AttributeError):
             if old.is_Indexed:
                 if old.base == self.base and old.indices == self.indices:
                     return new
                 else:
                     return self
-        except AttributeError:
-            pass
         return super()._subs(old, new, **hints)
 
     def _translate(self, mapper=None):
@@ -2002,10 +1956,8 @@ def _translate(self, mapper=None):
         mapper = mapper or {self.dimensions[-1]: 1}
 
         if any(d not in mapper for d in self.dimensions):
-            raise ValueError(
-                f"Cannot translate {self} with mapper {mapper} since not "
-                "all dimensions are covered"
-            )
+            raise ValueError(f"Cannot translate {self} with mapper {mapper} since not "
+                             "all dimensions are covered")
 
         translations = [mapper.get(d, 0) for d in self.dimensions]
         indices = [sum(i) for i in zip(self.indices, translations, strict=True)]
diff --git a/devito/types/sparse.py b/devito/types/sparse.py
index 6e913646a3..fd81d53a37 100644
--- a/devito/types/sparse.py
+++ b/devito/types/sparse.py
@@ -25,13 +25,8 @@
 from devito.types.equation import Eq, Inc
 from devito.types.utils import DimensionTuple, IgnoreDimSort
 
-__all__ = [
-    'MatrixSparseTimeFunction',
-    'PrecomputedSparseFunction',
-    'PrecomputedSparseTimeFunction',
-    'SparseFunction',
-    'SparseTimeFunction',
-]
+__all__ = ['MatrixSparseTimeFunction', 'PrecomputedSparseFunction',
+           'PrecomputedSparseTimeFunction', 'SparseFunction', 'SparseTimeFunction']
 
 
 _interpolators = {'linear': LinearInterpolator, 'sinc': SincInterpolator,
@@ -96,8 +91,7 @@ def __indices_setup__(cls, *args, **kwargs):
 
         if args:
             return tuple(dimensions), tuple(args)
-        else:
-            return dimensions, dimensions
+        return dimensions, dimensions
 
     @classmethod
     def __shape_setup__(cls, **kwargs):
@@ -183,10 +177,7 @@ def __subfunc_setup__(self, suffix, keys, dtype=None, inkwargs=False, **kwargs):
             d = self.indices[self._sparse_position]
             if d in key.indices:
                 # Can use as is, dimension already matches
-                if self.alias:
-                    return key._rebuild(alias=self.alias, name=name)
-                else:
-                    return key
+                return key._rebuild(alias=self.alias, name=name) if self.alias else key
             else:
                 # Need to rebuild so the dimensions match the parent
                 # SparseFunction, for example we end up here via `.subs(d, new_d)`
@@ -209,19 +200,14 @@ def __subfunc_setup__(self, suffix, keys, dtype=None, inkwargs=False, **kwargs):
             # Fallback to default behaviour
             dtype = dtype or self.dtype
         else:
-            if shape != key.shape and \
-               key.shape != (shape[1],) and \
-               self._distributor.nprocs == 1:
-                raise ValueError(
-                    f'Incompatible shape for {suffix}, `{key.shape[:2]}`;'
-                    f'expected `{shape}`'
-                )
+            if (shape != key.shape and key.shape != (shape[1],) and
+                    self._distributor.nprocs == 1):
+                raise ValueError(f'Incompatible shape for {suffix}, `{key.shape[:2]}`;'
+                                 f'expected `{shape}`')
 
             # Infer dtype
-            if np.issubdtype(key.dtype.type, np.integer):
-                dtype = dtype or np.int32
-            else:
-                dtype = dtype or self.dtype
+            dtype = dtype or (np.int32 if np.issubdtype(key.dtype.type, np.integer)
+                              else self.dtype)
 
         # Whether to initialize the subfunction with the provided data
         # Useful when rebuilding with a placeholder array only used to
@@ -231,11 +217,9 @@ def __subfunc_setup__(self, suffix, keys, dtype=None, inkwargs=False, **kwargs):
         # Complex coordinates are not valid, so fall back to corresponding
         # real floating point type if dtype is complex.
         dtype = dtype(0).real.__class__
-        sf = SparseSubFunction(
-            name=name, dtype=dtype, dimensions=dimensions,
-            shape=shape, space_order=0, alias=self.alias,
-            distributor=self._distributor, parent=self, **init
-        )
+        sf = SparseSubFunction(name=name, dtype=dtype, dimensions=dimensions, shape=shape,
+                               space_order=0, alias=self.alias,
+                               distributor=self._distributor, parent=self, **init)
 
         if self.npoint == 0:
             # This is a corner case -- we might get here, for example, when
@@ -270,8 +254,7 @@ def _grid_map(self):
 
     @cached_property
     def origin(self):
-        return DimensionTuple(*[0]*len(self.dimensions),
-                              getters=self.dimensions)
+        return DimensionTuple(*[0]*len(self.dimensions), getters=self.dimensions)
 
     @property
     def _mpitype(self):
@@ -291,13 +274,12 @@ def _comm(self):
     def _coords_indices(self):
         if self.gridpoints_data is not None:
             return self.gridpoints_data
-        else:
-            if self.coordinates_data is None:
-                raise ValueError("No coordinates or gridpoints attached"
-                                 "to this SparseFunction")
-            return (
-                np.floor((self.coordinates_data - self.grid.origin) / self.grid.spacing)
-            ).astype(int)
+        if self.coordinates_data is None:
+            raise ValueError("No coordinates or gridpoints attached"
+                             "to this SparseFunction")
+        return (
+            np.floor((self.coordinates_data - self.grid.origin) / self.grid.spacing)
+        ).astype(int)
 
     @property
     def _support(self):
@@ -399,17 +381,10 @@ def _position_map(self, shifts=None):
         is None, only the grid origin is subtracted.
         """
         shifts = shifts or (0,) * len(self.grid.dimensions)
-        return OrderedDict([
-            ((c - o - s)/d.spacing, p)
-            for p, c, d, o, s in zip(
-                self._pos_symbols(shifts=shifts),
-                self._coordinate_symbols,
-                self.grid.dimensions,
-                self.grid.origin_symbols,
-                shifts,
-                strict=True
-            )
-        ])
+        return OrderedDict([((c - o - s)/d.spacing, p) for p, c, d, o, s in zip(
+            self._pos_symbols(shifts=shifts), self._coordinate_symbols,
+            self.grid.dimensions, self.grid.origin_symbols, shifts, strict=True
+        )])
 
     @cached_property
     def dist_origin(self):
@@ -473,18 +448,13 @@ def guard(self, expr=None):
         temps = self.interpolator._positions(self.dimensions)
 
         # Create positions and indices temporaries/indirections
-        for d, pos in zip(
-            self.grid.dimensions,
-            pmap.values(),
-            strict=True
-        ):
+        for d, pos in zip(self.grid.dimensions, pmap.values(), strict=True):
             # Add conditional to avoid OOB
             lb = sympy.And(pos >= d.symbolic_min, evaluate=False)
             ub = sympy.And(pos <= d.symbolic_max, evaluate=False)
             conditions[d] = sympy.And(lb, ub, evaluate=False)
         condition = sympy.And(*conditions.values(), evaluate=False)
-        cd = ConditionalDimension(self._sparse_dim.name,
-                                  self._sparse_dim,
+        cd = ConditionalDimension(self._sparse_dim.name, self._sparse_dim,
                                   condition=condition, indirect=True)
 
         if expr is None:
@@ -769,9 +739,8 @@ def _arg_values(self, estimate_memory=False, **kwargs):
                 # Pure-data replacement (ndarray). Re-derive full defaults so
                 # any interpolator-owned SubFunctions get rebuilt alongside
                 # the scattered data.
-                values = self._arg_defaults(
-                    alias=self, estimate_memory=estimate_memory
-                ).reduce_all()
+                values = self._arg_defaults(alias=self,
+                                            estimate_memory=estimate_memory).reduce_all()
                 for k, v in self._dist_scatter(data=new).items():
                     values[k.name] = v
                     for i, s in zip(k.indices, v.shape, strict=True):
@@ -852,8 +821,7 @@ def __indices_setup__(cls, *args, **kwargs):
 
         if args:
             return tuple(dimensions), tuple(args)
-        else:
-            return dimensions, dimensions
+        return dimensions, dimensions
 
     @property
     def nt(self):
@@ -995,8 +963,7 @@ def __interp_setup__(self, interpolation='linear', r=None, **kwargs):
     def _coordinate_symbols(self):
         """Symbol representing the coordinate values in each Dimension."""
         d_dim = self.coordinates.dimensions[1]
-        return tuple([self.coordinates._subs(d_dim, i)
-                      for i in range(self.grid.dim)])
+        return tuple([self.coordinates._subs(d_dim, i) for i in range(self.grid.dim)])
 
     @cached_property
     def _decomposition(self):
@@ -1009,8 +976,7 @@ def _arg_defaults(self, alias=None, estimate_memory=False):
             return defaults
         key = alias or self
         coords = defaults.get(key.coordinates.name, self.coordinates.data)
-        defaults.update(key.interpolator._arg_defaults(coords=coords,
-                                                       sfunc=self))
+        defaults.update(key.interpolator._arg_defaults(coords=coords, sfunc=self))
         return defaults
 
     def _arg_values(self, estimate_memory=False, **kwargs):
@@ -1029,9 +995,8 @@ def _arg_values(self, estimate_memory=False, **kwargs):
         origin = tuple(kwargs.get(n, o) for n, o in
                        zip(onames, self.grid.origin, strict=True))
         coords = values.get(self.coordinates.name, self.coordinates.data)
-        values.update(self.interpolator._arg_defaults(
-            coords=coords, sfunc=self, origin=origin
-        ))
+        values.update(self.interpolator._arg_defaults(coords=coords, sfunc=self,
+                                                      origin=origin))
 
         return values
 
@@ -1247,8 +1212,7 @@ class PrecomputedSparseFunction(AbstractSparseFunction):
     _sub_functions = ('gridpoints', 'coordinates', 'interpolation_coeffs')
 
     __rkwargs__ = (AbstractSparseFunction.__rkwargs__ +
-                   ('r', 'gridpoints', 'coordinates',
-                    'interpolation_coeffs'))
+                   ('r', 'gridpoints', 'coordinates', 'interpolation_coeffs'))
 
     def __init_finalize__(self, *args, **kwargs):
         super().__init_finalize__(*args, **kwargs)
@@ -1261,14 +1225,12 @@ def __init_finalize__(self, *args, **kwargs):
         # Subfunctions setup
         self._dist_origin = {}
         dtype = kwargs.pop('dtype', self.grid.dtype)
-        self._gridpoints = self.__subfunc_setup__('gridpoints',
-                                                  ('gridpoints', 'gridpoints_data'),
-                                                  inkwargs=True,
-                                                  dtype=np.int32, **kwargs)
-        self._coordinates = self.__subfunc_setup__('coords',
-                                                   ('coordinates', 'coordinates_data'),
-                                                   inkwargs=self._gridpoints is not None,
-                                                   dtype=dtype, **kwargs)
+        self._gridpoints = self.__subfunc_setup__(
+            'gridpoints', ('gridpoints', 'gridpoints_data'), inkwargs=True,
+            dtype=np.int32, **kwargs)
+        self._coordinates = self.__subfunc_setup__(
+            'coords', ('coordinates', 'coordinates_data'),
+            inkwargs=self._gridpoints is not None, dtype=dtype, **kwargs)
 
         if self._coordinates is not None:
             self._dist_origin.update({self._coordinates: self.grid.origin_offset})
@@ -1293,10 +1255,8 @@ def __init_finalize__(self, *args, **kwargs):
                 if nr == r:
                     r = r // 2
                 else:
-                    raise ValueError(
-                        f'Interpolation coefficients shape {r} do not match'
-                        f'specified radius {nr}'
-                    )
+                    raise ValueError(f'Interpolation coefficients shape {r} do not match'
+                                     f'specified radius {nr}')
         self._radius = r
         self._dist_origin.update({self._interpolation_coeffs: None})
 
@@ -1316,18 +1276,11 @@ def _coordinate_symbols(self):
         """Symbol representing the coordinate values in each Dimension."""
         if self.gridpoints is not None:
             d_dim = self.gridpoints.dimensions[1]
-            return tuple([
-                self.gridpoints._subs(d_dim, di) * d.spacing + o
-                for ((di, d), o) in zip(
-                    enumerate(self.grid.dimensions),
-                    self.grid.origin,
-                    strict=True
-                )
-            ])
-        else:
-            d_dim = self.coordinates.dimensions[1]
-            return tuple([self.coordinates._subs(d_dim, i)
-                          for i in range(self.grid.dim)])
+            return tuple([self.gridpoints._subs(d_dim, di) * d.spacing + o
+                          for ((di, d), o) in zip(enumerate(self.grid.dimensions),
+                                                  self.grid.origin, strict=True)])
+        d_dim = self.coordinates.dimensions[1]
+        return tuple([self.coordinates._subs(d_dim, i) for i in range(self.grid.dim)])
 
     @memoized_meth
     def _position_map(self, shifts=None):
@@ -1346,16 +1299,10 @@ def _position_map(self, shifts=None):
         """
         if self.gridpoints_data is not None:
             ddim = self.gridpoints.dimensions[-1]
-            return OrderedDict(
-                (self.gridpoints._subs(ddim, di), p)
-                for (di, p) in zip(
-                    range(self.grid.dim),
-                    self._pos_symbols(shifts=shifts),
-                    strict=True
-                )
-            )
-        else:
-            return super()._position_map(shifts=shifts)
+            return OrderedDict((self.gridpoints._subs(ddim, di), p) for (di, p) in zip(
+                range(self.grid.dim), self._pos_symbols(shifts=shifts), strict=True
+            ))
+        return super()._position_map(shifts=shifts)
 
 
 class PrecomputedSparseTimeFunction(AbstractSparseTimeFunction,
@@ -1578,12 +1525,9 @@ def __init_finalize__(self, *args, **kwargs):
         locdim = Dimension(f'loc_{self.name}')
 
         self._gridpoints = SubFunction(
-            name=f"{self.name}_gridpoints",
-            dtype=np.int32,
-            dimensions=(locdim, ddim),
-            shape=(nloc, self.grid.dim),
-            allocator=self._allocator,
-            space_order=0, parent=self)
+            name=f"{self.name}_gridpoints", dtype=np.int32, dimensions=(locdim, ddim),
+            shape=(nloc, self.grid.dim), allocator=self._allocator, space_order=0,
+            parent=self)
 
         # There is a coefficient array per grid Dimension
         # I could pack these into one array but that seems less readable?
@@ -1592,10 +1536,8 @@ def __init_finalize__(self, *args, **kwargs):
         self.rdims = []
         for d in self.grid.dimensions:
             if self._radius[d] is not None:
-                rdim = DefaultDimension(
-                    name=f'r{d.name}_{self.name}',
-                    default_value=self._radius[d]
-                )
+                rdim = DefaultDimension(name=f'r{d.name}_{self.name}',
+                                        default_value=self._radius[d])
                 self.rdims.append(rdim)
                 coeff_dim = rdim
                 coeff_shape = self._radius[d]
@@ -1604,18 +1546,13 @@ def __init_finalize__(self, *args, **kwargs):
                 coeff_shape = self.grid.size_map[d].glb
 
             self.interpolation_coefficients[d] = SubFunction(
-                name=f"{self.name}_coefficients_{d.name}",
-                dtype=self.dtype,
-                dimensions=(locdim, coeff_dim),
-                shape=(nloc, coeff_shape),
-                allocator=self._allocator,
-                space_order=0, parent=self)
+                name=f"{self.name}_coefficients_{d.name}", dtype=self.dtype,
+                dimensions=(locdim, coeff_dim), shape=(nloc, coeff_shape),
+                allocator=self._allocator, space_order=0, parent=self)
 
             # For the _sub_functions, these must be named attributes of
             # this SparseFunction object
-            setattr(
-                self, f"coefficients_{d.name}",
-                self.interpolation_coefficients[d])
+            setattr(self, f"coefficients_{d.name}", self.interpolation_coefficients[d])
 
         # We also need arrays to represent the sparse matrix map
         # The shapes are bogus; these are really only used when
@@ -1632,32 +1569,14 @@ def __init_finalize__(self, *args, **kwargs):
             nnz_size = 1
 
         self._mrow = DynamicSubFunction(
-            name=f'mrow_{self.name}',
-            dtype=np.int32,
-            dimensions=(self.nnzdim,),
-            shape=(nnz_size,),
-            space_order=0,
-            parent=self,
-            allocator=self._allocator,
-        )
+            name=f'mrow_{self.name}', dtype=np.int32, dimensions=(self.nnzdim,),
+            shape=(nnz_size,), space_order=0, parent=self, allocator=self._allocator)
         self._mcol = DynamicSubFunction(
-            name=f'mcol_{self.name}',
-            dtype=np.int32,
-            dimensions=(self.nnzdim,),
-            shape=(nnz_size,),
-            space_order=0,
-            parent=self,
-            allocator=self._allocator,
-        )
+            name=f'mcol_{self.name}', dtype=np.int32, dimensions=(self.nnzdim,),
+            shape=(nnz_size,), space_order=0, parent=self, allocator=self._allocator)
         self._mval = DynamicSubFunction(
-            name=f'mval_{self.name}',
-            dtype=self.dtype,
-            dimensions=(self.nnzdim,),
-            shape=(nnz_size,),
-            space_order=0,
-            parent=self,
-            allocator=self._allocator,
-        )
+            name=f'mval_{self.name}', dtype=self.dtype, dimensions=(self.nnzdim,),
+            shape=(nnz_size,), space_order=0, parent=self, allocator=self._allocator)
 
         # This loop maintains a map of nnz indices which touch each
         # coordinate of the parallelised injection Dimension
@@ -1668,32 +1587,20 @@ def __init_finalize__(self, *args, **kwargs):
         # This map acts as an indirect sort of the sources according to their
         # position along the parallelisation dimension
         self._par_dim_to_nnz_map = DynamicSubFunction(
-            name=f'par_dim_to_nnz_map_{self.name}',
-            dtype=np.int32,
+            name=f'par_dim_to_nnz_map_{self.name}', dtype=np.int32,
             dimensions=(self.par_dim_to_nnz_dim,),
             # shape is unknown at this stage
-            shape=(1,),
-            space_order=0,
-            parent=self,
-        )
+            shape=(1,), space_order=0, parent=self)
         self._par_dim_to_nnz_m = DynamicSubFunction(
-            name=f'par_dim_to_nnz_m_{self.name}',
-            dtype=np.int32,
+            name=f'par_dim_to_nnz_m_{self.name}', dtype=np.int32,
             dimensions=(self._par_dim,),
             # shape is unknown at this stage
-            shape=(1,),
-            space_order=0,
-            parent=self,
-        )
+            shape=(1,), space_order=0, parent=self)
         self._par_dim_to_nnz_M = DynamicSubFunction(
-            name=f'par_dim_to_nnz_M_{self.name}',
-            dtype=np.int32,
+            name=f'par_dim_to_nnz_M_{self.name}', dtype=np.int32,
             dimensions=(self._par_dim,),
             # shape is unknown at this stage
-            shape=(1,),
-            space_order=0,
-            parent=self,
-        )
+            shape=(1,), space_order=0, parent=self)
 
         if self._distributor.nprocs == 1:
             self._mrow.data[:] = m_coo.row
@@ -1904,26 +1811,15 @@ def inject(self, field, expr, u_t=None, p_t=None):
 
         rhs = prod(coeffs) * expr
         field = field.subs(dim_subs)
-        out = [
-            Eq(
-                par_dim_to_nnz_dim.symbolic_min,
-                self._par_dim_to_nnz_m,
-                implicit_dims=tuple(implicit_dims_for_range)
-            ),
-            Eq(
-                par_dim_to_nnz_dim.symbolic_max,
-                self._par_dim_to_nnz_M,
-                implicit_dims=tuple(implicit_dims_for_range)
-            ),
-            Inc(
-                field,
-                rhs.subs(dim_subs),
-                implicit_dims=IgnoreDimSort(implicit_dims_for_inject),
-            ),
+        return [
+            Eq(par_dim_to_nnz_dim.symbolic_min, self._par_dim_to_nnz_m,
+               implicit_dims=tuple(implicit_dims_for_range)),
+            Eq(par_dim_to_nnz_dim.symbolic_max, self._par_dim_to_nnz_M,
+               implicit_dims=tuple(implicit_dims_for_range)),
+            Inc(field, rhs.subs(dim_subs),
+                implicit_dims=IgnoreDimSort(implicit_dims_for_inject)),
         ]
 
-        return out
-
     @classmethod
     def __shape_setup__(cls, **kwargs):
         # This happens before __init__, so we have to get 'npoint'
@@ -1991,11 +1887,9 @@ def _rank_to_points(self):
                 dim_r = self.grid.size_map[dim].glb
 
             # Define the split
-            dim_breaks[:-2:2] = [
-                decomp_part[0] - self.r + 1 for decomp_part in decomp]
+            dim_breaks[:-2:2] = [decomp_part[0] - self.r + 1 for decomp_part in decomp]
             dim_breaks[-2] = decomp[-1][-1] + 1 - self.r + 1
-            dim_breaks[1:-1:2] = [
-                decomp_part[0] for decomp_part in decomp]
+            dim_breaks[1:-1:2] = [decomp_part[0] for decomp_part in decomp]
             dim_breaks[-1] = decomp[-1][-1] + 1
 
             # Handle the radius is None case by ensuring we treat
@@ -2005,12 +1899,10 @@ def _rank_to_points(self):
                 gridpoints_dim = np.zeros_like(gridpoints_dim)
 
             try:
-                binned_gridpoints[:, idim] = np.digitize(
-                    gridpoints_dim, dim_breaks)
+                binned_gridpoints[:, idim] = np.digitize(gridpoints_dim, dim_breaks)
             except ValueError as e:
-                raise ValueError(
-                    "decomposition failed!  Are some ranks too skinny?"
-                ) from e
+                raise ValueError("decomposition failed!  Are some ranks too skinny?") \
+                    from e
 
             this_group_rank_map = {
                 0: {None},
@@ -2025,21 +1917,16 @@ def _rank_to_points(self):
         # This allows the points to be grouped into non-overlapping sets
         # based on their bin in each Dimension.  For each set we build a list
         # of points.
-        bins, inverse, counts = np.unique(
-            binned_gridpoints,
-            return_inverse=True,
-            return_counts=True,
-            axis=0)
+        bins, inverse, counts = np.unique(binned_gridpoints, return_inverse=True,
+                                          return_counts=True, axis=0)
 
         # inverse is now a "unique bin number" for each point gridpoints
         # we want to turn that into a list of points for each bin
         # so we argsort
         inverse_argsort = np.argsort(inverse).astype(np.int32)
         cumulative_counts = np.cumsum(counts)
-        gp_map = {
-            tuple(bi): inverse_argsort[cci-ci:cci]
-            for bi, cci, ci in zip(bins, cumulative_counts, counts, strict=True)
-        }
+        gp_map = {tuple(bi): inverse_argsort[cci-ci:cci]
+                  for bi, cci, ci in zip(bins, cumulative_counts, counts, strict=True)}
 
         # the result is now going to be a concatenation of these lists
         # for each of the output ranks
@@ -2054,30 +1941,22 @@ def _rank_to_points(self):
 
         global_rank_to_bins = {}
 
-        from itertools import product
         for bi in bins:
             # This is a list of sets for the Dimension-specific rank
-            dim_rank_sets = [
-                dgdr[bii]
-                for dgdr, bii in zip(dim_group_dim_rank, bi, strict=True)
-            ]
+            dim_rank_sets = [dgdr[bii]
+                             for dgdr, bii in zip(dim_group_dim_rank, bi, strict=True)]
 
             # Convert these to an absolute rank
             # This is where we will throw a KeyError if there are points OOB
             for dim_ranks in product(*dim_rank_sets):
                 global_rank = dim_ranks_to_glb[tuple(dim_ranks)]
-                global_rank_to_bins\
-                    .setdefault(global_rank, set())\
-                    .add(tuple(bi))
+                global_rank_to_bins.setdefault(global_rank, set()).add(tuple(bi))
 
         empty = np.array([], dtype=np.int32)
 
-        return [
-            np.concatenate(
-                (empty, *[gp_map[bi] for bi in global_rank_to_bins.get(rank, [])])
-            )
-            for rank in range(distributor.comm.Get_size())
-        ]
+        return [np.concatenate((empty, *[gp_map[bi]
+                                         for bi in global_rank_to_bins.get(rank, [])]))
+                for rank in range(distributor.comm.Get_size())]
 
     def _build_par_dim_to_nnz(self, active_gp, active_mrow):
         # The case where we parallelise over a non-local index is suboptimal, but
@@ -2092,12 +1971,10 @@ def _build_par_dim_to_nnz(self, active_gp, active_mrow):
             nnz_M = active_mrow.size - 1
             return {
                 self._par_dim_to_nnz_map: np.arange(active_mrow.size, dtype=np.int32),
-                self._par_dim_to_nnz_m: np.zeros(
-                    (self.grid.shape_local[pardim_index],), dtype=np.int32
-                ),
-                self._par_dim_to_nnz_M: np.full(
-                    (self.grid.shape_local[pardim_index],), nnz_M, dtype=np.int32
-                ),
+                self._par_dim_to_nnz_m: np.zeros((self.grid.shape_local[pardim_index],),
+                                                 dtype=np.int32),
+                self._par_dim_to_nnz_M: np.full((self.grid.shape_local[pardim_index],),
+                                                nnz_M, dtype=np.int32),
             }
 
         # Get the radius along the parallel Dimension
@@ -2143,9 +2020,7 @@ def manual_scatter(self, *, data_all_zero=False):
             self.scattered_data = self.data
             self.scatter_result = {
                 self: self.data,
-                **{
-                    getattr(self, k): getattr(self, k).data for k in self._sub_functions
-                },
+                **{getattr(self, k): getattr(self, k).data for k in self._sub_functions},
                 self.mrow: self.mrow.data,
                 self.mcol: self.mcol.data,
                 self.mval: self.mval.data,
@@ -2175,11 +2050,8 @@ def manual_scatter(self, *, data_all_zero=False):
         r_tuple = tuple(self.r[dim] for dim in self.grid.dimensions)
 
         npoint, nloc, nnz, ndim, r_tuple_bcast, nt = distributor.comm.bcast(
-            (self.npoint,
-             self._gridpoints.data.shape[0],
-             m_coo.nnz,
-             self._gridpoints.data.shape[-1],
-             r_tuple,
+            (self.npoint, self._gridpoints.data.shape[0], m_coo.nnz,
+             self._gridpoints.data.shape[-1], r_tuple,
              self.data.shape[self._time_position]), root=0)
 
         # important that all ranks have the same ndims and same r
@@ -2189,20 +2061,15 @@ def manual_scatter(self, *, data_all_zero=False):
         # handle None radius
         r_tuple_no_none = tuple(
             ri if ri is not None else self.grid.size_map[d].glb
-            for ri, d in zip(r_tuple, self.grid.dimensions, strict=True)
-        )
+            for ri, d in zip(r_tuple, self.grid.dimensions, strict=True))
 
         # now all ranks can allocate the buffers to receive into
         if distributor.myrank != 0:
-            if data_all_zero:
-                scattered_data = np.zeros([nt, npoint], dtype=self.dtype)
-            else:
-                scattered_data = np.empty([nt, npoint], dtype=self.dtype)
+            scattered_data = (np.zeros if data_all_zero else np.empty)([nt, npoint],
+                                                                       dtype=self.dtype)
             scattered_gp = np.empty([nloc, ndim], dtype=np.int32)
-            scattered_coeffs = [
-                np.empty([nloc, r_tuple_no_none[idim]], dtype=self.dtype)
-                for idim in range(ndim)
-            ]
+            scattered_coeffs = [np.empty([nloc, r_tuple_no_none[idim]], dtype=self.dtype)
+                                for idim in range(ndim)]
             scattered_mrow = np.empty([nnz], dtype=np.int32)
             scattered_mcol = np.empty([nnz], dtype=np.int32)
             scattered_mval = np.empty([nnz], dtype=self.dtype)
@@ -2211,9 +2078,8 @@ def manual_scatter(self, *, data_all_zero=False):
 
             # These are copies because we mess with them down below
             scattered_gp = self._gridpoints.data.copy()
-            scattered_coeffs = [
-                self.interpolation_coefficients[d].data.copy()
-                for d in self.grid.dimensions]
+            scattered_coeffs = [self.interpolation_coefficients[d].data.copy()
+                                for d in self.grid.dimensions]
             scattered_mrow = m_coo.row.copy()
             scattered_mcol = m_coo.col.copy()
             scattered_mval = m_coo.data.copy()
@@ -2250,9 +2116,8 @@ def manual_scatter(self, *, data_all_zero=False):
                 effective_gridpoints = np.zeros_like(effective_gridpoints)
 
             # rewrite the matrix to remove the rows in groups 0 and 4
-            mask = (
-                (effective_gridpoints >= _left - this_dim_r + 1)
-                & (effective_gridpoints < _right))
+            mask = ((effective_gridpoints >= _left - this_dim_r + 1)
+                    & (effective_gridpoints < _right))
 
             which = np.nonzero(mask)
             active_mrow = active_mrow[which]
@@ -2297,10 +2162,8 @@ def manual_scatter(self, *, data_all_zero=False):
         self.scatter_result = {
             self: scattered_data,
             self.gridpoints: scattered_gp,
-            **{
-                self.interpolation_coefficients[d]: scattered_coeffs[idim]
-                for idim, d in enumerate(self.grid.dimensions)
-            },
+            **{self.interpolation_coefficients[d]: scattered_coeffs[idim]
+               for idim, d in enumerate(self.grid.dimensions)},
             self.mrow: active_mrow,
             self.mcol: active_mcol,
             self.mval: active_mval,
@@ -2335,19 +2198,12 @@ def manual_gather(self):
 
         # This relies on all ranks having a copy of all data. Which feels "bad".
         if distributor.myrank != 0:
-            distributor.comm.Reduce(
-                self.scattered_data,
-                None,
-                op=MPI.SUM,
-                root=0
-            )
+            distributor.comm.Reduce(self.scattered_data, None, op=MPI.SUM, root=0)
         else:
             distributor.comm.Reduce(
                 MPI.IN_PLACE,
                 self.scattered_data,  # Note: on rank 0 data === scattered_data.
-                op=MPI.SUM,
-                root=0
-            )
+                op=MPI.SUM, root=0)
 
     def _dist_gather(self, data):
         pass

From 305c6b50b575104f01dac261295ae528969c9aa2 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 15:57:14 -0300
Subject: [PATCH 3/9] misc: Reduce types/dense

---
 devito/types/dense.py | 160 ++++++++++++++----------------------------
 1 file changed, 52 insertions(+), 108 deletions(-)

diff --git a/devito/types/dense.py b/devito/types/dense.py
index 30f543a47c..de44417ac9 100644
--- a/devito/types/dense.py
+++ b/devito/types/dense.py
@@ -114,9 +114,8 @@ def __init_finalize__(self, *args, function=None, **kwargs):
                 # case `self._data is None`
                 _ = self.data
         else:
-            raise ValueError(
-                f'`initializer` must be callable or buffer, not {type(initializer)}'
-            )
+            raise ValueError('`initializer` must be callable or buffer, not '
+                             f'{type(initializer)}')
 
     _subs = Differentiable._subs
 
@@ -170,10 +169,7 @@ def __dtype_setup__(cls, **kwargs):
         dtype = kwargs.get('dtype')
         if dtype is not None:
             return dtype
-        elif grid is not None:
-            return grid.dtype
-        else:
-            return np.float32
+        return grid.dtype if grid is not None else np.float32
 
     def __coefficients_setup__(self, **kwargs):
         """
@@ -184,10 +180,8 @@ def __coefficients_setup__(self, **kwargs):
             if coeffs == 'symbolic':
                 _ = deprecations.symbolic_warn
             else:
-                raise ValueError(
-                    f'coefficients must be one of {str(fd_weights_registry)}'
-                    f' not {coeffs}'
-                )
+                raise ValueError(f'coefficients must be one of {str(fd_weights_registry)}'
+                                 f' not {coeffs}')
         return coeffs
 
     @cached_property
@@ -253,10 +247,8 @@ def shape_with_halo(self):
         the outhalo of boundary ranks contains a number of elements depending
         on the rank position in the decomposed grid (corner, side, ...).
         """
-        return tuple(
-            j + i + k
-            for i, (j, k) in zip(self.shape, self._size_outhalo, strict=True)
-        )
+        return tuple(j + i + k
+                     for i, (j, k) in zip(self.shape, self._size_outhalo, strict=True))
 
     @cached_property
     def _shape_with_inhalo(self):
@@ -271,10 +263,7 @@ def _shape_with_inhalo(self):
         Typically, this property won't be used in user code, but it may come
         in handy for testing or debugging
         """
-        return tuple(
-            j + i + k
-            for i, (j, k) in zip(self.shape, self._halo, strict=True)
-        )
+        return tuple(j + i + k for i, (j, k) in zip(self.shape, self._halo, strict=True))
 
     @cached_property
     def shape_allocated(self):
@@ -286,13 +275,10 @@ def shape_allocated(self):
         -----
         In an MPI context, this is the *local* with_halo region shape.
         """
-        return DimensionTuple(
-            *[
-                j + i + k
-                for i, (j, k) in zip(self._shape_with_inhalo, self._padding, strict=True)
-            ],
-            getters=self.dimensions
-        )
+        return DimensionTuple(*[
+            j + i + k
+            for i, (j, k) in zip(self._shape_with_inhalo, self._padding, strict=True)
+        ], getters=self.dimensions)
 
     @cached_property
     def shape_global(self):
@@ -319,12 +305,8 @@ def shape_global(self):
 
     @property
     def symbolic_shape(self):
-        return DimensionTuple(
-            *[
-                self._C_get_field(FULL, d).size for d in self.dimensions
-            ],
-            getters=self.dimensions
-        )
+        return DimensionTuple(*[self._C_get_field(FULL, d).size for d in self.dimensions],
+                              getters=self.dimensions)
 
     # `dimension_shape` exposes the per-Dimension symbolic size. For dense
     # Functions it coincides with `symbolic_shape`; the alias exists because
@@ -358,16 +340,11 @@ def _size_outhalo(self):
             # and inhalo correspond
             return self._size_inhalo
 
-        left = [
-            abs(min(i.loc_abs_min-i.glb_min-j, 0))
-            if i and not i.loc_empty else 0
-            for i, j in zip(self._decomposition, self._size_inhalo.left, strict=True)
-        ]
-        right = [
-            max(i.loc_abs_max+j-i.glb_max, 0)
-            if i and not i.loc_empty else 0
-            for i, j in zip(self._decomposition, self._size_inhalo.right, strict=True)
-        ]
+        left = [abs(min(i.loc_abs_min-i.glb_min-j, 0)) if i and not i.loc_empty else 0
+                for i, j in zip(self._decomposition, self._size_inhalo.left, strict=True)]
+        right = [max(i.loc_abs_max+j-i.glb_max, 0) if i and not i.loc_empty else 0
+                 for i, j in zip(self._decomposition, self._size_inhalo.right,
+                                 strict=True)]
 
         sizes = tuple(Size(i, j) for i, j in zip(left, right, strict=True))
 
@@ -382,23 +359,13 @@ def _size_outhalo(self):
                 if not self._distributor.is_boundary_rank:
                     warning(' '.join(wrap(warning_msg)))
                 else:
-                    left_dist = [
-                        i
-                        for i, d in zip(left, self.dimensions, strict=True)
-                        if d in self._distributor.dimensions
-                    ]
-                    right_dist = [
-                        i
-                        for i, d in zip(right, self.dimensions, strict=True)
-                        if d in self._distributor.dimensions
-                    ]
-                    for i, j, k, l in zip(
-                        left_dist,
-                        right_dist,
-                        self._distributor.mycoords,
-                        self._distributor.topology,
-                        strict=False
-                    ):
+                    left_dist = [i for i, d in zip(left, self.dimensions, strict=True)
+                                 if d in self._distributor.dimensions]
+                    right_dist = [i for i, d in zip(right, self.dimensions, strict=True)
+                                  if d in self._distributor.dimensions]
+                    for i, j, k, l in zip(left_dist, right_dist,
+                                          self._distributor.mycoords,
+                                          self._distributor.topology, strict=False):
                         if l > 1 and ((j > 0 and k == 0) or (i > 0 and k == l-1)):
                             warning(' '.join(wrap(warning_msg)))
                             break
@@ -424,26 +391,20 @@ def _mask_modulo(self):
     @cached_property
     def _mask_domain(self):
         """Slice-based mask to access the domain region of the allocated data."""
-        return tuple(
-            slice(i, j)
-            for i, j in zip(self._offset_domain, self._offset_halo.right, strict=True)
-        )
+        return tuple(slice(i, j) for i, j in zip(self._offset_domain,
+                                                 self._offset_halo.right, strict=True))
 
     @cached_property
     def _mask_inhalo(self):
         """Slice-based mask to access the domain+inhalo region of the allocated data."""
-        return tuple(
-            slice(i.left, i.right + j.right)
-            for i, j in zip(self._offset_inhalo, self._size_inhalo, strict=True)
-        )
+        return tuple(slice(i.left, i.right + j.right)
+                     for i, j in zip(self._offset_inhalo, self._size_inhalo, strict=True))
 
     @cached_property
     def _mask_outhalo(self):
         """Slice-based mask to access the domain+outhalo region of the allocated data."""
-        return tuple(
-            slice(i.start - j.left, i.stop and i.stop + j.right or None)
-            for i, j in zip(self._mask_domain, self._size_outhalo, strict=True)
-        )
+        return tuple(slice(i.start - j.left, i.stop and i.stop + j.right or None)
+                     for i, j in zip(self._mask_domain, self._size_outhalo, strict=True))
 
     @cached_property
     def _decomposition(self):
@@ -464,11 +425,8 @@ def _decomposition_outhalo(self):
         """
         if self._distributor is None:
             return (None,)*self.ndim
-        return tuple(
-            v.reshape(*self._size_inhalo[d])
-            if v is not None else v
-            for d, v in zip(self.dimensions, self._decomposition, strict=True)
-        )
+        return tuple(v.reshape(*self._size_inhalo[d]) if v is not None else v
+                     for d, v in zip(self.dimensions, self._decomposition, strict=True))
 
     @property
     def data(self):
@@ -720,8 +678,7 @@ def local_indices(self):
     def initializer(self):
         if isinstance(self._data, np.ndarray):
             return self.data_with_halo.view(np.ndarray)
-        else:
-            return self._initializer
+        return self._initializer
 
     _C_structname = 'dataobj'
     _C_field_data = 'data'
@@ -760,9 +717,7 @@ def _C_make_dataobj(self, alias=None, **args):
 
         # MPI-related fields
         dataobj._obj.npsize = (c_ulong*self.ndim)(*[
-            i - sum(j)
-            for i, j in zip(data.shape, self._size_padding, strict=True)
-        ])
+            i - sum(j) for i, j in zip(data.shape, self._size_padding, strict=True)])
         dataobj._obj.dsize = (c_ulong*self.ndim)(*self._size_domain)
         dataobj._obj.hsize = (c_int*(self.ndim*2))(*flatten(self._size_halo))
         dataobj._obj.hofs = (c_int*(self.ndim*2))(*flatten(self._offset_halo))
@@ -834,10 +789,8 @@ def _C_get_field(self, region, dim, side=None):
 
     def _halo_exchange(self):
         """Perform the halo exchange with the neighboring processes."""
-        if not MPI.Is_initialized() or \
-                MPI.COMM_WORLD.size == 1 or \
-                not configuration['mpi'] or \
-                self.grid is None:
+        if (not MPI.Is_initialized() or MPI.COMM_WORLD.size == 1 or
+                not configuration['mpi'] or self.grid is None):
             # Nothing to do
             return
         if MPI.COMM_WORLD.size > 1 and self._distributor is None:
@@ -1192,14 +1145,11 @@ def __indices_setup__(cls, *args, **kwargs):
         if args:
             assert len(args) == len(dimensions)
             staggered_indices = tuple(args)
+        elif not staggered:
+            staggered_indices = dimensions
         else:
-            if not staggered:
-                staggered_indices = dimensions
-            else:
-                staggered_indices = (
-                    d + i * d.spacing / 2
-                    for d, i in zip(dimensions, staggered, strict=True)
-                )
+            staggered_indices = (d + i * d.spacing / 2
+                                 for d, i in zip(dimensions, staggered, strict=True))
         return tuple(dimensions), tuple(staggered_indices)
 
     @property
@@ -1269,8 +1219,7 @@ def __halo_setup__(self, **kwargs):
                    len(space_halo) != len(self.space_dimensions):
                     raise TypeError("Invalid `space_order`")
                 v = list(space_halo)
-                halo = [v.pop(0) if i.is_Space else (0, 0)
-                        for i in self.dimensions]
+                halo = [v.pop(0) if i.is_Space else (0, 0) for i in self.dimensions]
 
             else:
                 raise TypeError("Invalid `space_order`")
@@ -1348,9 +1297,8 @@ def _arg_check(self, args, intervals, **kwargs):
         """
         data = args[self.name]
 
-        if args.options['index-mode'] == 'int32' and \
-           args.options['linearize'] and \
-           data.size - 1 >= np.iinfo(np.int32).max:
+        if (args.options['index-mode'] == 'int32' and args.options['linearize'] and
+                data.size - 1 >= np.iinfo(np.int32).max):
             raise InvalidArgument(f"`{self.name}`, with its {data.size} elements, is too "
                                   "big for int32 pointer arithmetic. Consider using the "
                                   "'index-mode=int64' option, the save=Buffer(..) "
@@ -1516,9 +1464,8 @@ def __indices_setup__(cls, *args, **kwargs):
             dimensions = list(Function.__indices_setup__(**kwargs)[0])
             dimensions.insert(cls._time_position, time_dim)
 
-        return Function.__indices_setup__(
-            *args, dimensions=dimensions, staggered=kwargs.get('staggered')
-        )
+        return Function.__indices_setup__(*args, dimensions=dimensions,
+                                          staggered=kwargs.get('staggered'))
 
     @classmethod
     def __shape_setup__(cls, **kwargs):
@@ -1548,9 +1495,7 @@ def __shape_setup__(cls, **kwargs):
             raise TypeError("`dimensions` required if both `grid` and "
                             "`shape` are provided")
         else:
-            shape = super().__shape_setup__(
-                grid=grid, shape=shape, dimensions=dimensions
-            )
+            shape = super().__shape_setup__(grid=grid, shape=shape, dimensions=dimensions)
 
         return tuple(shape)
 
@@ -1648,14 +1593,13 @@ def _halo_exchange(self):
 
     def _arg_values(self, estimate_memory=False, **kwargs):
         if self._parent is not None and self.parent.name not in kwargs:
-            return self._parent._arg_defaults(
-                alias=self._parent, estimate_memory=estimate_memory
-            ).reduce_all()
+            return self._parent._arg_defaults(alias=self._parent,
+                                              estimate_memory=estimate_memory
+                                              ).reduce_all()
         elif self.name in kwargs:
             raise RuntimeError(f"`{self.name}` is a SubFunction, so it can't be assigned "
                                "a value dynamically")
-        else:
-            return self._arg_defaults(alias=self, estimate_memory=estimate_memory)
+        return self._arg_defaults(alias=self, estimate_memory=estimate_memory)
 
     def _arg_apply(self, *args, **kwargs):
         if self._parent is not None:

From 194002b265fee8b0faed92c5b6fee71fd1a693f4 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 16:02:45 -0300
Subject: [PATCH 4/9] misc: Reduce operator

---
 devito/operator/operator.py | 146 +++++++++++-------------------------
 1 file changed, 44 insertions(+), 102 deletions(-)

diff --git a/devito/operator/operator.py b/devito/operator/operator.py
index c97870f075..1ab228cee8 100644
--- a/devito/operator/operator.py
+++ b/devito/operator/operator.py
@@ -256,10 +256,8 @@ def _build(cls, expressions, **kwargs):
 
         # References to local or external routines
         op._func_table = OrderedDict()
-        op._func_table.update(OrderedDict([(i, MetaCall(None, False))
-                                           for i in profiler._ext_calls]))
-        op._func_table.update(OrderedDict([(i.root.name, i)
-                                           for i in byproduct.funcs]))
+        op._func_table.update((i, MetaCall(None, False)) for i in profiler._ext_calls)
+        op._func_table.update((i.root.name, i) for i in byproduct.funcs)
 
         # Internal mutable state to store information about previous runs,
         # autotuning reports, etc
@@ -382,9 +380,7 @@ def _lower_exprs(cls, expressions, **kwargs):
         # in particular uniqueness across expressions is ensured
         expressions = concretize_subdims(expressions, **kwargs)
 
-        processed = [LoweredEq(i) for i in expressions]
-
-        return processed
+        return [LoweredEq(i) for i in expressions]
 
     # Compilation -- Cluster level
 
@@ -459,9 +455,7 @@ def _lower_stree(cls, clusters, **kwargs):
         # Build a ScheduleTree from a sequence of Clusters
         stree = stree_build(clusters, **kwargs)
 
-        stree = cls._specialize_stree(stree)
-
-        return stree
+        return cls._specialize_stree(stree)
 
     # Compilation -- Iteration/Expression tree level
 
@@ -547,9 +541,7 @@ def dimensions(self):
         dimensions = FindSymbols('dimensions').visit(self)
         ret.update(d for d in dimensions if d.is_PerfKnob)
 
-        ret = tuple(sorted(ret, key=attrgetter('name')))
-
-        return ret
+        return tuple(sorted(ret, key=attrgetter('name')))
 
     @cached_property
     def input(self):
@@ -593,8 +585,6 @@ def _prepare_arguments(self, autotune=None, estimate_memory=False, **kwargs):
                 if k not in self._known_arguments:
                     raise InvalidArgument(f"Unrecognized argument `{k}={v}`")
 
-        overrides, defaults = split(self.input, lambda p: p.name in kwargs)
-
         # DiscreteFunctions may be created from CartesianDiscretizations, which in
         # turn could be Grids or SubDomains. Both may provide arguments
         discretizations = {getattr(kwargs.get(p.name, p), 'grid', None)
@@ -619,17 +609,14 @@ def _prepare_arguments(self, autotune=None, estimate_memory=False, **kwargs):
         # Pre-process Dimension overrides. This may help ruling out ambiguities
         # when processing the `defaults` arguments. A topological sorting is used
         # as DerivedDimensions may depend on their parents
-        edges = [(i, i.parent) for i in nodes
-                 if i.is_Derived and i.parent in nodes]
+        edges = [(i, i.parent) for i in nodes if i.is_Derived and i.parent in nodes]
         toposort = DAG(nodes, edges).topological_sort()
 
         # Prepare to process data-carriers
         args = kwargs['args'] = ReducerMap()
 
-        kwargs['metadata'] = {'language': self._language,
-                              'platform': self._platform,
-                              'transients': self.transients,
-                              **self.threads_info}
+        kwargs['metadata'] = {'language': self._language, 'platform': self._platform,
+                              'transients': self.transients, **self.threads_info}
 
         overrides, defaults = split(self.input, lambda p: p.name in kwargs)
 
@@ -640,9 +627,8 @@ def _prepare_arguments(self, autotune=None, estimate_memory=False, **kwargs):
                 args.reduce_inplace()
             except ValueError as e:
                 v = [i for i in overrides if i.name in args]
-                raise InvalidArgument(
-                    f"Override `{p}` is incompatible with overrides `{v}`"
-                ) from e
+                raise InvalidArgument(f"Override `{p}` is incompatible with overrides "
+                                      f"`{v}`") from e
 
         # Process data-carrier defaults
         for p in defaults:
@@ -705,16 +691,13 @@ def _prepare_arguments(self, autotune=None, estimate_memory=False, **kwargs):
 
         # Sanity check
         for p in self.parameters:
-            p._arg_check(args, self._dspace[p], am=self._access_modes.get(p),
-                         **kwargs)
+            p._arg_check(args, self._dspace[p], am=self._access_modes.get(p), **kwargs)
         for d in self.dimensions:
-            try:
+            with suppress(AttributeError):
                 if d.is_Space and any(self._dspace[d].offsets):
                     warn(f"Shrinking bounds (`{d.min_name}`, `{d.max_name}`); "
                          f"some `{d}` points will not be computed. Likely "
                          "insufficient space_order for the derivatives.")
-            except AttributeError:
-                pass
             if d.is_Derived:
                 d._arg_check(args)
 
@@ -753,8 +736,7 @@ def _postprocess_errors(self, retval, comm=None):
             raise ExecutionError(
                 "Kernel launch failed due to insufficient resources. This may be "
                 "due to excessive register pressure in one of the Operator "
-                "kernels. Try supplying a smaller `par-tile` value."
-            )
+                "kernels. Try supplying a smaller `par-tile` value.")
         elif retval == error_mapper['KernelLaunchClusterConfig']:
             raise ExecutionError(
                 "Kernel launch failed due to an invalid thread block cluster "
@@ -762,8 +744,7 @@ def _postprocess_errors(self, retval, comm=None):
                 "does not perfectly divide the number of blocks launched for a "
                 "kernel. This is a known, strong limitation which effectively "
                 "prevents the use of `tbc-tile` in realistic scenarios, but it "
-                "will be removed in future versions."
-            )
+                "will be removed in future versions.")
         elif retval == error_mapper['KernelLaunchUnknown']:
             raise ExecutionError(
                 "Kernel launch failed due to an unknown error. This might "
@@ -951,8 +932,7 @@ def estimate_memory(self, **kwargs):
         memreport = {'host': mem[host_layer], 'device': mem[device_layer]}
 
         # Extra information for enriched Operators
-        extras = self._enrich_memreport(args)
-        memreport.update(extras)
+        memreport.update(self._enrich_memreport(args))
 
         return MemoryEstimate(memreport, name=self.name)
 
@@ -1042,8 +1022,7 @@ def apply(self, **kwargs):
                 argnum = int(e.args[0][9:].split(':')[0]) - 1
                 raise ctypes.ArgumentError(
                     f"error in argument '{self.parameters[argnum].name}' with value"
-                    f" '{arg_values[argnum]}': {e.args[0]}"
-                ) from e
+                    f" '{arg_values[argnum]}': {e.args[0]}") from e
             else:
                 raise
 
@@ -1160,10 +1139,7 @@ def lower_perfentry(v):
             if v.gpointss:
                 values.append(f"{fround(v.gpointss):.2f} GPts/s")
 
-            if values:
-                return f"[{', '.join(values)}]"
-            else:
-                return ""
+            return f"[{', '.join(values)}]" if values else ""
 
         for k, v in summary.items():
             rank = f"[rank{k.rank}]" if k.rank is not None else ''
@@ -1252,8 +1228,7 @@ def __setstate__(self, state):
             self._lib.name = soname
 
         self._allocator = default_allocator(
-            f'{type(self._compiler).__name__}.{self._language}.{self._platform}'
-        )
+            f'{type(self._compiler).__name__}.{self._language}.{self._platform}')
 
 
 # *** Recursive compilation ("rcompile") machinery
@@ -1291,11 +1266,8 @@ def compile(self, **kwargs):
 # (because once, during the main compilation phase, is simply enough), but also
 # dangerous as some of them (the minority) might break in some circumstances
 # if applied in cascade (e.g., `linearization` on top of `linearization`)
-rcompile_registry = {
-    'avoid_denormals': False,
-    'linearize': False,
-    'place-transfers': False
-}
+rcompile_registry = {'avoid_denormals': False, 'linearize': False,
+                     'place-transfers': False}
 
 
 def rcompile(expressions, kwargs, options, target=None):
@@ -1322,9 +1294,7 @@ def rcompile(expressions, kwargs, options, target=None):
     irs, byproduct0 = RCompiles(expressions, cls).compile(**kwargs)
 
     key = lambda i: isinstance(i, (EntryFunction, DeviceFunction))
-    byproduct = byproduct0.filter(key)
-
-    return irs, byproduct
+    return irs, byproduct0.filter(key)
 
 
 # *** Misc helpers
@@ -1351,8 +1321,7 @@ def opkwargs(self):
         temp_registry = {v: k for k, v in compiler_registry.items()}
         compiler = temp_registry[self.compiler.__class__]
 
-        return {'platform': self.platform.name,
-                'compiler': compiler,
+        return {'platform': self.platform.name, 'compiler': compiler,
                 'language': self.language}
 
     @property
@@ -1381,15 +1350,11 @@ def saved_mapper(self):
         The number of saved TimeFunctions in the Operator, grouped by
         memory hierarchy layer.
         """
-        key0 = lambda f: (f.is_TimeFunction and
-                          f.save is not None and
+        key0 = lambda f: (f.is_TimeFunction and f.save is not None and
                           not isinstance(f.save, Buffer))
         functions = [f for f in self.op.input if key0(f)]
 
-        key1 = lambda f: f.layer
-        mapper = as_mapper(functions, key1)
-
-        return mapper
+        return as_mapper(functions, lambda f: f.layer)
 
     @cached_property
     def _op_symbols(self):
@@ -1402,10 +1367,8 @@ def _op_functions(self):
         return [i for i in self._op_symbols if i.is_DiscreteFunction and not i.alias]
 
     def _apply_override(self, i):
-        try:
-            return self.get(i.name, i)._obj
-        except AttributeError:
-            return self.get(i.name, i)
+        obj = self.get(i.name, i)
+        return getattr(obj, '_obj', obj)
 
     def _get_nbytes(self, i):
         """
@@ -1489,27 +1452,19 @@ def nbytes_avail_mapper(self):
                 # Since might not have this layer in the mapper
                 mapper[layer] -= self.nbytes_consumed_operator.get(layer, 0)
 
-        mapper = {k: int(v) for k, v in mapper.items()}
-
-        return mapper
+        return {k: int(v) for k, v in mapper.items()}
 
     @cached_property
     def nbytes_consumed(self):
         """Memory consumed by all objects in the Operator"""
-        mem_locations = (
-            self.nbytes_consumed_functions,
-            self.nbytes_consumed_arrays,
-            self.nbytes_consumed_memmapped
-        )
+        mem_locations = (self.nbytes_consumed_functions, self.nbytes_consumed_arrays,
+                         self.nbytes_consumed_memmapped)
         return {layer: sum(loc[layer] for loc in mem_locations) for layer in _layers}
 
     @cached_property
     def nbytes_consumed_operator(self):
         """Memory consumed by objects allocated within the Operator"""
-        mem_locations = (
-            self.nbytes_consumed_arrays,
-            self.nbytes_consumed_memmapped
-        )
+        mem_locations = (self.nbytes_consumed_arrays, self.nbytes_consumed_memmapped)
         return {layer: sum(loc[layer] for loc in mem_locations) for layer in _layers}
 
     @cached_property
@@ -1546,8 +1501,7 @@ def nbytes_consumed_arrays(self):
         # Temporaries such as Arrays are allocated and deallocated on-the-fly
         # while in C land, so they need to be accounted for as well
         for i in self._op_symbols:
-            if not i.is_Array or not i._mem_heap or i.alias \
-               or not i.is_regular:
+            if not i.is_Array or not i._mem_heap or i.alias or not i.is_regular:
                 continue
 
             nbytes = i.nbytes if i.is_regular else i.nbytes_max
@@ -1584,11 +1538,9 @@ def nbytes_consumed_memmapped(self):
             for i in self.op.input:
                 if not is_on_device(i, self.options['gpu-fit']):
                     continue
-                try:
+                with suppress(AttributeError):
                     if i._mem_mapped:
                         device += self._get_nbytes(i)
-                except AttributeError:
-                    pass
 
         return {disk_layer: 0, host_layer: 0, device_layer: device}
 
@@ -1597,14 +1549,12 @@ def nbytes_snapshots(self):
         disk = 0
         # Layers are sometimes aliases, so include aliases here
         for i in self._op_symbols:
-            try:
+            with suppress(AttributeError):
                 if i._child is None and i.alias is not True:
                     # Use only the "innermost" layer to avoid counting snapshots
                     # twice. This layer will have no child.
                     v = self._apply_override(i)
                     disk += v.size_snapshot*v._time_size_ideal*np.dtype(v.dtype).itemsize
-            except AttributeError:
-                pass
 
         return {disk_layer: disk, host_layer: 0, device_layer: 0}
 
@@ -1614,8 +1564,7 @@ def parse_kwargs(**kwargs):
     Parse keyword arguments provided to an Operator.
     """
     # `dse` -- deprecated, dropped
-    dse = kwargs.pop("dse", None)
-    if dse is not None:
+    if kwargs.pop("dse", None) is not None:
         warning("The `dse` argument is deprecated. "
                 "The optimization level is now controlled via the `opt` argument")
 
@@ -1675,11 +1624,9 @@ def parse_kwargs(**kwargs):
     # Handle deprecations
     deprecated_options = ('cire-mincost-inv', 'cire-mincost-sops', 'cire-maxalias')
     for i in deprecated_options:
-        try:
+        with suppress(KeyError):
             options.pop(i)
             warning(f"Ignoring deprecated optimization option `{i}`")
-        except KeyError:
-            pass
     kwargs['options'] = options
 
     # `opt`, mode
@@ -1714,8 +1661,7 @@ def parse_kwargs(**kwargs):
         kwargs['language'] = language
     elif kwopenmp is not None:
         # Handle deprecated `openmp` kwarg for backward compatibility
-        omp = {'C': 'openmp', 'CXX': 'CXXopenmp'}.get(configuration['language'],
-                                                      'openmp')
+        omp = {'C': 'openmp', 'CXX': 'CXXopenmp'}.get(configuration['language'], 'openmp')
         kwargs['language'] = omp if openmp else 'C'
     else:
         kwargs['language'] = configuration['language']
@@ -1727,15 +1673,13 @@ def parse_kwargs(**kwargs):
             raise ValueError("Argument `compiler` should be a `str`")
         if compiler not in configuration._accepted['compiler']:
             raise InvalidOperator(f"Illegal `compiler={str(compiler)}`")
-        kwargs['compiler'] = compiler_registry[compiler](platform=kwargs['platform'],
-                                                         language=kwargs['language'],
-                                                         mpi=configuration['mpi'],
-                                                         name=compiler)
+        kwargs['compiler'] = compiler_registry[compiler](
+            platform=kwargs['platform'], language=kwargs['language'],
+            mpi=configuration['mpi'], name=compiler)
     elif any([platform, language]):
-        kwargs['compiler'] =\
-            configuration['compiler'].__new_with__(platform=kwargs['platform'],
-                                                   language=kwargs['language'],
-                                                   mpi=configuration['mpi'])
+        kwargs['compiler'] = configuration['compiler'].__new_with__(
+            platform=kwargs['platform'], language=kwargs['language'],
+            mpi=configuration['mpi'])
     else:
         kwargs['compiler'] = configuration['compiler'].__new_with__()
 
@@ -1748,10 +1692,8 @@ def parse_kwargs(**kwargs):
 
     # `allocator`
     kwargs['allocator'] = default_allocator(
-        f"{kwargs['compiler'].__class__.__name__}"
-        f".{kwargs['language']}"
-        f".{kwargs['platform']}"
-    )
+        f"{kwargs['compiler'].__class__.__name__}.{kwargs['language']}"
+        f".{kwargs['platform']}")
 
     # Normalize `subs`, if any
     kwargs['subs'] = {k: sympify(v) for k, v in kwargs.get('subs', {}).items()}

From 06b49dee177d892aaa5d0e4664f7c9e129ba15f6 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 16:19:01 -0300
Subject: [PATCH 5/9] misc: Grid and differentiable reductions

---
 devito/finite_differences/differentiable.py | 147 ++++++--------------
 devito/types/grid.py                        | 145 +++++++------------
 2 files changed, 88 insertions(+), 204 deletions(-)

diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py
index d86734a454..7b000ef39e 100644
--- a/devito/finite_differences/differentiable.py
+++ b/devito/finite_differences/differentiable.py
@@ -1,4 +1,5 @@
 from collections import ChainMap
+from contextlib import suppress
 from functools import cached_property, singledispatch
 from itertools import product
 
@@ -28,16 +29,8 @@
 from devito.types.basic import AbstractFunction, Indexed
 
 __all__ = [
-    'Conj',
-    'DiffDerivative',
-    'Differentiable',
-    'EvalDerivative',
-    'Imag',
-    'IndexDerivative',
-    'IndexDerivativeProperty',
-    'LocalSum',
-    'Real',
-    'Weights',
+    'Conj', 'DiffDerivative', 'Differentiable', 'EvalDerivative', 'Imag',
+    'IndexDerivative', 'IndexDerivativeProperty', 'LocalSum', 'Real', 'Weights',
 ]
 
 
@@ -88,10 +81,7 @@ def grid(self):
         grids = {g.root for g in grids}
         if len(grids) > 1:
             warning("Expression contains multiple grids, returning first found")
-        try:
-            return grids.pop()
-        except KeyError:
-            return None
+        return grids.pop() if grids else None
 
     @cached_property
     def dtype(self):
@@ -154,12 +144,10 @@ def _fd(self):
         # Filter out all args with fd order too high
         fd_args = []
         for f in self._args_diff:
-            try:
+            with suppress(AttributeError):
                 if f.space_order <= self.space_order and \
                         (not f.is_TimeDependent or f.time_order <= self.time_order):
                     fd_args.append(f)
-            except AttributeError:
-                pass
         return dict(ChainMap(*[getattr(i, '_fd', {}) for i in fd_args]))
 
     @cached_property
@@ -168,10 +156,7 @@ def _symbolic_functions(self):
 
     @cached_property
     def function(self):
-        if len(self._functions) == 1:
-            return set(self._functions).pop()
-        else:
-            return None
+        return set(self._functions).pop() if len(self._functions) == 1 else None
 
     @cached_property
     def _uses_symbolic_coefficients(self):
@@ -198,10 +183,8 @@ def _subs(self, old, new, **hints):
             return self
         args = list(self.args)
         for i, arg in enumerate(args):
-            try:
+            with suppress(AttributeError):
                 args[i] = arg._subs(old, new, **hints)
-            except AttributeError:
-                continue
         return self.func(*args, evaluate=False)
 
     @property
@@ -290,10 +273,7 @@ def __rfloordiv__(self, other):
         return floor(other / self)
 
     def _inv(self, ref, safe=False):
-        if safe:
-            return SafeInv(self, ref or self)
-        else:
-            return 1 / self
+        return SafeInv(self, ref or self) if safe else 1 / self
 
     def __mod__(self, other):
         return Mod(self, other)
@@ -468,9 +448,8 @@ def has(self, *pattern):
         """
         for p in pattern:
             # Following sympy convention, return True if any is found
-            if isinstance(p, type) \
-                    and issubclass(p, sympy.Symbol) \
-                    and any(isinstance(i, p) for i in self.free_symbols):
+            if (isinstance(p, type) and issubclass(p, sympy.Symbol) and
+                    any(isinstance(i, p) for i in self.free_symbols)):
                 # Symbols (and subclasses) are the leaves of an expression, and they
                 # are promptly available via `free_symbols`. So this is super quick
                 return True
@@ -500,8 +479,7 @@ def deep_priority(expr):
     generic value, so `mu*tau_xx` reports .75 rather than `tau_xx`'s 2.1.
     """
     prio = getattr(expr, '_fd_priority', 0)
-    return max([prio] + [deep_priority(i)
-                         for i in getattr(expr, '_args_diff', ())])
+    return max([prio] + [deep_priority(i) for i in getattr(expr, '_args_diff', ())])
 
 
 def highest_priority(diff_op, candidates=None):
@@ -568,26 +546,9 @@ def _gather_for_diff(self):
     def _eval_is_even(self):
         return None
 
-    def _eval_is_odd(self):
-        return None
-
-    def _eval_is_integer(self):
-        return None
-
-    def _eval_is_negative(self):
-        return None
-
-    def _eval_is_extended_negative(self):
-        return None
-
-    def _eval_is_positive(self):
-        return None
-
-    def _eval_is_extended_positive(self):
-        return None
-
-    def _eval_is_zero(self):
-        return None
+    _eval_is_odd = _eval_is_integer = _eval_is_negative = _eval_is_even
+    _eval_is_extended_negative = _eval_is_positive = _eval_is_even
+    _eval_is_extended_positive = _eval_is_zero = _eval_is_even
 
 
 class DifferentiableFunction(DifferentiableOp):
@@ -656,14 +617,12 @@ def __new__(cls, *args, **kwargs):
         # which would destroy `EvalDerivative`s if present. So here we perform
         # a similar thing, but cautiously construct an evaluated Add, which
         # will preserve the integrity of `EvalDerivative`s, if any
-        try:
+        with suppress(AttributeError, ValueError):
             a, b = args
             if a.is_Rational:
                 r, b = b.as_coeff_Mul()
                 if r is sympy.S.One and type(b) is Add:
                     return Add(*[_keep_coeff(a, bi) for bi in b.args], evaluate=False)
-        except (AttributeError, ValueError):
-            pass
 
         return super().__new__(cls, *args, **kwargs)
 
@@ -686,9 +645,8 @@ def _gather_for_diff(self):
         derivs, other = split(self.args, lambda a: isinstance(a, sympy.Derivative))
         if len(derivs) == 0:
             return self._eval_at(highest_priority(self))
-        else:
-            other = self.func(*other)._eval_at(highest_priority(self))
-            return self.func(other, *derivs)
+        other = self.func(*other)._eval_at(highest_priority(self))
+        return self.func(other, *derivs)
 
     @classmethod
     def _off_func(cls, a, func):
@@ -777,8 +735,7 @@ def _eval_at(self, func, interp_mode='direct', **kwargs):
                                    if dim in func.indices_ref.getters})
             else:
                 source = a.indices_ref
-            new_factors.append(interp_at(a, source, block_indices,
-                                         self.interp_order))
+            new_factors.append(interp_at(a, source, block_indices, self.interp_order))
 
         # Final I from block's location to func
         return interp_at(self.func(*new_factors), block_indices,
@@ -836,8 +793,7 @@ class RealComplexPart(ComplexPart):
 
     @cached_property
     def dtype(self):
-        dtype = extract_dtype(self)
-        return dtype(0).real.__class__
+        return extract_dtype(self)(0).real.__class__
 
 
 class Real(RealComplexPart):
@@ -871,11 +827,9 @@ def __new__(cls, expr, dimensions, **kwargs):
         if not dimensions:
             return expr
         for d in dimensions:
-            try:
+            with suppress(AttributeError):
                 if d.is_Dimension and is_integer(d.symbolic_size):
                     continue
-            except AttributeError:
-                pass
             raise ValueError("Expected Dimension with numeric size, "
                              f"got `{d}` instead")
 
@@ -899,11 +853,8 @@ def __new__(cls, expr, dimensions, **kwargs):
         return obj
 
     def __repr__(self):
-        return "{}({}, ({}))".format(
-            self.__class__.__name__,
-            self.expr,
-            ', '.join(d.name for d in self.dimensions)
-        )
+        dims = ', '.join(d.name for d in self.dimensions)
+        return f"{self.__class__.__name__}({self.expr}, ({dims}))"
 
     __str__ = __repr__
 
@@ -935,11 +886,8 @@ def _evaluate(self, **kwargs):
             return self._rebuild(expr)
 
         values = product(*[list(d.range) for d in self.dimensions])
-        terms = []
-        for i in values:
-            mapper = dict(zip(self.dimensions, i, strict=True))
-            terms.append(expr.xreplace(mapper))
-        return sum(terms)
+        return sum([expr.xreplace(dict(zip(self.dimensions, i, strict=True)))
+                    for i in values])
 
     @property
     def bound_symbols(self):
@@ -1109,17 +1057,14 @@ def _xreplace(self, rule):
             return rule[self], True
         elif not rule:
             return self, False
-        else:
-            try:
-                weights, flags = zip(
-                    *[i._xreplace(rule) for i in self.weights], strict=True
-                )
-                if any(flags):
-                    return self.func(initvalue=weights, function=None), True
-            except AttributeError:
-                # `float` weights
-                pass
-            return super()._xreplace(rule)
+        try:
+            weights, flags = zip(*[i._xreplace(rule) for i in self.weights], strict=True)
+            if any(flags):
+                return self.func(initvalue=weights, function=None), True
+        except AttributeError:
+            # `float` weights
+            pass
+        return super()._xreplace(rule)
 
     @cached_property
     def _npweights(self):
@@ -1136,8 +1081,7 @@ def value(self, idx):
             v = self._npweights[idx]
         if v.is_Number or v.is_Indexed:
             return sympy.sympify(v)
-        else:
-            return self[idx]
+        return self[idx]
 
 
 class IndexDerivativeProperty(Tag):
@@ -1196,8 +1140,7 @@ def compare(self, other):
             return (self.weights.compare(other.weights) or
                     self.base.compare(other.base) or
                     (p1 > p2) - (p1 < p2))
-        else:
-            return super().compare(other)
+        return super().compare(other)
 
     @cached_property
     def base(self):
@@ -1250,8 +1193,7 @@ def _subs(self, old, new, **hints):
         # may fail to identify the sub-expression to be replaced (note: if
         # `a/b/c` are atoms or Indexeds, it's generally fine)
 
-        if not old.is_Mul or \
-           old is not self.base:
+        if not old.is_Mul or old is not self.base:
             return super()._subs(old, new, **hints)
 
         return self._rebuild(new * self.weights)
@@ -1300,7 +1242,6 @@ def __new__(cls, *args, base=None, **kwargs):
             # story: a zero-order derivative whose weights collapse to one is
             # the identity, so it comes back as the sum it was applied to.
             assert len(args) <= 1
-            return obj
 
         return obj
 
@@ -1339,8 +1280,7 @@ class diffify:
 
     def __new__(cls, obj):
         args = [diffify._doit(i) for i in obj.args]
-        obj = diffify._doit(obj, args)
-        return obj
+        return diffify._doit(obj, args)
 
     def _doit(obj, args=None):
         cls = diffify._cls(obj)
@@ -1401,8 +1341,7 @@ def _diff2sympy(obj):
 
         # Handle special objects
         if isinstance(obj, DiffDerivative):
-            return IndexDerivative(*args, obj.mapper,
-                                   deriv_order=obj.deriv_order,
+            return IndexDerivative(*args, obj.mapper, deriv_order=obj.deriv_order,
                                    properties=obj.properties), True
 
         # Handle generic objects such as arithmetic operations
@@ -1421,8 +1360,7 @@ def _diff2sympy(obj):
                 # In case of indices using other Function, evaluate
                 # may not be a supported argument.
                 return obj.func(*args), True
-        else:
-            return obj, False
+        return obj, False
 
     return _diff2sympy(expr)[0]
 
@@ -1460,8 +1398,7 @@ def _(expr, x0, **kwargs):
     def test0(a):
         return all(a.indices[d] is i for d, i in x0.items() if d in a.dimensions)
 
-    oa, ia = split(expr._args_diff,
-                   lambda a: isinstance(a, sympy.Derivative) or test0(a))
+    oa, ia = split(expr._args_diff, lambda a: isinstance(a, sympy.Derivative) or test0(a))
     oa = oa + tuple(a for a in expr.args if a not in expr._args_diff)
 
     # Interpolate the necessary args
@@ -1476,8 +1413,7 @@ def test0(a):
 def _(expr, x0, **kwargs):
     if expr.args:
         return expr.func(*[interp_for_fd(i, x0, **kwargs) for i in expr.args])
-    else:
-        return expr
+    return expr
 
 
 @interp_for_fd.register(AbstractFunction)
@@ -1486,5 +1422,4 @@ def _(expr, x0, **kwargs):
                and expr.indices.get(d, v) is not v}
     if x0_expr:
         return expr.subs({expr.indices[d]: v for d, v in x0_expr.items()})
-    else:
-        return expr
+    return expr
diff --git a/devito/types/grid.py b/devito/types/grid.py
index 5884b0fc30..7af2442282 100644
--- a/devito/types/grid.py
+++ b/devito/types/grid.py
@@ -68,6 +68,14 @@ def dtype(self):
     def root(self):
         return self
 
+    @property
+    def size_map(self):
+        """Map between SpaceDimensions and their global/local size."""
+        return {
+            d: GlobalLocal(g, l)
+            for d, g, l in zip(self.dimensions, self.shape, self.shape_local, strict=True)
+        }
+
 
 class Grid(CartesianDiscretization, ArgProvider):
 
@@ -159,14 +167,10 @@ def __init__(self, shape, extent=None, origin=None, dimensions=None,
             ndim = len(shape)
             assert ndim <= 3
             dim_names = self._default_dimensions[:ndim]
-            dim_spacing = tuple(
-                Spacing(name=f'h_{n}', dtype=dtype, is_const=True)
-                for n in dim_names
-            )
-            dimensions = tuple(
-                SpaceDimension(name=n, spacing=s)
-                for n, s in zip(dim_names, dim_spacing, strict=True)
-            )
+            dim_spacing = tuple(Spacing(name=f'h_{n}', dtype=dtype, is_const=True)
+                                for n in dim_names)
+            dimensions = tuple(SpaceDimension(name=n, spacing=s)
+                               for n, s in zip(dim_names, dim_spacing, strict=True))
         else:
             for d in dimensions:
                 if not d.is_Space:
@@ -175,22 +179,18 @@ def __init__(self, shape, extent=None, origin=None, dimensions=None,
                 if d.is_Derived and not d.is_Conditional:
                     raise ValueError(f"Cannot create Grid with derived Dimension `{d}` "
                                      f"of type `{type(d)}`")
-            dimensions = dimensions
 
         super().__init__(shape, dimensions, dtype)
 
         # Create a Distributor, used internally to implement domain decomposition
         # by all Functions defined on this Grid
         topology = topology or configuration['topology']
-        if topology:
-            if len(topology) == len(self.shape):
-                self._topology = topology
-            else:
-                warning(f"Ignoring the provided topology `{topology}` as it "
-                        f"is incompatible with the grid shape `{self.shape}`")
-                self._topology = None
-        else:
-            self._topology = None
+        self._topology = None
+        if topology and len(topology) == len(self.shape):
+            self._topology = topology
+        elif topology:
+            warning(f"Ignoring the provided topology `{topology}` as it "
+                    f"is incompatible with the grid shape `{self.shape}`")
         self._distributor = Distributor(shape, dimensions, comm, self._topology)
 
         # The physical extent and grid spacing
@@ -271,12 +271,8 @@ def origin_ioffset(self):
     def origin_offset(self):
         """Physical offset of the local (per-process) origin from the domain origin."""
         return DimensionTuple(
-            *[
-                i*h
-                for i, h in zip(self.origin_ioffset, self.spacing, strict=True)
-            ],
-            getters=self.dimensions
-        )
+            *[i*h for i, h in zip(self.origin_ioffset, self.spacing, strict=True)],
+            getters=self.dimensions)
 
     @property
     def time_dim(self):
@@ -331,9 +327,8 @@ def spacing_map(self):
                 # the SpaceDimensions
                 mapper[d.spacing] = s
             else:
-                raise AssertionError(
-                    'Cannot map between spacing symbol for SpaceDimension'
-                )
+                raise AssertionError('Cannot map between spacing symbol for '
+                                     'SpaceDimension')
 
         return mapper
 
@@ -342,14 +337,6 @@ def shape_local(self):
         """Shape of the local (per-process) physical domain."""
         return self._distributor.shape
 
-    @property
-    def size_map(self):
-        """Map between SpaceDimensions and their global/local size."""
-        return {
-            d: GlobalLocal(g, l)
-            for d, g, l in zip(self.dimensions, self.shape, self.shape_local, strict=True)
-        }
-
     @property
     def topology(self):
         """The topology used for decomposing the CartesianDiscretization."""
@@ -374,8 +361,7 @@ def is_distributed(self, dim):
 
     @cached_property
     def _arg_names(self):
-        ret = []
-        ret.append(self.time_dim.spacing.name)
+        ret = [self.time_dim.spacing.name]
         ret.extend([i.name for i in self.origin_map])
         for i in self.spacing_map:
             try:
@@ -541,26 +527,22 @@ def is_distributed(self, dim):
         True if `dim` is a distributed Dimension for this CartesianDiscretization,
         False otherwise.
         """
-        if self.grid:
-            return any(dim is d for d in self.distributor.dimensions)
-        return False
+        return bool(self.grid) and any(dim is d for d in self.distributor.dimensions)
 
     @property
     def comm(self):
         """The MPI communicator inherited from the distributor."""
         if self.grid:
             return self.grid.comm
-        raise ValueError(
-            f'`SubDomain` {self.name} has no `Grid` attached and thus no `comm`'
-        )
+        raise ValueError(f'`SubDomain` {self.name} has no `Grid` attached and thus '
+                         'no `comm`')
 
     def _arg_values(self, **kwargs):
         try:
             return self.grid._arg_values(**kwargs)
         except AttributeError as e:
-            raise AttributeError(
-                f'{self} is not attached to a Grid and has no _arg_values'
-            ) from e
+            raise AttributeError(f'{self} is not attached to a Grid and has no '
+                                 '_arg_values') from e
 
 
 class SubDomain(AbstractSubDomain):
@@ -628,12 +610,9 @@ def __subdomain_finalize_legacy__(self, grid):
         # Create the SubDomain's SubDimensions
         sub_dimensions = []
         sdshape = []
-        for k, v, s in zip(
-            self.define(grid.dimensions).keys(),
-            self.define(grid.dimensions).values(),
-            grid.shape,
-            strict=True
-        ):
+        for k, v, s in zip(self.define(grid.dimensions).keys(),
+                           self.define(grid.dimensions).values(),
+                           grid.shape, strict=True):
             if isinstance(v, Dimension):
                 sub_dimensions.append(v)
                 sdshape.append(s)
@@ -645,22 +624,18 @@ def __subdomain_finalize_legacy__(self, grid):
                         raise ValueError(f"Expected side 'middle', not `{side}`")
                     sub_dimensions.append(SubDimension.middle(f'i{k.name}',
                                                               k, ltkn, rtkn))
-                    thickness = s-ltkn-rtkn
-                    sdshape.append(thickness)
+                    sdshape.append(s-ltkn-rtkn)
                 except ValueError:
                     side, thickness = v
                     constructor = {'left': SubDimension.left,
                                    'right': SubDimension.right}.get(side)
                     if constructor is None:
-                        raise ValueError(
-                            f"Expected sides 'left|right', not `{side}`"
-                        ) from None
+                        raise ValueError(f"Expected sides 'left|right', not "
+                                         f"`{side}`") from None
 
                     if s - thickness < 0:
-                        raise ValueError(
-                            f"Maximum thickness of dimension {k.name} "
-                            f"is {s}, not {thickness}"
-                        ) from None
+                        raise ValueError(f"Maximum thickness of dimension {k.name} "
+                                         f"is {s}, not {thickness}") from None
                     sub_dimensions.append(constructor(f'i{k.name}', k, thickness))
                     sdshape.append(thickness)
 
@@ -672,14 +647,6 @@ def __subdomain_finalize_legacy__(self, grid):
     def shape_local(self):
         return self._shape_local
 
-    @property
-    def size_map(self):
-        """Map between SpaceDimensions and their global/local size."""
-        return {
-            d: GlobalLocal(g, l)
-            for d, g, l in zip(self.dimensions, self.shape, self.shape_local, strict=True)
-        }
-
     def define(self, dimensions):
         """
         Parametrically describe the SubDomain w.r.t. a generic Grid.
@@ -696,17 +663,14 @@ def _arg_names(self):
         try:
             ret = self.grid._arg_names
         except AttributeError as e:
-            raise AttributeError(
-                f'{self} is not attached to a Grid and has no _arg_names'
-            ) from e
+            raise AttributeError(f'{self} is not attached to a Grid and has no '
+                                 '_arg_names') from e
 
         # Names for SubDomain thicknesses
         thickness_names = tuple([k.name for k in d._thickness_map]
                                 for d in self.dimensions if d.is_Sub)
 
-        ret += tuple(thickness_names)
-
-        return ret
+        return ret + tuple(thickness_names)
 
     def __getstate__(self):
         state = self.__dict__.copy()
@@ -859,23 +823,16 @@ def __init__(self, **kwargs):
         self._global_bounds = kwargs.get('bounds')
         super().__init__(**kwargs)
 
-        try:
-            _ = self.implicit_dimension
+        if hasattr(self, 'implicit_dimension'):
             warning("`implicit_dimension` is deprecated. You may safely remove it "
                     "from the class definition")
-        except AttributeError:
-            pass
 
     def __subdomain_finalize_core__(self, grid):
         self._dtype = grid.dtype
 
         # Compute the SubDomainSet shapes
-        global_bounds = []
-        for i in self._global_bounds:
-            if isinstance(i, int):
-                global_bounds.append(np.full(self._n_domains, i, dtype=np.int32))
-            else:
-                global_bounds.append(i)
+        global_bounds = [np.full(self._n_domains, i, dtype=np.int32)
+                         if isinstance(i, int) else i for i in self._global_bounds]
         d_m = global_bounds[0::2]
         d_M = global_bounds[1::2]
         shapes = []
@@ -1138,8 +1095,7 @@ def _build_domains(self, grid: Grid) -> tuple[int, tuple[np.ndarray]]:
         """
         if self.corners == 'overlap':
             return self._build_domains_overlap(grid)
-        else:
-            return self._build_domains_nooverlap(grid)
+        return self._build_domains_nooverlap(grid)
 
     def _build_domains_overlap(self, grid: Grid) -> tuple[int, tuple[np.ndarray]]:
 
@@ -1177,13 +1133,8 @@ def _build_domains_nooverlap(self, grid: Grid) -> tuple[int, tuple[np.ndarray]]:
         # Unpack the user-provided specification into a set of sides (on which
         # a cartesian product is taken) and a mapper from those sides to a set of
         # bounds for each dimension.
-        for d, s, b, i in zip(
-            grid.dimensions,
-            grid.shape,
-            self.border,
-            self.inset,
-            strict=True
-        ):
+        for d, s, b, i in zip(grid.dimensions, grid.shape, self.border, self.inset,
+                              strict=True):
             if d in self.border_dims:
                 side = self.border_dims[d]
 
@@ -1218,10 +1169,8 @@ def _build_domains_nooverlap(self, grid: Grid) -> tuple[int, tuple[np.ndarray]]:
                 # Don't add any domains that are completely centered
                 if self.corners != 'nocorners' or any(i is CENTER for i in d):
                     # Don't add corners if 'no corners' option selected
-                    domains.append([
-                        interval_map[dim][dom]
-                        for (dim, dom) in zip(grid.dimensions, d, strict=True)
-                    ])
+                    domains.append([interval_map[dim][dom] for (dim, dom)
+                                    in zip(grid.dimensions, d, strict=True)])
 
         domains = np.array(domains)
 

From f75d3d62db6dda57797a6746fc35108bcf6515a4 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 16:24:39 -0300
Subject: [PATCH 6/9] misc: Reduce data structures

---
 devito/tools/data_structures.py | 96 +++++++++------------------------
 1 file changed, 24 insertions(+), 72 deletions(-)

diff --git a/devito/tools/data_structures.py b/devito/tools/data_structures.py
index ced3fd2d9b..f13daa52ea 100644
--- a/devito/tools/data_structures.py
+++ b/devito/tools/data_structures.py
@@ -11,20 +11,9 @@
 from devito.tools.algorithms import toposort
 from devito.tools.utils import as_tuple, filter_ordered, humanbytes
 
-__all__ = [
-    'DAG',
-    'Bunch',
-    'DefaultFrozenDict',
-    'DefaultOrderedDict',
-    'EnrichedTuple',
-    'MemoryEstimate',
-    'OrderedSet',
-    'Ordering',
-    'ReducerMap',
-    'UnboundTuple',
-    'UnboundedMultiTuple',
-    'frozendict',
-]
+__all__ = ['DAG', 'Bunch', 'DefaultFrozenDict', 'DefaultOrderedDict', 'EnrichedTuple',
+           'MemoryEstimate', 'OrderedSet', 'Ordering', 'ReducerMap', 'UnboundTuple',
+           'UnboundedMultiTuple', 'frozendict']
 
 
 class Bunch:
@@ -43,9 +32,7 @@ def __init__(self, **kwargs):
         self.__dict__.update(kwargs)
 
     def __repr__(self):
-        return "Bunch({})".format(
-            ", ".join(["{}={}".format(*i) for i in self.__dict__.items()])
-        )
+        return f"Bunch({', '.join(f'{k}={v}' for k, v in self.__dict__.items())})"
 
     def __iter__(self):
         yield from self.__dict__.values()
@@ -69,8 +56,7 @@ def __new__(cls, *items, getters=None, **kwargs):
 
     def _rebuild(self, *args, **kwargs):
         # Need to explicitly apply any additional attributes
-        _kwargs = dict(self.__dict__)
-        _kwargs.update(**kwargs)
+        _kwargs = {**self.__dict__, **kwargs}
 
         return super()._rebuild(*args, **_kwargs)
 
@@ -87,8 +73,7 @@ def __getitem__(self, key):
             kwargs = dict(self.__dict__)
             kwargs['getters'] = list(self.getters)[start:stop]
             return EnrichedTuple(*items, **kwargs)
-        else:
-            return self.__getitem_hook__(key)
+        return self.__getitem_hook__(key)
 
     def __getitem_hook__(self, key):
         return self.getters[key]
@@ -160,17 +145,13 @@ def compare_to_first(v):
             if isinstance(first, np.ndarray) or isinstance(v, np.ndarray):
                 return (first == v).all()
             elif isinstance(v, Set):
-                if isinstance(first, Set):
-                    return not v.isdisjoint(first)
-                else:
-                    return first in v
+                return not v.isdisjoint(first) if isinstance(first, Set) else first in v
             elif isinstance(first, Set):
                 return v in first
             elif isinstance(v, range):
                 if isinstance(first, range):
                     return first.stop > v.start or v.stop > first.start
-                else:
-                    return first >= v.start and first < v.stop
+                return first >= v.start and first < v.stop
             elif isinstance(first, range):
                 return v >= first.start and v < first.stop
             else:
@@ -184,10 +165,8 @@ def compare_to_first(v):
                 if not isinstance(c, range):
                     return c
             return candidates[0]
-        else:
-            raise ValueError(
-                f'Unable to find unique value for key {key}, candidates: {candidates}'
-            )
+        raise ValueError(f'Unable to find unique value for key {key}, candidates: '
+                         f'{candidates}')
 
     def reduce(self, key, op=None):
         """
@@ -209,8 +188,7 @@ def reduce(self, key, op=None):
         if op is None:
             # Return a unique value if it exists
             return self.unique(key)
-        else:
-            return reduce(op, self.getall(key))
+        return reduce(op, self.getall(key))
 
     def reduce_all(self):
         """Returns a dictionary with reduced/unique values for all keys."""
@@ -228,8 +206,7 @@ def reduce_inplace(self):
 class DefaultOrderedDict(OrderedDict):
     # Source: http://stackoverflow.com/a/6190500/562769
     def __init__(self, default_factory=None, *a, **kw):
-        if (default_factory is not None and
-           not isinstance(default_factory, Callable)):
+        if default_factory is not None and not isinstance(default_factory, Callable):
             raise TypeError('first argument must be callable')
         OrderedDict.__init__(self, *a, **kw)
         self.default_factory = default_factory
@@ -435,10 +412,7 @@ def roots(self):
 
     @property
     def edges(self):
-        ret = []
-        for k, v in self.graph.items():
-            ret.extend([(k, i) for i in v])
-        return tuple(ret)
+        return tuple((k, i) for k, v in self.graph.items() for i in v)
 
     @property
     def size(self):
@@ -532,8 +506,7 @@ def all_downstreams(self, node):
                     nodes_seen.add(downstream_node)
                     nodes.append(downstream_node)
             i += 1
-        return list(filter(lambda node: node in nodes_seen,
-                           self.topological_sort()))
+        return list(filter(lambda node: node in nodes_seen, self.topological_sort()))
 
     def topological_sort(self, choose_element=None):
         """
@@ -580,8 +553,7 @@ def topological_sort(self, choose_element=None):
 
         if len(l) == len(self.graph):
             return l
-        else:
-            raise ValueError('graph is not acyclic')
+        raise ValueError('graph is not acyclic')
 
     def connected_components(self, enumerated=False):
         """
@@ -599,12 +571,8 @@ def connected_components(self, enumerated=False):
                 groups.append(found)
 
         if enumerated:
-            mapper = OrderedDict()
-            for n, g in enumerate(groups):
-                mapper.update({i: n for i in g})
-            return mapper
-        else:
-            return tuple(groups)
+            return OrderedDict((i, n) for n, g in enumerate(groups) for i in g)
+        return tuple(groups)
 
     def find_paths(self, node):
         if node not in self.graph:
@@ -697,10 +665,7 @@ def __getitem__(self, key):
             if self._default is self._sentinel:
                 raise
 
-            if callable(self._default):
-                return self._default()
-            else:
-                return self._default
+            return self._default() if callable(self._default) else self._default
 
     def get(self, key, default=None):
         return self._dict.get(key, default)
@@ -779,14 +744,9 @@ class UnboundTuple(tuple):
     """
 
     def __new__(cls, *items, **kwargs):
-        nitems = []
-        for i in as_tuple(items):
-            if isinstance(i, UnboundTuple):
-                nitems.append(i)
-            elif isinstance(i, Iterable):
-                nitems.append(UnboundTuple(*i))
-            else:
-                nitems.append(i)
+        nitems = [UnboundTuple(*i)
+                  if isinstance(i, Iterable) and not isinstance(i, UnboundTuple) else i
+                  for i in as_tuple(items)]
 
         obj = super().__new__(cls, tuple(nitems))
         obj.last = len(nitems)
@@ -819,8 +779,7 @@ def __len__(self):
         return self.last
 
     def __repr__(self):
-        sitems = [s.__repr__() for s in self]
-        return "{}({})".format(self.__class__.__name__, ", ".join(sitems))
+        return f"{self.__class__.__name__}({', '.join(repr(s) for s in self)})"
 
     def __getitem__(self, idx):
         if not self:
@@ -833,10 +792,7 @@ def __getitem__(self, idx):
             step = idx.step or 1
             return UnboundTuple(*[self[i] for i in range(start, stop, step)])
         try:
-            if idx >= self.last-1:
-                return super().__getitem__(self.last-1)
-            else:
-                return super().__getitem__(idx)
+            return super().__getitem__(min(idx, self.last-1))
         except TypeError:
             # Slice, ...
             return UnboundTuple(self[i] for i in idx)
@@ -905,12 +861,8 @@ def index(self, item):
         return self.index(item)
 
     def iter(self):
-        if self.current is None:
-            self.current = 0
-        else:
-            self.current = min(self.current + 1, self.last - 1)
+        self.current = 0 if self.current is None else min(self.current + 1, self.last - 1)
         self[self.current].reset()
-        return
 
     def next(self):
         if not self:

From 422277d31c8679a58dec6d31f38b1cc0cc6f806e Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 16:38:12 -0300
Subject: [PATCH 7/9] misc: Reduce interpolators

---
 devito/operations/interpolators.py | 139 +++++++++++------------------
 1 file changed, 51 insertions(+), 88 deletions(-)

diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py
index ff4e14b1de..aa1e2b1514 100644
--- a/devito/operations/interpolators.py
+++ b/devito/operations/interpolators.py
@@ -57,8 +57,7 @@ def wrapper(interp, *args, **kwargs):
         for f in a_sfuncs:
             for s in f._sub_functions:
                 if getattr(f, s, None) not in subfuncs:
-                    raise ValueError(f"Interpolation/injection with {sfunc}"
-                                     f"requires {f} "
+                    raise ValueError(f"Interpolation/injection with {sfunc}requires {f} "
                                      f"to use the same {s} as {sfunc}")
 
         return func(interp, *args, **kwargs)
@@ -79,9 +78,7 @@ def _extract_subdomain(variables):
     if len(sdms) > 1:
         raise NotImplementedError("Sparse operation on multiple Functions defined on"
                                   " different SubDomains currently unsupported")
-    elif len(sdms) == 1:
-        return sdms.pop()
-    return None
+    return sdms.pop() if sdms else None
 
 
 class UnevaluatedSparseOperation(sympy.Expr, Evaluable, Pickable):
@@ -263,8 +260,7 @@ def _gdims(self):
     @property
     def _cdim(self):
         """Base CustomDimensions used to construct _rdim"""
-        dims = [self.sfunction._crdim(d) for d in self._gdims]
-        return dims
+        return [self.sfunction._crdim(d) for d in self._gdims]
 
     def _field_shifts(self, field):
         """
@@ -311,6 +307,7 @@ def _rdim(self, subdomain=None, shifts=None):
         return DimensionTuple(*rdims, getters=gdims)
 
     def _augment_implicit_dims(self, implicit_dims, extras=None):
+        extra = ()
         if extras is not None:
             # If variables are defined on a SubDomain of the Grid, then omit the
             # dimensions of that SubDomain from any extra dimensions found
@@ -322,12 +319,9 @@ def _augment_implicit_dims(self, implicit_dims, extras=None):
                                       if d.is_Sub and d.root in self._gdims])
 
             gdims = filter_ordered(edims + list(self._gdims))
-            extra = filter_ordered([i for v in extras for i in v.dimensions
-                                    if i not in gdims and
-                                    i not in self.sfunction.dimensions])
-            extra = tuple(extra)
-        else:
-            extra = tuple()
+            extra = tuple(filter_ordered([i for v in extras for i in v.dimensions
+                                          if i not in gdims and
+                                          i not in self.sfunction.dimensions]))
 
         if self.sfunction._sparse_position == -1:
             idims = self.sfunction.dimensions + as_tuple(implicit_dims) + extra
@@ -351,13 +345,11 @@ def _generate_gridpoints(self, key):
 
         name = f'{self.sfunction.name}_gp{_shift_tag(as_list(key))}'
         sfdim = self.sfunction._sparse_dim
-        ddim = CustomDimension(f'{name}d', 0, self.grid.dim - 1,
-                               self.grid.dim, sfdim)
+        ddim = CustomDimension(f'{name}d', 0, self.grid.dim - 1, self.grid.dim, sfdim)
         return Gridpoints(name=name, dtype=np.int32,
                           shape=(self.sfunction.npoint, self.grid.dim),
                           dimensions=(sfdim, ddim), space_order=0,
-                          alias=self.sfunction.alias,
-                          parent=self.sfunction)
+                          alias=self.sfunction.alias, parent=self.sfunction)
 
     def _gridpoints(self, shifts=None):
         return self._generate_gridpoints(tuple(shifts) if shifts else None)
@@ -378,8 +370,7 @@ def _positions(self, implicit_dims, shifts=None):
         gp = self._gridpoints(shifts=shifts)
         ddim = gp.dimensions[-1]
         return [Eq(p, gp._subs(ddim, di), implicit_dims=implicit_dims)
-                for (di, p) in enumerate(
-                    self.sfunction._pos_symbols(shifts=shifts))]
+                for (di, p) in enumerate(self.sfunction._pos_symbols(shifts=shifts))]
 
     def _coeff_data(self, coords, grid, shifts, spacing, origin):
         """
@@ -401,8 +392,7 @@ def _arg_defaults(self, coords=None, sfunc=None, origin=None):
         # Fp64 grid geometry -- avoids fp32 rounding on cell boundaries.
         grid = sfunc.grid
         spacing = np.array([as_fp64_decimal(h) for h in grid.spacing])
-        origin = np.array([as_fp64_decimal(o)
-                           for o in (origin or grid.origin)])
+        origin = np.array([as_fp64_decimal(o) for o in (origin or grid.origin)])
 
         args = {}
         for key in self._shifts_used or {None}:
@@ -420,8 +410,7 @@ def _arg_defaults(self, coords=None, sfunc=None, origin=None):
 
         return args
 
-    def _interp_idx(self, variables, implicit_dims=None, subdomain=None,
-                    shifts=None):
+    def _interp_idx(self, variables, implicit_dims=None, subdomain=None, shifts=None):
         """
         Generate interpolation indices for the DiscreteFunctions in `variables`.
 
@@ -443,11 +432,8 @@ def _interp_idx(self, variables, implicit_dims=None, subdomain=None,
         mapper = self._rdim(subdomain=subdomain, shifts=shifts).getters
 
         # Index substitution to make in variables
-        subs = {
-            ki: c + p
-            for ((k, c), p) in zip(mapper.items(), pos, strict=True)
-            for ki in {k, k.root}
-        }
+        subs = {ki: c + p for ((k, c), p) in zip(mapper.items(), pos, strict=True)
+                for ki in {k, k.root}}
 
         idx_subs = {v: v.subs(subs) for v in variables}
 
@@ -538,9 +524,8 @@ def _interpolate(self, expr, increment=False, self_subs=None, implicit_dims=None
         # Write/Incr `self`
         lhs = self.sfunction.subs(self_subs)
         ecls = Inc if increment else Eq
-        last = [ecls(lhs, rhs, implicit_dims=implicit_dims)]
 
-        return temps + last
+        return temps + [ecls(lhs, rhs, implicit_dims=implicit_dims)]
 
     def _inject(self, field, expr, increment=True, implicit_dims=None):
         """
@@ -597,9 +582,8 @@ def _inject(self, field, expr, increment=True, implicit_dims=None):
             # Can only be done for inject as interpolation needs a summing temp
             # that wouldn't allow collapsing
             with suppress(AttributeError):
-                implicit_dims = implicit_dims + tuple(r.parent for r in
-                                                      self._rdim(subdomain=subdomain,
-                                                                 shifts=shifts))
+                implicit_dims = implicit_dims + tuple(
+                    r.parent for r in self._rdim(subdomain=subdomain, shifts=shifts))
 
             # List of indirection indices for all adjacent grid points
             idx_subs, _temps = self._interp_idx(list(g_fields) + variables,
@@ -627,8 +611,7 @@ def _shift_values(shifts, grid, spacing):
     """Physical half-cell offsets for each grid dim, as fp64."""
     if not shifts:
         return np.zeros(grid.dim, dtype=np.float64)
-    subs = {d.spacing: float(h)
-            for d, h in zip(grid.dimensions, spacing, strict=True)}
+    subs = {d.spacing: float(h) for d, h in zip(grid.dimensions, spacing, strict=True)}
     return np.array([float(sympy.sympify(s).xreplace(subs)) for s in shifts])
 
 
@@ -710,7 +693,21 @@ def _sinc_weights(coords, grid, shifts, j, dtype, spacing, origin, r, b):
     return data
 
 
-class LinearInterpolator(WeightedInterpolator):
+class _TabulatedInterpolator(WeightedInterpolator):
+    """Shared plumbing for schemes whose weights are tabulated on the host."""
+
+    def _coeffs(self, shifts=None):
+        return self._generate_coeffs(tuple(shifts) if shifts else None)
+
+    @memoized_meth
+    def _weights(self, subdomain=None, shifts=None):
+        rdims = self._rdim(subdomain=subdomain, shifts=shifts)
+        coeffs = self._coeffs(shifts=shifts)
+        return Mul(*[w._subs(rd, rd - rd.parent.symbolic_min)
+                     for (rd, w) in zip(rdims, coeffs, strict=True)])
+
+
+class LinearInterpolator(_TabulatedInterpolator):
     """
     Linear (bilinear/trilinear) interpolator.
 
@@ -736,33 +733,15 @@ def _generate_coeffs(self, key):
         sfdim = self.sfunction._sparse_dim
 
         # Per-dim linear weights: `(npoint, 2)` holding `(1 - frac, frac)`.
-        return tuple(
-            Coeffs(name=f'{sfname}_w{d.name}{tag}',
-                   dtype=self._coeff_dtype,
-                   shape=(self.sfunction.npoint, 2),
-                   dimensions=(sfdim, r), space_order=0,
-                   alias=self.sfunction.alias,
-                   parent=self.sfunction)
-            for d, r in zip(self._gdims, self._cdim, strict=True)
-        )
-
-    def _coeffs(self, shifts=None):
-        return self._generate_coeffs(tuple(shifts) if shifts else None)
-
-    @memoized_meth
-    def _weights(self, subdomain=None, shifts=None):
-        rdims = self._rdim(subdomain=subdomain, shifts=shifts)
-        coeffs = self._coeffs(shifts=shifts)
-        return Mul(*[
-            w._subs(rd, rd - rd.parent.symbolic_min)
-            for (rd, w) in zip(rdims, coeffs, strict=True)
-        ])
+        return tuple(Coeffs(name=f'{sfname}_w{d.name}{tag}', dtype=self._coeff_dtype,
+                            shape=(self.sfunction.npoint, 2), dimensions=(sfdim, r),
+                            space_order=0, alias=self.sfunction.alias,
+                            parent=self.sfunction)
+                     for d, r in zip(self._gdims, self._cdim, strict=True))
 
     def _coeff_data(self, coords, grid, shifts, spacing, origin):
-        return {
-            w: _linear_weights(coords, grid, shifts, i, w.dtype, spacing, origin)
-            for i, w in enumerate(self._coeffs(shifts=shifts))
-        }
+        return {w: _linear_weights(coords, grid, shifts, i, w.dtype, spacing, origin)
+                for i, w in enumerate(self._coeffs(shifts=shifts))}
 
 
 class NearestInterpolator(LinearInterpolator):
@@ -810,10 +789,9 @@ def _positions(self, implicit_dims, shifts=None):
             # Only the coordinates are known, and the user-provided coefficients
             # are tied to the cell index the kernel derives from them
             return self._floor_positions(implicit_dims, shifts=shifts)
-        else:
-            # No position temp as we have directly the gridpoints
-            return[Eq(p, k, implicit_dims=implicit_dims)
-                   for (k, p) in self.sfunction._position_map(shifts=shifts).items()]
+        # No position temp as we have directly the gridpoints
+        return[Eq(p, k, implicit_dims=implicit_dims)
+               for (k, p) in self.sfunction._position_map(shifts=shifts).items()]
 
     def _arg_defaults(self, **kwargs):
         # Gridpoints and coefficients are user-provided SubFunctions of the
@@ -830,11 +808,10 @@ def _weights(self, subdomain=None, shifts=None):
         mappers = [{ddim: ri, cdim: rd-rd.parent.symbolic_min}
                    for (ri, rd) in enumerate(self._rdim(subdomain=subdomain,
                                                         shifts=shifts))]
-        return Mul(*[self.interpolation_coeffs.subs(mapper)
-                     for mapper in mappers])
+        return Mul(*[self.interpolation_coeffs.subs(mapper) for mapper in mappers])
 
 
-class SincInterpolator(WeightedInterpolator):
+class SincInterpolator(_TabulatedInterpolator):
     """
     Hicks windowed sinc interpolation scheme.
 
@@ -851,8 +828,7 @@ class SincInterpolator(WeightedInterpolator):
     _name = 'sinc'
 
     # Table 1
-    _b_table = {2: 2.94, 3: 4.53,
-                4: 4.14, 5: 5.26, 6: 6.40,
+    _b_table = {2: 2.94, 3: 4.53, 4: 4.14, 5: 5.26, 6: 6.40,
                 7: 7.51, 8: 8.56, 9: 9.56, 10: 10.64}
 
     def __init__(self, sfunction, shifts=()):
@@ -875,24 +851,11 @@ def _generate_coeffs(self, key):
         tag = _shift_tag(as_list(key))
         shape = (self.sfunction.npoint, 2 * self.r)
 
-        return tuple(
-            Coeffs(name=f'wsinc{r.name}{tag}', dtype=self._coeff_dtype,
-                   shape=shape, dimensions=(self.sfunction._sparse_dim, r),
-                   space_order=0, alias=self.sfunction.alias,
-                   parent=self.sfunction)
-            for r in self._cdim
-        )
-
-    def _coeffs(self, shifts=None):
-        return self._generate_coeffs(tuple(shifts) if shifts else None)
-
-    @memoized_meth
-    def _weights(self, subdomain=None, shifts=None):
-        rdims = self._rdim(subdomain=subdomain, shifts=shifts)
-        return Mul(*[
-            w._subs(rd, rd-rd.parent.symbolic_min)
-            for (rd, w) in zip(rdims, self._coeffs(shifts=shifts), strict=True)
-        ])
+        return tuple(Coeffs(name=f'wsinc{r.name}{tag}', dtype=self._coeff_dtype,
+                            shape=shape, dimensions=(self.sfunction._sparse_dim, r),
+                            space_order=0, alias=self.sfunction.alias,
+                            parent=self.sfunction)
+                     for r in self._cdim)
 
     def _coeff_data(self, coords, grid, shifts, spacing, origin):
         b = self._b_table[self.r]

From 4d844ce607f77fa6f506bd604cf4c34dc6cc4bd3 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 16:56:18 -0300
Subject: [PATCH 8/9] misc: Improve readability

---
 devito/finite_differences/differentiable.py | 34 +++++++++++++++++----
 devito/tools/data_structures.py             | 11 +++++--
 devito/types/basic.py                       |  1 +
 devito/types/dense.py                       |  6 ++--
 devito/types/dimension.py                   |  6 +++-
 devito/types/sparse.py                      |  6 ++--
 6 files changed, 49 insertions(+), 15 deletions(-)

diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py
index 7b000ef39e..143b91bba1 100644
--- a/devito/finite_differences/differentiable.py
+++ b/devito/finite_differences/differentiable.py
@@ -546,9 +546,26 @@ def _gather_for_diff(self):
     def _eval_is_even(self):
         return None
 
-    _eval_is_odd = _eval_is_integer = _eval_is_negative = _eval_is_even
-    _eval_is_extended_negative = _eval_is_positive = _eval_is_even
-    _eval_is_extended_positive = _eval_is_zero = _eval_is_even
+    def _eval_is_odd(self):
+        return None
+
+    def _eval_is_integer(self):
+        return None
+
+    def _eval_is_negative(self):
+        return None
+
+    def _eval_is_extended_negative(self):
+        return None
+
+    def _eval_is_positive(self):
+        return None
+
+    def _eval_is_extended_positive(self):
+        return None
+
+    def _eval_is_zero(self):
+        return None
 
 
 class DifferentiableFunction(DifferentiableOp):
@@ -827,9 +844,11 @@ def __new__(cls, expr, dimensions, **kwargs):
         if not dimensions:
             return expr
         for d in dimensions:
-            with suppress(AttributeError):
+            try:
                 if d.is_Dimension and is_integer(d.symbolic_size):
                     continue
+            except AttributeError:
+                pass
             raise ValueError("Expected Dimension with numeric size, "
                              f"got `{d}` instead")
 
@@ -886,8 +905,11 @@ def _evaluate(self, **kwargs):
             return self._rebuild(expr)
 
         values = product(*[list(d.range) for d in self.dimensions])
-        return sum([expr.xreplace(dict(zip(self.dimensions, i, strict=True)))
-                    for i in values])
+        terms = []
+        for i in values:
+            mapper = dict(zip(self.dimensions, i, strict=True))
+            terms.append(expr.xreplace(mapper))
+        return sum(terms)
 
     @property
     def bound_symbols(self):
diff --git a/devito/tools/data_structures.py b/devito/tools/data_structures.py
index f13daa52ea..b98e4f1b6e 100644
--- a/devito/tools/data_structures.py
+++ b/devito/tools/data_structures.py
@@ -744,9 +744,14 @@ class UnboundTuple(tuple):
     """
 
     def __new__(cls, *items, **kwargs):
-        nitems = [UnboundTuple(*i)
-                  if isinstance(i, Iterable) and not isinstance(i, UnboundTuple) else i
-                  for i in as_tuple(items)]
+        nitems = []
+        for i in as_tuple(items):
+            if isinstance(i, UnboundTuple):
+                nitems.append(i)
+            elif isinstance(i, Iterable):
+                nitems.append(UnboundTuple(*i))
+            else:
+                nitems.append(i)
 
         obj = super().__new__(cls, tuple(nitems))
         obj.last = len(nitems)
diff --git a/devito/types/basic.py b/devito/types/basic.py
index 699b644cca..ec2520c2bd 100644
--- a/devito/types/basic.py
+++ b/devito/types/basic.py
@@ -1166,6 +1166,7 @@ def dmap(self):
             return DeviceMap(f'd_{self.name}', shape=self._shape, function=self.function)
         elif self._mem_local:
             return self.indexed
+        return None
 
     @property
     def size(self):
diff --git a/devito/types/dense.py b/devito/types/dense.py
index de44417ac9..c2899d4e45 100644
--- a/devito/types/dense.py
+++ b/devito/types/dense.py
@@ -1593,9 +1593,9 @@ def _halo_exchange(self):
 
     def _arg_values(self, estimate_memory=False, **kwargs):
         if self._parent is not None and self.parent.name not in kwargs:
-            return self._parent._arg_defaults(alias=self._parent,
-                                              estimate_memory=estimate_memory
-                                              ).reduce_all()
+            return self._parent._arg_defaults(
+                alias=self._parent, estimate_memory=estimate_memory
+            ).reduce_all()
         elif self.name in kwargs:
             raise RuntimeError(f"`{self.name}` is a SubFunction, so it can't be assigned "
                                "a value dynamically")
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index 945b4f9a6b..d6aa433fef 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -1335,7 +1335,11 @@ def _rebuild_hierarchy(self, callback=None, step=None):
 
         name0 = pp.name
 
-        name1 = p.name if callback is None else callback(f'{callback(name0)}_blk')
+        if callback is None:
+            name1 = p.name
+        else:
+            base = callback(name0)
+            name1 = callback(f'{base}_blk')
 
         bd = p._rebuild(name1, pp, step=step or p.step)
 
diff --git a/devito/types/sparse.py b/devito/types/sparse.py
index fd81d53a37..b7b98c833d 100644
--- a/devito/types/sparse.py
+++ b/devito/types/sparse.py
@@ -2065,8 +2065,10 @@ def manual_scatter(self, *, data_all_zero=False):
 
         # now all ranks can allocate the buffers to receive into
         if distributor.myrank != 0:
-            scattered_data = (np.zeros if data_all_zero else np.empty)([nt, npoint],
-                                                                       dtype=self.dtype)
+            if data_all_zero:
+                scattered_data = np.zeros([nt, npoint], dtype=self.dtype)
+            else:
+                scattered_data = np.empty([nt, npoint], dtype=self.dtype)
             scattered_gp = np.empty([nloc, ndim], dtype=np.int32)
             scattered_coeffs = [np.empty([nloc, r_tuple_no_none[idim]], dtype=self.dtype)
                                 for idim in range(ndim)]

From e2105ac97224d294798f4453a0189cab3007ebd7 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Tue, 22 Sep 2026 17:01:29 -0300
Subject: [PATCH 9/9] misc: Remove underscores

---
 devito/operations/interpolators.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py
index aa1e2b1514..43ec2b7e3f 100644
--- a/devito/operations/interpolators.py
+++ b/devito/operations/interpolators.py
@@ -693,7 +693,7 @@ def _sinc_weights(coords, grid, shifts, j, dtype, spacing, origin, r, b):
     return data
 
 
-class _TabulatedInterpolator(WeightedInterpolator):
+class TabulatedInterpolator(WeightedInterpolator):
     """Shared plumbing for schemes whose weights are tabulated on the host."""
 
     def _coeffs(self, shifts=None):
@@ -707,7 +707,7 @@ def _weights(self, subdomain=None, shifts=None):
                      for (rd, w) in zip(rdims, coeffs, strict=True)])
 
 
-class LinearInterpolator(_TabulatedInterpolator):
+class LinearInterpolator(TabulatedInterpolator):
     """
     Linear (bilinear/trilinear) interpolator.
 
@@ -811,7 +811,7 @@ def _weights(self, subdomain=None, shifts=None):
         return Mul(*[self.interpolation_coeffs.subs(mapper) for mapper in mappers])
 
 
-class SincInterpolator(_TabulatedInterpolator):
+class SincInterpolator(TabulatedInterpolator):
     """
     Hicks windowed sinc interpolation scheme.