From fe8961cf6d6d9b5e49d25f40effceb3076bf3d9b Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 12 Aug 2026 15:21:43 -0400 Subject: [PATCH 1/7] dsl: fix processing of tensor args --- devito/types/basic.py | 23 ++++++++++++++++------- devito/types/tensor.py | 9 ++++++--- tests/test_tensors.py | 20 ++++++++++++++++++++ 3 files changed, 42 insertions(+), 10 deletions(-) diff --git a/devito/types/basic.py b/devito/types/basic.py index 81d20bb670..3f8e0f6ec0 100644 --- a/devito/types/basic.py +++ b/devito/types/basic.py @@ -4,6 +4,7 @@ from contextlib import contextmanager, suppress from ctypes import POINTER, Structure, _Pointer, c_char, c_char_p from functools import cached_property, reduce +from numbers import Number from operator import mul import numpy as np @@ -1539,7 +1540,10 @@ def _new(cls, *args, **kwargs): # Filter grid and dimensions grid, dimensions = newobj._infer_dims() if grid is None and dimensions is None: - return sympy.ImmutableDenseMatrix(*args) + # Downgrade to a plain Matrix, reusing the representation rather + # than rebuilding from `args`, as the latter would sympify the + # entries with sympy's `sympify` instead of `cls._sympify` + return sympy.ImmutableDenseMatrix._fromrep(newobj._rep) # Initialized with constructed object newobj.__init_finalize__(newobj.rows, newobj.cols, newobj.flat(), grid=grid, dimensions=dimensions) @@ -1581,13 +1585,18 @@ def __subfunc_setup__(cls, *args, **kwargs): @classmethod def _sympify(cls, arg): # This is used internally by sympy to process arguments at rebuilt. And since - # some of our properties are non-sympyfiable we need to have a fallback. - # `strict` so that strings are left alone rather than parsed into Symbols, - # while plain numbers are turned into `Expr` as sympy expects (a Matrix - # holding non-`Expr` entries, such as a plain `int` 0, is deprecated) + # some of our properties are non-sympyfiable we need to have a fallback + if isinstance(arg, Number): + # Plain numbers must be sympified, as sympy assigns the `EXRAW` domain + # to a Matrix holding non-`Expr` entries such as a plain `int` 0 + return sympy.sympify(arg) try: - return sympy.sympify(arg, strict=True) - except sympy.SympifyError: + # Pure sympy object + return arg._sympy_() + except AttributeError: + # Anything else, such as a `Staggering`, is passed through untouched. + # Note that sympifying is not an option here, as it would convert + # away the type, `Staggering` being a `tuple` for example return arg @classmethod diff --git a/devito/types/tensor.py b/devito/types/tensor.py index 3e5586eb32..08cb95fe63 100644 --- a/devito/types/tensor.py +++ b/devito/types/tensor.py @@ -24,12 +24,15 @@ def staggering(stagg, i, j, d, dims): if stagg is None: # No input return NODE if i == j else (d, dims[j]) + elif isinstance(stagg, MatrixBase): + # From rebuild/tensor property. Indexed as a sympy Matrix. Note that this + # may be a plain Matrix rather than an AbstractTensor, as rebuilding a + # tensor component-wise downgrades it when the components aren't Devito + # objects, which is the case for a Matrix of `Staggering` + return stagg[i, j] elif isinstance(stagg, (tuple, list)): # User input as list or tuple return stagg[i][j] - elif isinstance(stagg, AbstractTensor): - # From rebuild/tensor property. Indexed as a sympy Matrix - return stagg[i, j] class TensorFunction(AbstractTensor): diff --git a/tests/test_tensors.py b/tests/test_tensors.py index 53725a493e..0b4ad4603a 100644 --- a/tests/test_tensors.py +++ b/tests/test_tensors.py @@ -12,6 +12,7 @@ ) from devito.symbolics import retrieve_derivatives from devito.types import NODE +from devito.types.utils import Staggering def dimify(dimensions): @@ -528,6 +529,25 @@ def test_diag_sympified_zeros(func1): assert all(isinstance(c, sympy.Expr) for c in f2.flat()) +@pytest.mark.parametrize('func1', [TensorFunction, TensorTimeFunction, + VectorFunction, VectorTimeFunction]) +def test_staggered_attribute_roundtrip(func1): + """ + Accessing an attribute rebuilds the tensor component-wise, which must not + sympify a `Staggering` away, otherwise it can no longer be fed back as the + `staggered` kwarg. + """ + grid = Grid(tuple([5]*3)) + f1 = func1(name="f1", grid=grid, time_order=1) + + stagg = f1.staggered + assert all(isinstance(s, Staggering) for s in stagg.flat()) + + f2 = func1(name="f2", grid=grid, time_order=1, staggered=stagg) + assert all(c1.staggered == c2.staggered + for c1, c2 in zip(f1.flat(), f2.flat(), strict=True)) + + def test_non_expr_components(): """ A tensor may legitimately hold non-`Expr` components, which sympy deprecates From 4f52a11850ecce6ea26edab869c3264393f05d94 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 12 Aug 2026 21:20:58 -0400 Subject: [PATCH 2/7] dsl: Tabulate sparse interpolation in the runtime grid frame An Operator compiled against one model and applied to another overrides every symbol, sparse ones included, and the kernel then runs in the override's frame. `_arg_values` took the coordinates from the override but the origin from `self`, tabulating real coordinates against the compile-time frame. Every point landed off by the difference between the two origins, silently, since sources and receivers shift together. --- devito/types/sparse.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/devito/types/sparse.py b/devito/types/sparse.py index d5952ef750..209dbf6991 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -1020,13 +1020,22 @@ def _arg_values(self, estimate_memory=False, **kwargs): # Resolve the runtime grid origin (honours `o_x`/`o_y`/... overrides) # and hand it to the interpolator so tables reflect the actual frame - # of reference used by the kernel. - onames = [o.name for o in self.grid.origin_symbols] + # of reference used by the kernel. A runtime override may carry its + # own Grid -- an Operator compiled on one model and applied to another + # overrides every symbol, ours included -- and the kernel then runs in + # *that* frame. The coordinates below already come from the override, + # so the origin has to as well; reading it off `self` tabulates real + # coordinates against the compile-time frame and places every point + # off by the difference between the two origins. + key = kwargs.get(self.name, self) + if not isinstance(key, AbstractSparseFunction): + key = self + onames = [o.name for o in key.grid.origin_symbols] 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) + zip(onames, key.grid.origin, strict=True)) + coords = values.get(self.coordinates.name, key.coordinates.data) values.update(self.interpolator._arg_defaults( - coords=coords, sfunc=self, origin=origin + coords=coords, sfunc=key, origin=origin )) return values From ccf9b73e4e0e2c06ddb5e5bda062d8f3c51c94f4 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 12 Aug 2026 21:20:59 -0400 Subject: [PATCH 3/7] tests: Cover sparse interpolation under a runtime grid override Injects through an Operator built on a grid whose origin differs from the one it is applied to. The linear path regressed here; sinc computes its positions from the runtime `o_x` and did not. --- tests/test_interpolation.py | 47 +++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index fdd69fc1d7..9e9bbf9cae 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -1510,3 +1510,50 @@ def test_inject_subdomain_mpi(self, mode, grid, coords): assert data1 == None # noqa assert data2 == None # noqa assert data3 == None # noqa + + +class TestRuntimeGridOverride: + """ + A runtime override may carry its own Grid: an Operator compiled against one + model and applied to another overrides every symbol, sparse ones included. + Tabulated interpolation data must then be built in the override's frame of + reference, not the compile-time one. + """ + + @pytest.mark.parametrize('interpolation,r', [('linear', 1), ('sinc', 4)]) + def test_inject_honors_override_origin(self, interpolation, r): + """ + Inject at a fixed physical coordinate through an Operator built on a + grid whose origin differs from the one it is applied to. The point + must land at the index the *runtime* origin implies. + + Reading the origin off the compile-time symbol instead tabulates real + coordinates against the wrong frame and displaces every point by the + difference between the two origins -- silently, since sources and + receivers shift together. + """ + shape, spacing, coord = (41, 41), (10., 10.), 120. + extent = tuple((s - 1) * h for s, h + in zip(shape, spacing, strict=True)) + kw = dict(interpolation=interpolation, r=r) + + def setup(origin): + grid = Grid(shape=shape, extent=extent, origin=origin) + u = TimeFunction(name='u', grid=grid, space_order=8) + src = SparseTimeFunction(name='src', grid=grid, npoint=1, nt=2, **kw) + src.coordinates.data[0, :] = coord + src.data[:] = 1. + return u, src + + # Compile against a grid whose origin is the physical one ... + u_build, src_build = setup((0., 0.)) + op = Operator(src_build.inject(field=u_build.forward, expr=src_build)) + + # ... and apply it to one shifted by a full absorbing layer. + shift = -100. + u, src = setup((shift, shift)) + op.apply(time_M=0, u=u, src=src) + + expected = tuple(int((coord - shift) / h) for h in spacing) + peak = np.unravel_index(np.argmax(np.abs(u.data)), u.data.shape)[1:] + assert peak == expected From 50dd4953a2f9174e5fe7cf49f9d1cd31d45a35e7 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 12 Aug 2026 22:30:02 -0400 Subject: [PATCH 4/7] dsl: Resolve the sparse interpolation frame through the alias key `_arg_defaults` is already handed the runtime object as `self` and the compile-time symbols as `alias`, like the rest of the `_arg_*` protocol, so the frame to tabulate in was there all along. Take the names from `alias` and the data and Grid from `self`, rather than looking the override up in kwargs from `_arg_values`. `_arg_values` now only rebuilds the tables for an explicit `o_x`/`o_y` override, instead of redoing on every call what `_arg_defaults` has already tabulated correctly. --- devito/types/sparse.py | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/devito/types/sparse.py b/devito/types/sparse.py index 209dbf6991..f74276f91c 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -1007,10 +1007,16 @@ def _arg_defaults(self, alias=None, estimate_memory=False): defaults = super()._arg_defaults(alias=alias, estimate_memory=estimate_memory) if estimate_memory: return defaults + # `alias` names the symbols the Operator was compiled with, while + # `self` carries the data and the Grid the kernel actually runs in. + # The two differ when an Operator built against one model is applied + # to another, and tabulated data belongs to the latter's frame of + # reference: tabulating `self`'s coordinates against `alias`'s origin + # displaces every point by the difference between the two. key = alias or self - coords = defaults.get(key.coordinates.name, key.coordinates.data) + coords = defaults.get(key.coordinates.name, self.coordinates.data) defaults.update(key.interpolator._arg_defaults(coords=coords, - sfunc=key)) + sfunc=self)) return defaults def _arg_values(self, estimate_memory=False, **kwargs): @@ -1018,24 +1024,19 @@ def _arg_values(self, estimate_memory=False, **kwargs): if estimate_memory: return values - # Resolve the runtime grid origin (honours `o_x`/`o_y`/... overrides) - # and hand it to the interpolator so tables reflect the actual frame - # of reference used by the kernel. A runtime override may carry its - # own Grid -- an Operator compiled on one model and applied to another - # overrides every symbol, ours included -- and the kernel then runs in - # *that* frame. The coordinates below already come from the override, - # so the origin has to as well; reading it off `self` tabulates real - # coordinates against the compile-time frame and places every point - # off by the difference between the two origins. - key = kwargs.get(self.name, self) - if not isinstance(key, AbstractSparseFunction): - key = self - onames = [o.name for o in key.grid.origin_symbols] + # `super` has already tabulated through `_arg_defaults`, in the frame + # of whichever object supplied the runtime values. Only an explicit + # `o_x`/`o_y`/... override moves that frame again, and the tables then + # have to be rebuilt against it. + onames = [o.name for o in self.grid.origin_symbols] + if not any(n in kwargs for n in onames): + return values + origin = tuple(kwargs.get(n, o) for n, o in - zip(onames, key.grid.origin, strict=True)) - coords = values.get(self.coordinates.name, key.coordinates.data) + 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=key, origin=origin + coords=coords, sfunc=self, origin=origin )) return values From 3101801595449212bee7890d92cd329630ca5e6f Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 12 Aug 2026 22:31:52 -0400 Subject: [PATCH 5/7] misc: Drop the redundant alias comment in sparse _arg_defaults --- devito/types/sparse.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/devito/types/sparse.py b/devito/types/sparse.py index f74276f91c..6e913646a3 100644 --- a/devito/types/sparse.py +++ b/devito/types/sparse.py @@ -1007,12 +1007,6 @@ def _arg_defaults(self, alias=None, estimate_memory=False): defaults = super()._arg_defaults(alias=alias, estimate_memory=estimate_memory) if estimate_memory: return defaults - # `alias` names the symbols the Operator was compiled with, while - # `self` carries the data and the Grid the kernel actually runs in. - # The two differ when an Operator built against one model is applied - # to another, and tabulated data belongs to the latter's frame of - # reference: tabulating `self`'s coordinates against `alias`'s origin - # displaces every point by the difference between the two. key = alias or self coords = defaults.get(key.coordinates.name, self.coordinates.data) defaults.update(key.interpolator._arg_defaults(coords=coords, From d65dab177a35ddc76869c5024bf797a811ef6c81 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 12 Aug 2026 22:47:58 -0400 Subject: [PATCH 6/7] tests: Fold the runtime grid override case into TestCustom It sits next to `test_position`, which covers the same origin shift through an explicit `o_x`, rather than carrying a class of its own. --- tests/test_interpolation.py | 78 +++++++++++++++---------------------- 1 file changed, 31 insertions(+), 47 deletions(-) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 9e9bbf9cae..635857bdea 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -1077,6 +1077,37 @@ def test_position(self, shape): assert(np.allclose(rec.data, rec1.data, atol=1e-5)) + @pytest.mark.parametrize('interpolation,r', [('linear', 1), ('sinc', 4)]) + def test_position_override_grid(self, interpolation, r): + """ + Inject through an Operator built on a grid whose origin differs from + the one it is applied to, as when an Operator compiled against one + model is applied to another. The point must land where the runtime + origin puts it, not the compile-time one. + """ + shape, spacing, coord = (41, 41), (10., 10.), 120. + extent = tuple((s - 1) * h for s, h in zip(shape, spacing, strict=True)) + kw = dict(interpolation=interpolation, r=r) + + def setup(origin): + grid = Grid(shape=shape, extent=extent, origin=origin) + u = TimeFunction(name='u', grid=grid, space_order=8) + src = SparseTimeFunction(name='src', grid=grid, npoint=1, nt=2, **kw) + src.coordinates.data[0, :] = coord + src.data[:] = 1. + return u, src + + u_build, src_build = setup((0., 0.)) + op = Operator(src_build.inject(field=u_build.forward, expr=src_build)) + + shift = -100. + u, src = setup((shift, shift)) + op.apply(time_M=0, u=u, src=src) + + expected = tuple(int((coord - shift) / h) for h in spacing) + peak = np.unravel_index(np.argmax(np.abs(u.data)), u.data.shape)[1:] + assert peak == expected + def test_sparse_first(self): """ Tests custom sprase function with sparse dimension as first index. @@ -1510,50 +1541,3 @@ def test_inject_subdomain_mpi(self, mode, grid, coords): assert data1 == None # noqa assert data2 == None # noqa assert data3 == None # noqa - - -class TestRuntimeGridOverride: - """ - A runtime override may carry its own Grid: an Operator compiled against one - model and applied to another overrides every symbol, sparse ones included. - Tabulated interpolation data must then be built in the override's frame of - reference, not the compile-time one. - """ - - @pytest.mark.parametrize('interpolation,r', [('linear', 1), ('sinc', 4)]) - def test_inject_honors_override_origin(self, interpolation, r): - """ - Inject at a fixed physical coordinate through an Operator built on a - grid whose origin differs from the one it is applied to. The point - must land at the index the *runtime* origin implies. - - Reading the origin off the compile-time symbol instead tabulates real - coordinates against the wrong frame and displaces every point by the - difference between the two origins -- silently, since sources and - receivers shift together. - """ - shape, spacing, coord = (41, 41), (10., 10.), 120. - extent = tuple((s - 1) * h for s, h - in zip(shape, spacing, strict=True)) - kw = dict(interpolation=interpolation, r=r) - - def setup(origin): - grid = Grid(shape=shape, extent=extent, origin=origin) - u = TimeFunction(name='u', grid=grid, space_order=8) - src = SparseTimeFunction(name='src', grid=grid, npoint=1, nt=2, **kw) - src.coordinates.data[0, :] = coord - src.data[:] = 1. - return u, src - - # Compile against a grid whose origin is the physical one ... - u_build, src_build = setup((0., 0.)) - op = Operator(src_build.inject(field=u_build.forward, expr=src_build)) - - # ... and apply it to one shifted by a full absorbing layer. - shift = -100. - u, src = setup((shift, shift)) - op.apply(time_M=0, u=u, src=src) - - expected = tuple(int((coord - shift) / h) for h in spacing) - peak = np.unravel_index(np.argmax(np.abs(u.data)), u.data.shape)[1:] - assert peak == expected From 8f120211ec25dcfa5e2c5e8847550e04ec7d82fb Mon Sep 17 00:00:00 2001 From: mloubout Date: Thu, 13 Aug 2026 07:30:20 -0400 Subject: [PATCH 7/7] examples: make rotated fd accept input b --- examples/seismic/tti/operators.py | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/examples/seismic/tti/operators.py b/examples/seismic/tti/operators.py index f6ec8d34bb..110ecb2c86 100644 --- a/examples/seismic/tti/operators.py +++ b/examples/seismic/tti/operators.py @@ -62,7 +62,7 @@ def trig_func(model): return costheta, sintheta -def Gzz_centered(model, field): +def Gzz_centered(model, field, b=None): """ 3D rotated second order derivative in the direction z. @@ -72,12 +72,17 @@ def Gzz_centered(model, field): Physical parameters model structure. field : Function Input for which the derivative is computed. + b : Function, optional + Buoyancy to build the operator with, defaulting to the model's. Since + the operator is linear in it, passing a perturbation here gives the + derivative of the operator with respect to the buoyancy in that + direction. Returns ------- Rotated second order derivative w.r.t. z. """ - b = getattr(model, 'b', 1) + b = getattr(model, 'b', 1) if b is None else b costheta, sintheta, cosphi, sinphi = trig_func(model) order1 = field.space_order // 2 @@ -99,7 +104,7 @@ def Gzz_centered(model, field): return Gzz -def Gzz_centered_2d(model, field): +def Gzz_centered_2d(model, field, b=None): """ 2D rotated second order derivative in the direction z. @@ -109,12 +114,17 @@ def Gzz_centered_2d(model, field): Physical parameters model structure. field : Function Input for which the derivative is computed. + b : Function, optional + Buoyancy to build the operator with, defaulting to the model's. Since + the operator is linear in it, passing a perturbation here gives the + derivative of the operator with respect to the buoyancy in that + direction. Returns ------- Rotated second order derivative w.r.t. z. """ - b = getattr(model, 'b', 1) + b = getattr(model, 'b', 1) if b is None else b costheta, sintheta = trig_func(model) order1 = field.space_order // 2 @@ -133,7 +143,7 @@ def Gzz_centered_2d(model, field): # Centered case produces directly Gxx + Gyy -def Gh_centered(model, field): +def Gh_centered(model, field, b=None): """ Sum of the 3D rotated second order derivative in the direction x and y. As the Laplacian is rotation invariant, it is computed as the conventional @@ -146,13 +156,19 @@ def Gh_centered(model, field): Physical parameters model structure. field : Function Input field. + b : Function, optional + Buoyancy to build the operator with, defaulting to the model's. See + :func:`Gzz_centered`. Returns ------- Sum of the 3D rotated second order derivative in the direction x and y. """ - Gzz = Gzz_centered(model, field) if model.dim == 3 else Gzz_centered_2d(model, field) - b = getattr(model, 'b', None) + b = getattr(model, 'b', None) if b is None else b + if model.dim == 3: # noqa: SIM108 + Gzz = Gzz_centered(model, field, b=b) + else: + Gzz = Gzz_centered_2d(model, field, b=b) if b is not None: _diff = lambda f, d: getattr(f, f'd{d.name}') so = field.space_order // 2