api: Fix mixed subdomain derivative by mloubout · Pull Request #3035 · devitocodes/devito

mloubout · 2026-09-28T19:57:52Z

No description provided.

codecov · 2026-09-28T20:16:50Z

Codecov Report

❌ Patch coverage is 99.29907% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 84.24%. Comparing base (64d15c5) to head (6573324).

Files with missing lines	Patch %	Lines
devito/passes/clusters/aliases.py	75.00%	1 Missing and 1 partial ⚠️
devito/types/grid.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3035      +/-   ##
==========================================
+ Coverage   84.12%   84.24%   +0.11%     
==========================================
  Files         258      258              
  Lines       56449    56854     +405     
  Branches     4817     4845      +28     
==========================================
+ Hits        47487    47894     +407     
+ Misses       8138     8136       -2     
  Partials      824      824

Flag	Coverage Δ
pytest-gpu-aomp-amdgpuX	`68.61% <58.10%> (-0.06%)`	⬇️
pytest-gpu-gcc-	`78.99% <99.29%> (+0.15%)`	⬆️
pytest-gpu-icx-	`78.93% <99.29%> (+0.17%)`	⬆️
pytest-gpu-nvc-nvidiaX	`69.14% <58.10%> (-0.05%)`	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

FabioLuporini · 2026-09-29T06:44:43Z

        obj._transpose = kwargs.get("transpose", direct)
        obj._method = kwargs.get("method", 'FD')
        obj._weights = cls._process_weights(**kwargs)
+        obj._halo = cls._validate_halo(kwargs.get("halo"))


would "support" be a better word?

I don't think so, since the "support region" conventionally refers to the footprint of the stencil itself, whilst this is for what happens at the edge of the domain

FabioLuporini · 2026-09-29T06:46:09Z

-        return self._subs(dim, dim + shift)
+        expr = self
+        for d in self.free_symbols:
+            if d is dim or not getattr(d, 'is_Space', False):


isn't it safer to check for if not isinstance(d, Dimensions): continue, then u don't need getattr?

in fact, u may also use retrieve_dimensions or something like that instead of free_symbols

FabioLuporini · 2026-09-29T06:55:46Z

        # A SubDistributor wraps an MPI communicator, which can't and shouldn't be pickled
        state.pop('_distributor', None)
+        # Memoized grown SubDomains and masks are rebuilt on demand
+        state.pop('_memoized_meth__cache_meth', None)


this looks fairly hacky

FabioLuporini · 2026-09-29T06:59:39Z

+                Max(0, Min(1, upper - point + 1)))
+
+
+class GrownSubDomain(SubDomain):


I don't see the need for this subclass

FabioLuporini · 2026-09-29T06:59:45Z

+        radius : frozendict of {Dimension: int}
+            Number of points to grow by, per root Dimension.
+        """
+        return GrownSubDomain(self, radius)


why do you need a subclass?

I agree - couldn't you just have a parent attribute on SubDomain or similar?

FabioLuporini · 2026-09-29T07:02:56Z

        # E.g., extract `1/h_x`
-        rule = lambda e: e.is_Pow and (not e.exp.is_Number or e.exp < 0)
+        rule = lambda e: (e.is_Pow and (not e.exp.is_Number or e.exp < 0) and
+                          not bare_stencil_dimensions(e))


same as before

FabioLuporini · 2026-09-29T07:03:41Z

+                pass
+            # Aliases are built by translating Indexeds; a term reading an
+            # iteration Dimension outside of any Indexed cannot be translated
+            if any(not d.is_Stencil for d in bare_dimensions(a)):


this is a correctness check , and so it doesn't belong here (see below...)

FabioLuporini · 2026-09-29T07:08:48Z

    return maybe_coeff, others


+def bare_dimensions(expr):


this thing you did with bare_dimensions is quite fishy, why do you need it ?

if the problem are unbound StencilDimensions , then you could use https://github.com/devitocodes/devito/blob/main/devito/ir/support/utils.py#L267

and the right place to discard alias candidates for correctness reasons is def _do_generate probably around line 290, when u see that ... & exclude

a check looking for unbound dimensions may have to be added there, and there only

FabioLuporini · 2026-09-29T07:11:12Z

+                           and e.halo is not None))
+
+    @staticmethod
+    def _extend(rhs, subdomain):


imho this doesn't belong here but should be part of the _eval_at machinery, where subdomain is now also available

rationale : we may want to be able to evaluate any expressions , not just Eq's, restricting to a subdomain

Seconded - the imports inside functions also implies that these are misplaced imo

FabioLuporini · 2026-09-29T07:14:10Z

        # E.g. `x + i*h_x` into `f(x)` s.t. `f(x + i*h_x)`
-        expr = expr._subs(dim, indices.expr)
+        expr = expr.shift(dim, indices.expr - dim)
+        if restricted:


potentially nitpicking: you could systematically pass a SubDomain , the whole grid if u don't have a SubDomain, with the grid's indicator being 1 everywhere; in that case, you would substantially reduce the special-casing

EdCaunt · 2026-09-29T09:55:21Z

+        """
+        return sympify((index - self.dim) / self.spacing)
+
+    @property


Cached property?

EdCaunt · 2026-09-29T09:57:37Z

-        rule = lambda e: e.is_Function or (e.is_Pow and e.exp.is_Number and 0 < e.exp < 1)
+        rule = lambda e: ((e.is_Function or
+                           (e.is_Pow and e.exp.is_Number and 0 < e.exp < 1)) and
+                          not bare_stencil_dimensions(e))


Readability is also becoming an issue here - a function would be better than a lambda. That being said, probably irrelevant given @FabioLuporini 's comments

EdCaunt · 2026-09-29T10:04:59Z

+        """True if the rhs has derivatives with halo=0."""
+        from devito.symbolics import search  # noqa
+
+        return bool(search(self.rhs, lambda e: getattr(e, 'is_Derivative', False)


Use retrieve_derivatives?

EdCaunt · 2026-09-29T10:06:07Z

+                           and e.halo is not None))
+
+    @staticmethod
+    def _extend(rhs, subdomain):


Seconded - the imports inside functions also implies that these are misplaced imo

EdCaunt · 2026-09-29T10:09:42Z

+        # The terms without indicators have no halo=0 derivative, so they are
+        # restricted to the SubDomain at the evaluation point
+        here = sympy.Mul(*[subdomain.indicator(d, 0, r) for d, r in radius.items()])
+        terms = [t if search(t, SubDomainIndicator) else here * t


search returns a set by default - is ordering determinism a problem here? Perhaps we should add a mode to search (ordered-unique) that returns an OrderedSet if so, since I'm pretty sure it will come in use

EdCaunt · 2026-09-29T10:10:52Z

 from itertools import product

 import numpy as np
+import sympy


as sp? Might as well remove the import below and use sympy.prod regardless

EdCaunt · 2026-09-29T10:28:06Z

+        radius : frozendict of {Dimension: int}
+            Number of points to grow by, per root Dimension.
+        """
+        return GrownSubDomain(self, radius)


I agree - couldn't you just have a parent attribute on SubDomain or similar?

EdCaunt · 2026-09-29T10:29:34Z

+        """
+        return GrownSubDomain(self, radius)
+
+    def indicator(self, dim, offset, radius):


Should this be a constructor class method? That being said, if this is the only way SubDomainIndicators are going to get constructed, then why not fold this into the __init__

EdCaunt · 2026-09-29T10:30:24Z

+    def radius(self):
+        return int(self.args[4])
+
+    @property


EdCaunt · 2026-09-29T10:31:41Z

+    def define(self, dimensions):
+        sizes = dict(zip(dimensions, self.subdomain.grid.shape, strict=True))
+        grown = {}
+        for d, v in self.subdomain.define(dimensions).items():
+            radius = self.radius.get(d, 0)
+            if isinstance(v, Dimension) or radius == 0:
+                grown[d] = v
+            elif v[0] == 'middle':
+                ltkn, rtkn = max(v[1] - radius, 0), max(v[2] - radius, 0)
+                grown[d] = d if ltkn == rtkn == 0 else ('middle', ltkn, rtkn)
+            else:
+                side, thickness = v
+                grown[d] = (side, min(thickness + radius, sizes[d]))
+        return grown


No need to go through define here - just construct the SubDimensions directly as needed. Probably irrelevant in the light of other comments however

FabioLuporini · 2026-09-30T08:33:37Z

+    @property
+    def subdomain(self):
+        """
+        With halo=0, the SubDomain outside of which the argument is treated as


just to be sure, do you actually need both halo and subdomain? or is subdomain just enough (essentially , when NOT None, it encodes the fact the users wants 0-valued taps outside of it) by any chance?

This is a good point

This only gets set at evaluate when it neeeds one. So conceptually it could but it would be very intricated and would need a lot of weird logic to use only one

I don't understand why that would be the case?

Because the subdomain cannot be constructed without knowing the radius

FabioLuporini · 2026-09-30T08:45:29Z

+    def _halo_radius(self):
+        """
+        Distance, per root Dimension, by which the evaluated derivatives with
+        halo=0 in the expression extend it past the SubDomain their argument is


a bit contrived explanation but maybe it's just me... could use an example

FabioLuporini · 2026-09-30T08:45:58Z


+    @cached_property
+    def _has_halo(self):
+        """True if the expression has derivatives with halo=0."""


nitpicking `halo=0`

FabioLuporini · 2026-09-30T08:46:45Z

+
+        expr = self
+        for d in retrieve_dimensions(self, mode='unique'):
+            if d is not dim and d.root is dim.root and (d.is_Sub or d is dim.root):


can u not simplify this with _defines somehow

_defines also contains ConditionalDimensions (xc._defines == {xc, x}), which must not be shifted, so the explicit check stays.

so to be sure the following would be borked:

if d in dim._defines and not d.is_NonlinearDerived

?

FabioLuporini · 2026-09-30T08:47:26Z

+        the argument as zero outside the equation's SubDomain) are supported.
+        """
+        if halo not in (None, 0):
+            raise ValueError(f"Expected halo=None or halo=0, got halo={halo}")


looks like this should be a boolean, not None or 0

It could but 0 is a better representation of what it does, which is to treat the halo as zeros.

it's just that u can't do not halo .. 😬 but OK

FabioLuporini · 2026-09-30T08:49:48Z

    is_commutative = True

-    __rkwargs__ = ('base',)
+    __rkwargs__ = ('base', '_halo_radius')


no prefix _

btw, u mean halo_map ?

since the need for this is essentially SubDomains, I'm not sure calling "halo-something" is actually informative, halo is such a generic and overloaded term... but again, it might just be me

subdomain_padding maybe?

It's not padding no, and halo is what is used throughout to define/deominate this part of a stencil.

FabioLuporini · 2026-09-30T08:52:25Z

+    The unbound StencilDimensions of `expr` appearing outside of any Indexed,
+    which translating `expr` into an alias would not shift.
+    """
+    return unbounded(expr) & set(retrieve_dimensions(expr, mode='unique'))


what's the & set(retrieve_dimensions(expr, mode='unique')) part for? isn't unbounded(expr) already a subset of retrieve_dimensions(expr) ?

Also retrieve_dimensions(expr, mode='unique') is already a set

isn't unbounded(expr) already a subset of retrieve_dimensions(expr)

No unbounded only look for IndexDerivative

FabioLuporini · 2026-09-30T09:00:50Z

        try:
            lhs = self.lhs._evaluate(**kwargs)
-            rhs = self.rhs._eval_at(self.lhs, **kwargs)._evaluate(**kwargs)
+            rhs = self.rhs._eval_at(self.lhs, **at, **kwargs)._evaluate(**kwargs)


I don't think u need this **at thing , in fact u don't need at at all -- just pass ..., subdomain=self.subdomain, **kwargs)....

FabioLuporini · 2026-09-30T09:02:21Z

        return eq

+    @cached_property
+    def _has_halo(self):


you also have _halo_halo in differentiable.py; if the recursion is done properly, you shouldn't need this method -- you can likely drop it; if you can't because something breaks, it's likely there's a deeper issue

anyway, I still think we can find a better alternative to "has_halo" , since it doesn't capture the SubDomain meaning, and halo is such an overloaded word

This is an Eq it's not differentiable

FabioLuporini · 2026-09-30T09:03:11Z

-        eq = self.func(lhs, rhs, subdomain=self.subdomain,
+
+        # ... and extend the rhs past it by the radius of their stencils
+        if restrict:


probably just if self._has_halo

EdCaunt · 2026-09-30T09:44:28Z

+    @property
+    def subdomain(self):
+        """
+        With halo=0, the SubDomain outside of which the argument is treated as


This is a good point

EdCaunt · 2026-09-30T10:23:32Z

    is_commutative = True

-    __rkwargs__ = ('base',)
+    __rkwargs__ = ('base', '_halo_radius')


subdomain_padding maybe?

EdCaunt · 2026-09-30T10:46:58Z

+    The unbound StencilDimensions of `expr` appearing outside of any Indexed,
+    which translating `expr` into an alias would not shift.
+    """
+    return unbounded(expr) & set(retrieve_dimensions(expr, mode='unique'))


Also retrieve_dimensions(expr, mode='unique') is already a set

EdCaunt · 2026-09-30T10:56:07Z

+        if self.parent is None:
+            return self.define(grid.dimensions)
+
+        sizes = dict(zip(grid.dimensions, grid.shape, strict=True))


grid.shape is already a DimensionTuple, so no need for this mapper

EdCaunt · 2026-09-30T11:05:41Z

+        Branch-free indicator of the point `offset` points away from the current
+        point along `dim` lying within this SubDomain: 1 inside, 0 outside.


This docstring is a little difficult to parse

EdCaunt · 2026-09-30T11:06:31Z

+        offset : expr-like
+            The offset, e.g. 2, or `i0` for a stencil in unexpanded form.
+        """
+        sub = self.dimension_map.get(dim, dim)


Nitpick: maybe_subdim =?

EdCaunt · 2026-09-30T11:14:48Z

+        for d in self.dimensions:
+            if d.is_Sub:
+                expr = expr * self.indicator(d.root, 0)


Is it worth having this product as a cached property? Since self.indicator = 1 when the dimension is not a SubDimension, is the special-casing needed?

Cache what product, this applies the indicator to the input expr, it's not a property.

EdCaunt · 2026-09-30T11:15:36Z

+        super().__init__(**kwargs)
+
+    def define(self, dimensions):
+        return {dimensions[0]: ('left', self.width)}


Does this not need to do anything for the other dimensions?

FabioLuporini · 2026-10-01T15:28:02Z

+
+    def __new__(cls, *args, halo_radius=None, **kwargs):
+        obj = super().__new__(cls, *args, **kwargs)
+        # With `halo=0`, the stencil radius (see `Differentiable.halo_radius`)


sorry for being pedantic, but imho this now wants a docstring. halo is undefined in this class scope, and having both halo and halo_radius doesn't help

FabioLuporini · 2026-10-01T15:34:07Z

+    does not look within Indexeds, so the intersection with `unbounded(expr)`
+    leaves out those appearing only in Indexeds, which are translated.
+    """
+    return unbounded(expr) & retrieve_dimensions(expr, mode='unique')


I'm not entirely sure if you're actually just trying to working around a bug inside unbounded

FabioLuporini · 2026-10-01T15:34:35Z

        if self._interp_mode is not None:
            kwargs['interp_mode'] = self._interp_mode
+
+        # Derivatives with `halo=0` treat their argument as zero outside the


this comment ain't useful here

FabioLuporini · 2026-10-01T15:39:17Z


                terms = cbk_compose(i)

+                # Aliases only translate Indexeds, so a term reading an unbound


can u give me an example of one such term?

…pressions A derivative of an expression mixing Grid Functions, e.g. f(x), with Functions defined on a SubDomain, e.g. p(ix), only substituted the derivative Dimension in the stencil. Functions indexed by the sibling Dimension stayed unshifted, so (f + p).dx silently dropped p and (f*p).dx never shifted p. Differentiable.shift now shifts every space Dimension sharing the root of the given Dimension, and make_derivative builds the stencil points with it, in both the expanded and unexpanded paths. Tests cover Add and Mul with both expansion modes, and a one-face CPML forward/adjoint dot test whose adjoint term mixes Grid and SubDomain Functions inside .dx.T.

CIRE builds aliases by translating Indexeds. An unbound StencilDimension appearing outside of any Indexed, e.g. `i0` in `x + i0` within the body of an IndexDerivative, is not translated, so the alias would be evaluated at the wrong stencil point. For example, the invariants pass hoisted such an expression into an array over `x`, dropping `i0`. `_do_generate` now leaves the terms reading such StencilDimensions out of the alias, or discards the candidate altogether.

…djoints `expr.dx(halo=0)` treats `expr` as zero outside the SubDomain of its equation, and `.T` keeps the flag. Hence `Inc(q_bar, out_bar.dx(halo=0).T, subdomain=S)` is the adjoint of `Eq(out, q.dx, subdomain=S)`, i.e. D^T R^T with R the restriction to S, without a user-side zero-padded work field. Equations without `halo=0` derivatives, or without a SubDomain, are evaluated as before. - An equation evaluates its rhs with `_eval_at(lhs, subdomain=S)`. A `halo=0` derivative records S, and multiplies each stencil tap by a branch-free integer MIN/MAX indicator of the grid point it reads being in S (SubDomain.indicator), which invariant hoisting computes once. Along the Dimensions it does not differentiate, it restricts its argument to S at the evaluation point. In a sum with such derivatives, the other terms are restricted to S at the evaluation point (SubDomain.restrict); factors are not. - The evaluated derivative records its stencil radius (halo_radius), and the equation iterates over S grown by the largest one (SubDomain.grow, a SubDomain with a parent). Tensor equations aggregate their components. - It works for expanded, unexpanded and staggered derivatives. Derivatives of `halo=0` derivatives, MultiSubDomains and methods other than FD raise NotImplementedError. Tests: dot tests of `Eq(out, g.dx, subdomain=S)` against `Inc(g_bar, out_bar.dx(halo=0).T, subdomain=S)` for left/right/middle SubDomains, first and second derivatives, orders 2/4/8, Eq and Inc, 2D and unexpanded forms, in 1D and 3D; staggered dot tests against `-out_bar.dx(halo=0)`, as in elastic adjoints, including a staggered vector equation; sums, factors and vector equations; SubDomain growth; a forward check; and a one-face CPML operator.

mloubout added API api (symbolics, types, ...) feature-request labels Sep 28, 2026

mloubout force-pushed the fix-mixed-subdomain-derivative branch 2 times, most recently from 3ae0cf3 to a56e2c6 Compare September 29, 2026 01:11

mloubout marked this pull request as ready for review September 29, 2026 01:37

mloubout force-pushed the fix-mixed-subdomain-derivative branch 2 times, most recently from 0707079 to 384bb84 Compare September 29, 2026 03:51

FabioLuporini reviewed Sep 29, 2026

View reviewed changes

EdCaunt reviewed Sep 29, 2026

View reviewed changes

mloubout force-pushed the fix-mixed-subdomain-derivative branch 9 times, most recently from 1682f1c to b02b62e Compare September 30, 2026 04:38

FabioLuporini reviewed Sep 30, 2026

View reviewed changes

EdCaunt reviewed Sep 30, 2026

View reviewed changes

mloubout force-pushed the fix-mixed-subdomain-derivative branch 2 times, most recently from d2f12f7 to d581d3f Compare October 1, 2026 14:10

FabioLuporini reviewed Oct 1, 2026

View reviewed changes

mloubout added 3 commits October 1, 2026 11:48

mloubout force-pushed the fix-mixed-subdomain-derivative branch from d581d3f to 6573324 Compare October 1, 2026 15:48

		Max(0, Min(1, upper - point + 1)))


		class GrownSubDomain(SubDomain):

		Branch-free indicator of the point `offset` points away from the current
		point along `dim` lying within this SubDomain: 1 inside, 0 outside.


		terms = cbk_compose(i)

		# Aliases only translate Indexeds, so a term reading an unbound

Conversation

mloubout commented Sep 28, 2026

Uh oh!

codecov Bot commented Sep 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

codecov Bot commented Sep 28, 2026 •

edited

Loading