From d0eba846eb919381a024eadacb288f4a304c93b9 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Tue, 3 Nov 2026 17:10:49 +0000 Subject: [PATCH 1/6] dsl: Fix up staggerd derivative evaluation for functions on subdomains --- devito/finite_differences/tools.py | 8 +++++-- tests/test_subdomains.py | 37 ++++++++++++++++++++++++++++-- 2 files changed, 41 insertions(+), 4 deletions(-) diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py index 4d816866fb..b7fe56f859 100644 --- a/devito/finite_differences/tools.py +++ b/devito/finite_differences/tools.py @@ -299,8 +299,12 @@ def generate_indices(expr, dim, order, side=None, matvec=None, x0=None, nweights iexpr = x0 + d * dim.spacing return IndexSet(dim, expr=iexpr), x0 - # Evaluation point relative to the expression's grid - mid = (x0 - expr.indices_ref[dim]).subs({dim: 0, dim.spacing: 1}) + # Evaluation point relative to the expression's grid, which may be indexed by a + # SubDimension of `dim` (e.g., a Function on a SubDomain) + rel = x0 - expr.indices_ref[dim] + mid = rel.subs({**{i: 0 for i in rel.free_symbols + if i.is_Dimension and i.root is dim.root}, + dim.spacing: 1}) # Shift for side side = side or centered diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py index eff6758636..48c9746448 100644 --- a/tests/test_subdomains.py +++ b/tests/test_subdomains.py @@ -6,8 +6,8 @@ from conftest import assert_structure, opts_tiling from devito import ( - Border, Buffer, ConditionalDimension, Constant, Eq, Function, Grid, Lt, Operator, - SparseFunction, SparseTimeFunction, SubDimension, SubDomain, SubDomainSet, + Border, Buffer, ConditionalDimension, Constant, Derivative, Eq, Function, Grid, Lt, + Operator, SparseFunction, SparseTimeFunction, SubDimension, SubDomain, SubDomainSet, TensorFunction, TimeFunction, VectorFunction, solve ) from devito.exceptions import InvalidArgument @@ -1601,6 +1601,39 @@ def define(self, dimensions): assert np.all(np.isclose(fdx.data[:], gdx.data[2:-2, 3:-1])) assert np.all(np.isclose(fdy.data[:], gdy.data[2:-2, 3:-1])) + @pytest.mark.parametrize('s_o', [2, 4]) + @pytest.mark.parametrize('axis', [0, 1]) + def test_derivatives_staggered(self, s_o, axis): + """ + Derivatives of Functions on SubDomains can be evaluated at points staggered + along the derivative Dimension + """ + + class Middle(SubDomain): + + name = 'middle' + + def define(self, dimensions): + x, y = dimensions + return {x: ('middle', 2, 2), y: ('middle', 3, 1)} + + grid = Grid(shape=(10, 10), extent=(9., 9.)) + d = grid.dimensions[axis] + mid = Middle(grid=grid) + + f = Function(name='f', grid=mid, space_order=s_o) + g = Function(name='g', grid=grid, space_order=s_o) + fd, gd = (Function(name=i, grid=grid, staggered=d) for i in ('fd', 'gd')) + + f.data[:] = np.random.default_rng(0).random(f.shape) + g.data[2:-2, 3:-1] = f.data + + Operator([Eq(fd, Derivative(f, d), subdomain=mid), + Eq(gd, Derivative(g, d), subdomain=mid)])() + + assert np.any(fd.data != 0) + assert np.allclose(fd.data, gd.data) + def test_operator_shortcuts(self): """ Test that shortcuts for Laplacian, divergence, etc all work with Functions From d5e87e37ebdcad5729cbc883eb0d5811e15ed57c Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Tue, 3 Nov 2026 17:34:13 +0000 Subject: [PATCH 2/6] dsl: Fixup halo=0 with cross-derivatives --- devito/finite_differences/finite_difference.py | 10 +++++++--- tests/test_derivatives.py | 2 +- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py index bb81e712a9..0f94304f97 100644 --- a/devito/finite_differences/finite_difference.py +++ b/devito/finite_differences/finite_difference.py @@ -6,7 +6,7 @@ from devito.logger import warning -from .differentiable import DiffDerivative, EvalDerivative, Weights +from .differentiable import DiffDerivative, EvalDerivative, Weights, merge_growth from .tools import ( centered, check_input, direct, fd_weights_registry, generate_indices, left, process_weights, right, transpose @@ -236,8 +236,12 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, coeffici # With a `subdomain`, `expr` is treated as zero outside it: every stencil tap # is multiplied by the indicator of the point it reads, 1 if `subdomain` # spans the whole of `dim`. The derivative then extends past `subdomain` by - # the stencil radius - growth = {dim.root: indices.radius} if subdomain is not None else {} + # the stencil radius, on top of the growth of `expr` itself (e.g., the inner + # derivative of a cross derivative) + if subdomain is not None: + growth = merge_growth([expr._growth, {dim.root: indices.radius}]) + else: + growth = {} if not expand and indices.expr is not None: weights = Weights(name='w', dimensions=indices.free_dim, diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py index ee93c00488..f15dfe78b0 100644 --- a/tests/test_derivatives.py +++ b/tests/test_derivatives.py @@ -1861,7 +1861,7 @@ def test_adjoint_unexpanded_3d(self): assert self.dot_test(grid, subdomain, 'dy', 8, opt=opt) < 1e-12 - @pytest.mark.parametrize('deriv', ['dx', 'dy']) + @pytest.mark.parametrize('deriv', ['dx', 'dy', 'dxdy']) def test_adjoint_2d(self, deriv): grid = Grid(shape=(20, 22), extent=(19., 21.), dtype=np.float64) subdomain = Region('s', {'x': ('left', 6), 'y': ('middle', 4, 5)}, grid=grid) From 27d703504a1e065d092ca93c7d2e07754896e147 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Wed, 4 Nov 2026 16:14:57 +0000 Subject: [PATCH 3/6] dsl: Normalise staggers and fixup test which passed due to stagger of function on subdomain being dropped --- devito/types/dense.py | 8 +++++--- tests/test_subdomains.py | 6 +++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/devito/types/dense.py b/devito/types/dense.py index 30f543a47c..6e18f2e382 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -1165,12 +1165,14 @@ def __staggered_setup__(cls, dimensions, staggered=None, **kwargs): processed = staggered else: # Staggering is not NODE or CELL or None - # therefore it's a tuple of dimensions + # therefore it's a tuple of dimensions. These may be the roots of + # `dimensions` (e.g., `staggered=x` for a Function on a SubDomain) + staggered = set(as_tuple(staggered)) processed = [] for d in dimensions: - if d in as_tuple(staggered): + if {d, d.root} & staggered: processed.append(sympy.S.One) - elif -d in as_tuple(staggered): + elif {-d, -d.root} & staggered: processed.append(sympy.S.NegativeOne) else: processed.append(sympy.S.Zero) diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py index 48c9746448..61ff2199fb 100644 --- a/tests/test_subdomains.py +++ b/tests/test_subdomains.py @@ -2015,14 +2015,14 @@ def define(self, dimensions): return {x: ('left', 2), y: y, z: z} sd1 = SD1(grid=grid) - f = TimeFunction(name='f', grid=grid, time_order=1, space_order=2, - staggered=(x,)) + f = TimeFunction(name='f', grid=grid, time_order=1, space_order=2) g = TimeFunction(name='g', grid=sd1, time_order=1, space_order=2, staggered=(x,)) eq = Eq(g, g + f.dx) eqe = eq.evaluate - assert eqe.rhs == g + f.dx(x0=x).evaluate._subs(x, g.dimensions[1]) + x0 = x + x.spacing/2 + assert eqe.rhs == g + f.dx(x0=x0).evaluate._subs(x, g.dimensions[1]) class TestSubDomainArguments: From ce049364860066a13de6fbcf2583e72df2e38ef3 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Thu, 5 Nov 2026 11:03:42 +0000 Subject: [PATCH 4/6] dsl: Express x0 in the same Dimension as the indices --- devito/finite_differences/tools.py | 4 ++++ tests/test_subdomains.py | 3 ++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py index b7fe56f859..f383702a50 100644 --- a/devito/finite_differences/tools.py +++ b/devito/finite_differences/tools.py @@ -323,6 +323,10 @@ def generate_indices(expr, dim, order, side=None, matvec=None, x0=None, nweights d = make_stencil_dimension(expr, o_min, o_max) iexpr = expr.indices_ref[dim] + d * dim.spacing + # Express `x0` in the same Dimension as the indices, or the FD weights retain + # e.g. `iy - y` and are not reduced to numbers + x0 = expr.indices_ref[dim] + mid * dim.spacing + return IndexSet(dim, expr=iexpr), x0 diff --git a/tests/test_subdomains.py b/tests/test_subdomains.py index 61ff2199fb..bcde422cba 100644 --- a/tests/test_subdomains.py +++ b/tests/test_subdomains.py @@ -1631,8 +1631,9 @@ def define(self, dimensions): Operator([Eq(fd, Derivative(f, d), subdomain=mid), Eq(gd, Derivative(g, d), subdomain=mid)])() + # Same stencil, so bitwise identical assert np.any(fd.data != 0) - assert np.allclose(fd.data, gd.data) + assert np.all(fd.data == gd.data) def test_operator_shortcuts(self): """ From d3318573b83e8c7fd82596893f92c2e08ae3fc71 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Thu, 5 Nov 2026 11:28:11 +0000 Subject: [PATCH 5/6] misc: Key staggering DimensionTuple off root dimensions --- devito/types/dense.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/devito/types/dense.py b/devito/types/dense.py index 6e18f2e382..fbfcd2dac2 100644 --- a/devito/types/dense.py +++ b/devito/types/dense.py @@ -1176,7 +1176,7 @@ def __staggered_setup__(cls, dimensions, staggered=None, **kwargs): processed.append(sympy.S.NegativeOne) else: processed.append(sympy.S.Zero) - return Staggering(*processed, getters=dimensions) + return Staggering(*processed, getters=[d.root for d in dimensions]) @classmethod def __indices_setup__(cls, *args, **kwargs): From c804e064d3e10951dd06583dcf64811f4bea7672 Mon Sep 17 00:00:00 2001 From: Edward Caunt Date: Fri, 6 Nov 2026 15:22:50 +0000 Subject: [PATCH 6/6] dsl: Normalise x0 --- devito/finite_differences/tools.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py index f383702a50..f517e4fd25 100644 --- a/devito/finite_differences/tools.py +++ b/devito/finite_differences/tools.py @@ -299,12 +299,12 @@ def generate_indices(expr, dim, order, side=None, matvec=None, x0=None, nweights iexpr = x0 + d * dim.spacing return IndexSet(dim, expr=iexpr), x0 - # Evaluation point relative to the expression's grid, which may be indexed by a - # SubDimension of `dim` (e.g., a Function on a SubDomain) - rel = x0 - expr.indices_ref[dim] - mid = rel.subs({**{i: 0 for i in rel.free_symbols - if i.is_Dimension and i.root is dim.root}, - dim.spacing: 1}) + # Normalize x0 `Dimension` to match index if the latter is a `SubDimension` of + # the former + x0 = x0.subs({d.root: d for d in expr.dimensions if d.root is dim and d.is_Sub}) + + # Evaluation point relative to the expression's grid + mid = (x0 - expr.indices_ref[dim]).subs({dim: 0, dim.spacing: 1}) # Shift for side side = side or centered @@ -323,10 +323,6 @@ def generate_indices(expr, dim, order, side=None, matvec=None, x0=None, nweights d = make_stencil_dimension(expr, o_min, o_max) iexpr = expr.indices_ref[dim] + d * dim.spacing - # Express `x0` in the same Dimension as the indices, or the FD weights retain - # e.g. `iy - y` and are not reduced to numbers - x0 = expr.indices_ref[dim] + mid * dim.spacing - return IndexSet(dim, expr=iexpr), x0