From a452ad2b00ae4d190abfddf5fb589f5034d7b9d0 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 11 Aug 2026 13:17:42 -0400 Subject: [PATCH 1/6] dsl: Let an Eq carry its own interp mode The symmetric placement costs an extra inline evaluation of the block, so it is wanted on the handful of equations that need it rather than on a whole operator. `Eq` takes an `interp_mode` and passes it down when it evaluates, leaving `sym_opt` as the default for everything else. --- devito/types/equation.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/devito/types/equation.py b/devito/types/equation.py index 5b709ee89e..42ddb07dae 100644 --- a/devito/types/equation.py +++ b/devito/types/equation.py @@ -36,6 +36,15 @@ class Eq(sympy.Eq, Evaluable, Pickable): An ordered list of Dimensions that do not explicitly appear in either the left-hand side or in the right-hand side, but that should be honored when constructing an Operator. + interp_mode : str, optional, default=None + Overrides the Operator's `sym_opt={'interp-mode': ...}` for this + equation only. An Operator generally wants one mode -- `'direct'`, which + keeps finite-difference stencils compact -- while a single equation in + it may need `'symmetric'`, which re-associates a product of staggered + operands so that the discretized operator is the transpose of itself. + A gradient accumulated alongside the propagation it is the adjoint of is + the typical case: the stencils must stay compact, but the accumulation + has to be an exact transpose or the gradient degrades to first order. Examples -------- @@ -60,10 +69,10 @@ class Eq(sympy.Eq, Evaluable, Pickable): is_Reduction = False __rargs__ = ('lhs', 'rhs') - __rkwargs__ = ('subdomain', 'coefficients', 'implicit_dims') + __rkwargs__ = ('subdomain', 'coefficients', 'implicit_dims', 'interp_mode') def __new__(cls, lhs, rhs=0, subdomain=None, coefficients=None, - implicit_dims=None, **kwargs): + implicit_dims=None, interp_mode=None, **kwargs): if coefficients is not None: _ = deprecations.coeff_warn kwargs['evaluate'] = False @@ -76,9 +85,15 @@ def __new__(cls, lhs, rhs=0, subdomain=None, coefficients=None, obj._subdomain = subdomain obj._substitutions = coefficients obj._implicit_dims = as_tuple(implicit_dims) + obj._interp_mode = interp_mode return obj + @property + def interp_mode(self): + """Per-equation override of the Operator's `interp-mode`, or None.""" + return self._interp_mode + @classmethod def _apply_coeffs(cls, expr, coefficients): """ @@ -108,6 +123,8 @@ def _evaluate(self, **kwargs): The RHS of the Equation is evaluated at the indices of the LHS if required. """ + if self._interp_mode is not None: + kwargs['interp_mode'] = self._interp_mode try: lhs = self.lhs._evaluate(**kwargs) rhs = self.rhs._eval_at(self.lhs, **kwargs)._evaluate(**kwargs) @@ -115,7 +132,8 @@ def _evaluate(self, **kwargs): lhs, rhs = self._evaluate_args(**kwargs) eq = self.func(lhs, rhs, subdomain=self.subdomain, coefficients=self.substitutions, - implicit_dims=self._implicit_dims) + implicit_dims=self._implicit_dims, + interp_mode=self._interp_mode) return eq From 9c692ecc27a8bb8bf42367a4ac73b7809ee4e78b Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 11 Aug 2026 13:17:43 -0400 Subject: [PATCH 2/6] dsl: Place a product by where its derivatives land Two changes to what `interp-mode='symmetric'` re-associates. A derivative is where its `x0` puts it, not where the field it differentiates lives. `div(v)` of a staggered velocity lands on the node, and reading the operand's staggering instead made a node-centred kernel look off-node, so it was re-associated onto the operands and its compact stencil replaced by an interpolated one twice as wide. A composite is off only if one of its own operands is, which is what makes a sum of derivatives -- a divergence, a trace -- read as the node-centred quantity it is. The block is then chosen among the factors that are actually off `func`, and it takes two of them for there to be anything to re-associate. With one, the single interpolation `direct` puts on it is already the transpose-consistent form -- it is what makes the `i, j` entry of a stiffness matrix the transpose of its `j, i` entry. With two, `direct` interpolates each separately, and `I(a)*I(b)` is not `I(a*b)`: the discretized operator stops being the transpose of itself, which is invisible in a forward simulation and shows up as a first-order gradient in an adjoint one. The mode also travels down the fallback now. A product that cannot itself be re-associated is routinely wrapped around one that can -- a `dt` scaling, a sum of per-component contractions -- and dropping it there evaluated all of that in `direct`. --- devito/finite_differences/differentiable.py | 94 +++++++++++++++------ 1 file changed, 70 insertions(+), 24 deletions(-) diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index b77d2c158d..89a49be865 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -15,7 +15,9 @@ # Moved in 1.13 from sympy.core.basic import ordering_of_classes -from devito.finite_differences.interpolation import interp_at, post_x0_indices +from devito.finite_differences.interpolation import ( + interp_at, interp_mapper, post_x0_indices +) from devito.finite_differences.tools import coeff_priority, make_shift_x0 from devito.logger import warning from devito.tools import ( @@ -487,8 +489,15 @@ def has_free(self, *patterns): return all(i in self.free_symbols for i in patterns) -def highest_priority(diff_op): - if not diff_op._args_diff: +def highest_priority(diff_op, candidates=None): + """ + The Function whose location a product should be evaluated at. + + `candidates` restricts the choice to a subset of the operands; without it + the whole expression's differentiable arguments are considered. + """ + args_diff = diff_op._args_diff if candidates is None else tuple(candidates) + if not args_diff: return diff_op # We want to get the object with highest priority @@ -496,7 +505,7 @@ def highest_priority(diff_op): # set of dimensions is used when multiple ones with the same # priority appear prio = lambda x: (getattr(x, '_fd_priority', 0), len(x.dimensions)) - prio_func = sorted(diff_op._args_diff, key=prio, reverse=True)[0] + prio_func = sorted(args_diff, key=prio, reverse=True)[0] # The highest priority must be a Function if not isinstance(prio_func, AbstractFunction): @@ -664,6 +673,27 @@ def _gather_for_diff(self): other = self.func(*other)._eval_at(highest_priority(self)) return self.func(other, *derivs) + @classmethod + def _off_func(cls, a, func): + """ + Whether evaluating `a` at `func`'s location takes an interpolation. + + A derivative is where its `x0` puts it, not where the field it + differentiates lives: `div(v)` of a staggered velocity lands on the + node, and reading its operand's staggering instead would re-associate a + product that is already node-centred and interpolate it for nothing. + A composite is off only if one of its own operands is, which is what + makes a sum of derivatives -- a divergence, a trace -- read as the + node-centred quantity it is. + """ + if isinstance(a, sympy.Derivative): + source = post_x0_indices(a, func) + elif isinstance(a, AbstractFunction) or not getattr(a, '_args_diff', ()): + source = a.indices_ref + else: + return any(cls._off_func(i, func) for i in a._args_diff) + return bool(interp_mapper(source, func.indices_ref, a.dimensions)) + def _eval_at(self, func, interp_mode='direct', **kwargs): """ Evaluate a Mul at the location of `func`. @@ -674,36 +704,52 @@ def _eval_at(self, func, interp_mode='direct', **kwargs): independently evaluated at `func`'s location via `Differentiable._eval_at`. - - `interp_mode='symmetric'`: when every Differentiable factor has a - staggering different from `func`'s, apply the `I * (a * I^T * b)` - form: - - 1. Pick a `block` location -- the highest-priority factor's - staggering (NODE is the highest priority, so coefficient-like - NODE factors win, as in the `I * C * I^T` elastic stiffness - pattern). Each factor not at the block is brought there via - `I^T` (an explicit 0-order FD interpolation operator). - Derivatives additionally set `x0` on their own derivative - dimensions to `func`'s indices. + - `interp_mode='symmetric'`: the product is formed *away* from `func` + and closed with a single interpolation, the `I * (a * I^T * b)` form: + + 1. Pick a `block` location -- the highest-priority staggering + among the factors that are not already at `func`'s (NODE is the + highest priority, so coefficient-like NODE factors win, as in + the `I * C * I^T` elastic stiffness pattern). Each factor not at + the block is brought there via `I^T` (an explicit 0-order FD + interpolation operator). Derivatives additionally set `x0` on + their own derivative dimensions to `func`'s indices. 2. The product is formed at `block`'s location. 3. The whole product is interpolated to `func` via `I` (an explicit 0-order FD operator). - When the trigger does not hold (e.g. some factor already matches - `func`'s staggering), we fall back to `direct`. + It takes *two* factors away from `func` to have something to + re-associate. With one, the single interpolation `direct` puts on it + is already the transpose-consistent form -- it is what makes the + `i, j` entry of a stiffness matrix the transpose of its `j, i` entry + -- so we fall back. With two, `direct` would interpolate each of them + separately, and `I(a)*I(b)` is not `I(a*b)`: the discretized operator + stops being the transpose of itself, which is invisible in a forward + simulation and shows up as a first-order gradient in an adjoint one. + + Which factors count is the other half of it. A derivative is where + its `x0` puts it, not where its operand lives, so a product around a + `div(v)` of a staggered velocity -- node-centred however its operands + are staggered -- is left alone rather than re-associated onto the + operands' location, which would replace a compact stencil with an + interpolated one twice as wide. """ if interp_mode != 'symmetric': return super()._eval_at(func, **kwargs) diff, other = split(self.args, lambda a: isinstance(a, Differentiable)) - # Symmetric form requires every Differentiable factor to differ from - # func; otherwise direct evaluation is cleaner and equivalent. - if len(diff) < 2 or \ - any(a.staggered == func.staggered for a in diff): - return super()._eval_at(func, **kwargs) - - block_indices = highest_priority(self).indices_ref + # A single factor cannot be re-associated, and with everything already + # on `func` there is no interpolation to place. The mode still has to + # travel down: a product that cannot itself be re-associated is + # routinely wrapped around one that can (a `dt` scaling, a sum of + # per-component contractions), and dropping the mode here would silently + # evaluate all of that in `direct`. + off_func = [a for a in diff if self._off_func(a, func)] + if len(off_func) < 2: + return super()._eval_at(func, interp_mode=interp_mode, **kwargs) + + block_indices = highest_priority(self, candidates=off_func).indices_ref # Bring each factor to block's location (I^T where needed) new_factors = list(other) From 94089ffbf1a7f099fa89beabd6851f7a26d783a8 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 11 Aug 2026 13:17:44 -0400 Subject: [PATCH 3/6] dsl: Accept the degenerate rebuild of a zero-order derivative A zero-order derivative whose weights collapse to one is the identity, and it comes back as whatever it was applied to. When that is a sum -- the engineering shear strain of a staggered velocity, say -- the rebuild asserted rather than accepting the degenerate result it already handles for a single argument. It is reachable through the symmetric placement of a stiffness contraction. --- devito/finite_differences/differentiable.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py index 89a49be865..88e3cb214e 100644 --- a/devito/finite_differences/differentiable.py +++ b/devito/finite_differences/differentiable.py @@ -1165,9 +1165,11 @@ def __new__(cls, *args, base=None, **kwargs): except AttributeError: # This might happen if e.g. one attempts a (re)construction with # one sole argument. The (re)constructed EvalDerivative degenerates - # to an object of different type, in classic SymPy style. That's fine + # to an object of different type, in classic SymPy style. That's + # fine -- and a single argument that is itself a sum is the same + # 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 - assert not obj.is_Add return obj return obj From eb91c3f60ca4ae4afae441572dfc06dfa523e237 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 11 Aug 2026 13:17:45 -0400 Subject: [PATCH 4/6] tests: Cover symmetric placement and per-equation interp mode Four of the five fail without the preceding commits. The fifth guards against over-reach and passes either way: a derivative already sitting on the target must be left alone. --- tests/test_differentiable.py | 121 ++++++++++++++++++++++++++++++++++- 1 file changed, 120 insertions(+), 1 deletion(-) diff --git a/tests/test_differentiable.py b/tests/test_differentiable.py index 9708f7c59e..444c8456e4 100644 --- a/tests/test_differentiable.py +++ b/tests/test_differentiable.py @@ -6,7 +6,7 @@ from devito import NODE, Differentiable, Eq, Function, Grid, Operator from devito.finite_differences.differentiable import ( - Add, Mul, Pow, SafeInv, diffify, interp_for_fd + Add, EvalDerivative, Mul, Pow, SafeInv, diffify, interp_for_fd ) @@ -534,3 +534,122 @@ def test_direct_breaks_adjoint(self): f' = {inner_e!r}, = {inner_t!r} ' f'(rel diff {rel:.3e})' ) + + +class TestPerEquationInterpMode: + """ + `Eq(..., interp_mode=...)` overrides the Operator's `sym_opt`. + + An Operator generally wants one mode for all of it -- `'direct'`, which + keeps finite-difference stencils compact -- while a single equation in it + needs the other. Accumulating a gradient alongside the propagation it is + the adjoint of is the case this exists for: the stencils have to stay + compact and the accumulation has to be an exact transpose. + """ + + @staticmethod + def _setup(): + grid = Grid((11, 11)) + x, y = grid.dimensions + return (grid, + Function(name='g', grid=grid, space_order=4, staggered=NODE), + Function(name='a', grid=grid, space_order=4, staggered=(x, y)), + Function(name='b', grid=grid, space_order=4, staggered=(x, y))) + + @staticmethod + def _lowered(eq, mode): + return str(Operator([eq], sym_opt={'interp-mode': mode})) + + def test_equation_overrides_operator(self): + """Either mode can be asked for on one equation of an Operator.""" + _, g, a, b = self._setup() + direct = self._lowered(Eq(g, a * b), 'direct') + symmetric = self._lowered(Eq(g, a * b), 'symmetric') + assert direct != symmetric + + # An Operator in 'direct' with the equation asking for 'symmetric' + assert self._lowered(Eq(g, a * b, interp_mode='symmetric'), + 'direct') == symmetric + # ... and the other way round + assert self._lowered(Eq(g, a * b, interp_mode='direct'), + 'symmetric') == direct + + def test_mode_survives_a_wrapper(self): + """ + The mode reaches a product wrapped in one that cannot be re-associated. + + A `dt` scaling or a sum of per-component contractions is exactly such a + wrapper, and dropping the mode there evaluates everything under it in + `'direct'` while reporting that `'symmetric'` was asked for. + """ + grid, g, a, b = self._setup() + dt = grid.stepping_dim.spacing + bare = self._lowered(Eq(g, a * b, interp_mode='symmetric'), 'direct') + scaled = self._lowered(Eq(g, dt * (a * b), interp_mode='symmetric'), + 'direct') + # the scaling is the only difference; the placement must not change + assert scaled.count('a[') == bare.count('a[') + assert scaled.count('b[') == bare.count('b[') + assert scaled != self._lowered(Eq(g, dt * (a * b)), 'direct') + + +class TestSymmetricPlacement: + """ + What `'symmetric'` re-associates, and what it deliberately leaves alone. + """ + + @staticmethod + def _both(eq): + return (str(Operator([eq], sym_opt={'interp-mode': 'direct'})), + str(Operator([eq], sym_opt={'interp-mode': 'symmetric'}))) + + def test_node_coefficient_does_not_block_it(self): + """ + A node-centred coefficient beside staggered operands is the norm for a + gradient -- the accumulator is node-centred and so is the model + parameter it belongs to -- so it must not stop the product being formed + where the operands are. + """ + grid = Grid((11, 11)) + x, y = grid.dimensions + g = Function(name='g', grid=grid, space_order=4, staggered=NODE) + c = Function(name='c', grid=grid, space_order=4, staggered=NODE) + u = Function(name='u', grid=grid, space_order=4, staggered=(x, y)) + v = Function(name='v', grid=grid, space_order=4, staggered=(x, y)) + direct, symmetric = self._both(Eq(g, c * u * v)) + assert direct != symmetric + + def test_a_derivative_on_the_target_is_left_alone(self): + """ + `div(v)` of a staggered velocity lands on the node. + + Re-associating a product around it would replace its compact stencil + with an interpolated one twice as wide -- a different and much more + dispersive discretization, not a re-association of the same one. + """ + grid = Grid((11, 11)) + x, y = grid.dimensions + g = Function(name='g', grid=grid, space_order=4, staggered=NODE) + c = Function(name='c', grid=grid, space_order=4, staggered=NODE) + vx = Function(name='vx', grid=grid, space_order=4, staggered=x) + vy = Function(name='vy', grid=grid, space_order=4, staggered=y) + # vx.dx and vy.dy both land on the node, so does their sum + direct, symmetric = self._both(Eq(g, c * (vx.dx + vy.dy))) + assert direct == symmetric + + def test_degenerate_interpolation_of_a_sum(self): + """ + A zero-order derivative whose weights collapse to one is the identity. + + It rebuilds as whatever it was applied to, and when that is a sum -- + the engineering shear strain of a staggered velocity, say -- the + rebuild used to assert rather than accept the degenerate result. It is + reachable through the symmetric placement of a stiffness contraction. + """ + grid = Grid((11, 11)) + x, y = grid.dimensions + vx = Function(name='vx', grid=grid, space_order=4, staggered=x) + vy = Function(name='vy', grid=grid, space_order=4, staggered=y) + shear = vx.dy + vy.dx + # the identity rebuild: one argument, and it is the sum itself + assert EvalDerivative(shear, base=shear) == shear From a5a789cdeaffcf507a9bfccff9e7841e0a6daae7 Mon Sep 17 00:00:00 2001 From: mloubout Date: Tue, 11 Aug 2026 13:17:46 -0400 Subject: [PATCH 5/6] examples: Show why symmetric interpolation is needed for gradients The staggered-interpolation notebook showed the two modes but not a case that needs the symmetric one. Add the gradient case it exists for: a nodal parameter read at a staggered location, and the matching accumulation back onto the node, where the standard placement is off by 4e-1 against 2e-7. --- .../userapi/08_staggered_interpolation.ipynb | 520 +++++++++++++++--- 1 file changed, 430 insertions(+), 90 deletions(-) diff --git a/examples/userapi/08_staggered_interpolation.ipynb b/examples/userapi/08_staggered_interpolation.ipynb index ec8795f487..64bb995a75 100644 --- a/examples/userapi/08_staggered_interpolation.ipynb +++ b/examples/userapi/08_staggered_interpolation.ipynb @@ -72,10 +72,10 @@ "id": "c2", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:29.927586Z", - "iopub.status.busy": "2026-05-12T14:40:29.927293Z", - "iopub.status.idle": "2026-05-12T14:40:31.299379Z", - "shell.execute_reply": "2026-05-12T14:40:31.298910Z" + "iopub.execute_input": "2026-08-11T20:10:32.563233Z", + "iopub.status.busy": "2026-08-11T20:10:32.562869Z", + "iopub.status.idle": "2026-08-11T20:10:33.750380Z", + "shell.execute_reply": "2026-08-11T20:10:33.749997Z" } }, "outputs": [], @@ -107,10 +107,10 @@ "id": "c4", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.301235Z", - "iopub.status.busy": "2026-05-12T14:40:31.300970Z", - "iopub.status.idle": "2026-05-12T14:40:31.303818Z", - "shell.execute_reply": "2026-05-12T14:40:31.303603Z" + "iopub.execute_input": "2026-08-11T20:10:33.752118Z", + "iopub.status.busy": "2026-08-11T20:10:33.751903Z", + "iopub.status.idle": "2026-08-11T20:10:33.754637Z", + "shell.execute_reply": "2026-08-11T20:10:33.754439Z" } }, "outputs": [ @@ -149,10 +149,10 @@ "id": "c6", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.305053Z", - "iopub.status.busy": "2026-05-12T14:40:31.304943Z", - "iopub.status.idle": "2026-05-12T14:40:31.310855Z", - "shell.execute_reply": "2026-05-12T14:40:31.310542Z" + "iopub.execute_input": "2026-08-11T20:10:33.755843Z", + "iopub.status.busy": "2026-08-11T20:10:33.755746Z", + "iopub.status.idle": "2026-08-11T20:10:33.760860Z", + "shell.execute_reply": "2026-08-11T20:10:33.760594Z" } }, "outputs": [ @@ -194,10 +194,10 @@ "id": "c8", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.312211Z", - "iopub.status.busy": "2026-05-12T14:40:31.312065Z", - "iopub.status.idle": "2026-05-12T14:40:31.405486Z", - "shell.execute_reply": "2026-05-12T14:40:31.405198Z" + "iopub.execute_input": "2026-08-11T20:10:33.762076Z", + "iopub.status.busy": "2026-08-11T20:10:33.761988Z", + "iopub.status.idle": "2026-08-11T20:10:33.859441Z", + "shell.execute_reply": "2026-08-11T20:10:33.859057Z" } }, "outputs": [ @@ -303,10 +303,10 @@ "id": "c12", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.406996Z", - "iopub.status.busy": "2026-05-12T14:40:31.406859Z", - "iopub.status.idle": "2026-05-12T14:40:31.579916Z", - "shell.execute_reply": "2026-05-12T14:40:31.579652Z" + "iopub.execute_input": "2026-08-11T20:10:33.861162Z", + "iopub.status.busy": "2026-08-11T20:10:33.861051Z", + "iopub.status.idle": "2026-08-11T20:10:34.043193Z", + "shell.execute_reply": "2026-08-11T20:10:34.042914Z" } }, "outputs": [ @@ -351,12 +351,12 @@ "$\\mathbf{C}$ a symmetric $6 \\times 6$ stiffness. On the standard staggered\n", "grid for elastodynamics:\n", "\n", - "| Voigt index | Field | Location | `staggered` |\n", - "|-------------|----------------|-----------|-------------|\n", - "| 1, 2, 3 | normal | cell centre | `NODE` |\n", - "| 4 | shear $yz$ | $yz$-edge | `(y, z)` |\n", - "| 5 | shear $xz$ | $xz$-edge | `(x, z)` |\n", - "| 6 | shear $xy$ | $xy$-edge | `(x, y)` |\n", + "| Voigt index | Field | Location | `staggered` |\n", + "|-------------|------------|-------------|-------------|\n", + "| 1, 2, 3 | normal | cell centre | `NODE` |\n", + "| 4 | shear $yz$ | $yz$-edge | `(y, z)` |\n", + "| 5 | shear $xz$ | $xz$-edge | `(x, z)` |\n", + "| 6 | shear $xy$ | $xy$-edge | `(x, y)` |\n", "\n", "and the stiffness coefficients $C_{ij}$ live at the cell centre.\n", "\n", @@ -377,10 +377,10 @@ "id": "c14", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.581219Z", - "iopub.status.busy": "2026-05-12T14:40:31.581128Z", - "iopub.status.idle": "2026-05-12T14:40:31.610846Z", - "shell.execute_reply": "2026-05-12T14:40:31.610567Z" + "iopub.execute_input": "2026-08-11T20:10:34.044598Z", + "iopub.status.busy": "2026-08-11T20:10:34.044509Z", + "iopub.status.idle": "2026-08-11T20:10:34.073772Z", + "shell.execute_reply": "2026-08-11T20:10:34.073459Z" } }, "outputs": [], @@ -419,18 +419,22 @@ "id": "c16", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.612403Z", - "iopub.status.busy": "2026-05-12T14:40:31.612312Z", - "iopub.status.idle": "2026-05-12T14:40:31.614405Z", - "shell.execute_reply": "2026-05-12T14:40:31.614157Z" + "iopub.execute_input": "2026-08-11T20:10:34.075152Z", + "iopub.status.busy": "2026-08-11T20:10:34.075045Z", + "iopub.status.idle": "2026-08-11T20:10:34.077115Z", + "shell.execute_reply": "2026-08-11T20:10:34.076890Z" } }, "outputs": [], "source": [ "def show(eq, mode):\n", " op = Operator(eq, sym_opt={'interp-mode': mode})\n", - " [u] = [n.expr for n in FindNodes(Expression).visit(op)\n", - " if n.expr.lhs.function is eq.lhs.function]\n", + " matches = [\n", + " n.expr\n", + " for n in FindNodes(Expression).visit(op)\n", + " if n.expr.lhs.function is eq.lhs.function\n", + " ]\n", + " [u] = matches\n", " return u\n" ] }, @@ -449,10 +453,10 @@ "id": "c18", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.615775Z", - "iopub.status.busy": "2026-05-12T14:40:31.615692Z", - "iopub.status.idle": "2026-05-12T14:40:31.680333Z", - "shell.execute_reply": "2026-05-12T14:40:31.679899Z" + "iopub.execute_input": "2026-08-11T20:10:34.078263Z", + "iopub.status.busy": "2026-08-11T20:10:34.078184Z", + "iopub.status.idle": "2026-08-11T20:10:34.160262Z", + "shell.execute_reply": "2026-08-11T20:10:34.159791Z" } }, "outputs": [ @@ -488,10 +492,10 @@ "id": "c20", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.682438Z", - "iopub.status.busy": "2026-05-12T14:40:31.682217Z", - "iopub.status.idle": "2026-05-12T14:40:31.783991Z", - "shell.execute_reply": "2026-05-12T14:40:31.783719Z" + "iopub.execute_input": "2026-08-11T20:10:34.162211Z", + "iopub.status.busy": "2026-08-11T20:10:34.162013Z", + "iopub.status.idle": "2026-08-11T20:10:34.266101Z", + "shell.execute_reply": "2026-08-11T20:10:34.265696Z" } }, "outputs": [ @@ -499,7 +503,13 @@ "name": "stdout", "output_type": "stream", "text": [ - "direct : Eq(s1[x + 4, y + 4, z + 4], (0.25*e4[x + 4, y + 3, z + 3] + 0.25*e4[x + 4, y + 4, z + 3] + 0.25*e4[x + 4, y + 3, z + 4] + 0.25*e4[x + 4, y + 4, z + 4])*C14[x + 4, y + 4, z + 4])\n", + "direct : Eq(s1[x + 4, y + 4, z + 4], (0.25*e4[x + 4, y + 3, z + 3] + 0.25*e4[x + 4, y + 4, z + 3] + 0.25*e4[x + 4, y + 3, z + 4] + 0.25*e4[x + 4, y + 4, z + 4])*C14[x + 4, y + 4, z + 4])\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ "symmetric: Eq(s1[x + 4, y + 4, z + 4], (0.25*e4[x + 4, y + 3, z + 3] + 0.25*e4[x + 4, y + 4, z + 3] + 0.25*e4[x + 4, y + 3, z + 4] + 0.25*e4[x + 4, y + 4, z + 4])*C14[x + 4, y + 4, z + 4])\n" ] } @@ -531,10 +541,10 @@ "id": "c22", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.785521Z", - "iopub.status.busy": "2026-05-12T14:40:31.785424Z", - "iopub.status.idle": "2026-05-12T14:40:31.904302Z", - "shell.execute_reply": "2026-05-12T14:40:31.904083Z" + "iopub.execute_input": "2026-08-11T20:10:34.269340Z", + "iopub.status.busy": "2026-08-11T20:10:34.269076Z", + "iopub.status.idle": "2026-08-11T20:10:34.403933Z", + "shell.execute_reply": "2026-08-11T20:10:34.403383Z" } }, "outputs": [ @@ -579,10 +589,10 @@ "id": "c24", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:31.905740Z", - "iopub.status.busy": "2026-05-12T14:40:31.905626Z", - "iopub.status.idle": "2026-05-12T14:40:32.058302Z", - "shell.execute_reply": "2026-05-12T14:40:32.058048Z" + "iopub.execute_input": "2026-08-11T20:10:34.405748Z", + "iopub.status.busy": "2026-08-11T20:10:34.405650Z", + "iopub.status.idle": "2026-08-11T20:10:34.574179Z", + "shell.execute_reply": "2026-08-11T20:10:34.573917Z" } }, "outputs": [ @@ -590,15 +600,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "direct :" + "direct : Eq(s4[x + 4, y + 4, z + 4], (0.25*C45[x + 4, y + 4, z + 4] + 0.25*C45[x + 4, y + 5, z + 4] + 0.25*C45[x + 4, y + 4, z + 5] + 0.25*C45[x + 4, y + 5, z + 5])*(0.25*e5[x + 3, y + 4, z + 4] + 0.25*e5[x + 4, y + 4, z + 4] + 0.25*e5[x + 3, y + 5, z + 4] + 0.25*e5[x + 4, y + 5, z + 4]))\n", + "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - " Eq(s4[x + 4, y + 4, z + 4], (0.25*C45[x + 4, y + 4, z + 4] + 0.25*C45[x + 4, y + 5, z + 4] + 0.25*C45[x + 4, y + 4, z + 5] + 0.25*C45[x + 4, y + 5, z + 5])*(0.25*e5[x + 3, y + 4, z + 4] + 0.25*e5[x + 4, y + 4, z + 4] + 0.25*e5[x + 3, y + 5, z + 4] + 0.25*e5[x + 4, y + 5, z + 4]))\n", - "\n", "symmetric: Eq(s4[x + 4, y + 4, z + 4], 0.5*(r0[z] + r0[z + 1]))\n" ] } @@ -627,10 +636,10 @@ "id": "c26", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:32.059688Z", - "iopub.status.busy": "2026-05-12T14:40:32.059602Z", - "iopub.status.idle": "2026-05-12T14:40:32.164584Z", - "shell.execute_reply": "2026-05-12T14:40:32.164357Z" + "iopub.execute_input": "2026-08-11T20:10:34.575703Z", + "iopub.status.busy": "2026-08-11T20:10:34.575584Z", + "iopub.status.idle": "2026-08-11T20:10:34.687709Z", + "shell.execute_reply": "2026-08-11T20:10:34.687433Z" } }, "outputs": [ @@ -638,7 +647,13 @@ "name": "stdout", "output_type": "stream", "text": [ - "direct : Eq(s4[x + 4, y + 4, z + 4], (0.25*C44[x + 4, y + 4, z + 4] + 0.25*C44[x + 4, y + 5, z + 4] + 0.25*C44[x + 4, y + 4, z + 5] + 0.25*C44[x + 4, y + 5, z + 5])*e4[x + 4, y + 4, z + 4])\n", + "direct : Eq(s4[x + 4, y + 4, z + 4], (0.25*C44[x + 4, y + 4, z + 4] + 0.25*C44[x + 4, y + 5, z + 4] + 0.25*C44[x + 4, y + 4, z + 5] + 0.25*C44[x + 4, y + 5, z + 5])*e4[x + 4, y + 4, z + 4])\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ "symmetric: Eq(s4[x + 4, y + 4, z + 4], (0.25*C44[x + 4, y + 4, z + 4] + 0.25*C44[x + 4, y + 5, z + 4] + 0.25*C44[x + 4, y + 4, z + 5] + 0.25*C44[x + 4, y + 5, z + 5])*e4[x + 4, y + 4, z + 4])\n" ] } @@ -675,10 +690,10 @@ "id": "c28", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:32.165895Z", - "iopub.status.busy": "2026-05-12T14:40:32.165823Z", - "iopub.status.idle": "2026-05-12T14:40:32.207425Z", - "shell.execute_reply": "2026-05-12T14:40:32.207187Z" + "iopub.execute_input": "2026-08-11T20:10:34.689191Z", + "iopub.status.busy": "2026-08-11T20:10:34.689086Z", + "iopub.status.idle": "2026-08-11T20:10:34.739170Z", + "shell.execute_reply": "2026-08-11T20:10:34.738896Z" } }, "outputs": [], @@ -714,11 +729,15 @@ "\n", "def run(mode):\n", " Operator(eqns, sym_opt={'interp-mode': mode}).apply()\n", - " lhs = sum(float(np.dot(e1[i].data.flatten(), e2[i].data.flatten()))\n", - " for i in range(1, 7))\n", - " rhs = sum(float(np.dot(t1[i].data.flatten(), t2[i].data.flatten()))\n", - " for i in range(1, 7))\n", - " return lhs, rhs" + " lhs = sum(\n", + " float(np.dot(e1[i].data.flatten(), e2[i].data.flatten()))\n", + " for i in range(1, 7)\n", + " )\n", + " rhs = sum(\n", + " float(np.dot(t1[i].data.flatten(), t2[i].data.flatten()))\n", + " for i in range(1, 7)\n", + " )\n", + " return lhs, rhs\n" ] }, { @@ -727,10 +746,10 @@ "id": "c29", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:32.208811Z", - "iopub.status.busy": "2026-05-12T14:40:32.208715Z", - "iopub.status.idle": "2026-05-12T14:40:35.866848Z", - "shell.execute_reply": "2026-05-12T14:40:35.866548Z" + "iopub.execute_input": "2026-08-11T20:10:34.740643Z", + "iopub.status.busy": "2026-08-11T20:10:34.740536Z", + "iopub.status.idle": "2026-08-11T20:10:37.633754Z", + "shell.execute_reply": "2026-08-11T20:10:37.633370Z" } }, "outputs": [ @@ -770,10 +789,10 @@ "id": "c30", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:35.868269Z", - "iopub.status.busy": "2026-05-12T14:40:35.868150Z", - "iopub.status.idle": "2026-05-12T14:40:35.870223Z", - "shell.execute_reply": "2026-05-12T14:40:35.870036Z" + "iopub.execute_input": "2026-08-11T20:10:37.635705Z", + "iopub.status.busy": "2026-08-11T20:10:37.635408Z", + "iopub.status.idle": "2026-08-11T20:10:37.638703Z", + "shell.execute_reply": "2026-08-11T20:10:37.638366Z" } }, "outputs": [ @@ -815,10 +834,10 @@ "id": "c32", "metadata": { "execution": { - "iopub.execute_input": "2026-05-12T14:40:35.871434Z", - "iopub.status.busy": "2026-05-12T14:40:35.871343Z", - "iopub.status.idle": "2026-05-12T14:40:35.872864Z", - "shell.execute_reply": "2026-05-12T14:40:35.872647Z" + "iopub.execute_input": "2026-08-11T20:10:37.640420Z", + "iopub.status.busy": "2026-08-11T20:10:37.640292Z", + "iopub.status.idle": "2026-08-11T20:10:37.642099Z", + "shell.execute_reply": "2026-08-11T20:10:37.641855Z" } }, "outputs": [], @@ -831,18 +850,339 @@ "cell_type": "markdown", "id": "c33", "metadata": {}, + "source": [ + "## Why a gradient needs it\n", + "\n", + "The stiffness above is one $\\mathbf{I}\\,\\mathbf{A}\\,\\mathbf{I}^{\\!\\top}$\n", + "operator among others. The case that forces the issue in practice is a\n", + "*gradient*: an adjoint-state method accumulates $\\partial F/\\partial m$ into a\n", + "field that lives where the model parameter lives -- the node -- out of\n", + "wavefields that live on faces and corners.\n", + "\n", + "That accumulation is a transpose, and it is the transpose of something the\n", + "forward already does. Where a node-centred parameter multiplies a staggered\n", + "quantity, the forward reads it through an interpolation,\n", + "\n", + "$$ v_i \\mathrel{+}= \\Delta t \\; (\\mathbf{I}\\,b)_i \\; (\\nabla\\!\\cdot\\tau)_i , $$\n", + "\n", + "so the derivative with respect to the *nodal* $b$ spreads the staggered\n", + "residual back over the nodes that fed it, $\\mathbf{I}^{\\!\\top}$. Forming the\n", + "product where its operands are and interpolating once is exactly that.\n", + "Interpolating the operands one at a time instead computes\n", + "$\\langle a\\rangle\\langle b\\rangle$ where the transpose needs\n", + "$\\langle ab\\rangle$ -- second-order accurate, and not a transpose. The\n", + "symptom is not a wrong-looking gradient; it is a gradient whose Taylor test\n", + "converges at first order instead of second.\n", + "\n", + "Here is that pattern on its own: a node-centred `m` read at a staggered\n", + "location on the way in, and the matching accumulation on the way back." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "c34", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T20:10:37.643370Z", + "iopub.status.busy": "2026-08-11T20:10:37.643243Z", + "iopub.status.idle": "2026-08-11T20:10:38.136764Z", + "shell.execute_reply": "2026-08-11T20:10:38.136481Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x-face direct relative adjoint error 4.06e-01\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x-face symmetric relative adjoint error 2.07e-07\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "corner direct relative adjoint error 8.79e-01\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "corner symmetric relative adjoint error 7.36e-07\n" + ] + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "import numpy as np\n", + "from devito import Grid, Function, Eq, Operator, NODE\n", + "\n", + "np.random.seed(0)\n", + "grid = Grid(shape=(24, 24), extent=(1., 1.))\n", + "x, y = grid.dimensions\n", + "\n", + "\n", + "def adjoint_test(loc, mode):\n", + " m = Function(name='m', grid=grid, space_order=4, staggered=NODE)\n", + " g = Function(name='g', grid=grid, space_order=4, staggered=NODE)\n", + " S = Function(name='S', grid=grid, space_order=4, staggered=loc)\n", + " t = Function(name='t', grid=grid, space_order=4, staggered=loc)\n", + " q = Function(name='q', grid=grid, space_order=4, staggered=loc)\n", + "\n", + " # support strictly inside, so no edge term enters the identity\n", + " inner = (slice(6, -6), slice(6, -6))\n", + " for f in (m, S, q):\n", + " f.data[:] = 0.\n", + " f.data[inner] = np.random.rand(*f.data[inner].shape) - .5\n", + "\n", + " # forward: the propagation, always the compact `'direct'` stencil\n", + " Operator([Eq(t, m * S)], sym_opt={'interp-mode': 'direct'}).apply()\n", + " # adjoint: the accumulation back onto the node, in the mode under test\n", + " Operator([Eq(g, q * S)], sym_opt={'interp-mode': mode}).apply()\n", + "\n", + " lhs = float(np.sum(np.array(t.data) * np.array(q.data))) # \n", + " rhs = float(np.sum(np.array(m.data) * np.array(g.data))) # \n", + " return abs(lhs - rhs) / max(abs(lhs), abs(rhs))\n", + "\n", + "\n", + "for name, loc in (('x-face', (x,)), ('corner', (x, y))):\n", + " for mode in ('direct', 'symmetric'):\n", + " print(\n", + " f\"{name:7s} {mode:9s} relative adjoint error \"\n", + " f\"{adjoint_test(loc, mode):.2e}\"\n", + " )\n" + ] + }, + { + "cell_type": "markdown", + "id": "c35", + "metadata": {}, + "source": [ + "Only the accumulation changes mode -- the forward keeps the compact\n", + "`'direct'` stencil throughout, as a propagation should. With that,\n", + "`'symmetric'` satisfies $\\langle \\mathbf{I}m\\,S, q\\rangle =\n", + "\\langle m, \\mathbf{I}^{\\!\\top}(qS)\\rangle$ to float precision on faces and on\n", + "corners alike, while `'direct'` is out by tens of per cent. That discrepancy\n", + "is what turns a second-order gradient into a first-order one." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "c36", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T20:10:38.138182Z", + "iopub.status.busy": "2026-08-11T20:10:38.138072Z", + "iopub.status.idle": "2026-08-11T20:10:38.659080Z", + "shell.execute_reply": "2026-08-11T20:10:38.658846Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + } + ], + "source": [ + "for loc in ((x,), (x, y)):\n", + " assert adjoint_test(loc, 'symmetric') < 1e-5\n", + " assert adjoint_test(loc, 'direct') > 1e-2\n" + ] + }, + { + "cell_type": "markdown", + "id": "c37", + "metadata": {}, + "source": [ + "### Asking for it on one equation\n", + "\n", + "An Operator usually wants `'direct'` throughout: it is the compact stencil,\n", + "and re-associating a product around a derivative that already sits on the\n", + "target would replace that derivative with an interpolated one twice as wide.\n", + "A gradient accumulated *inside* the time loop of the propagation it is the\n", + "adjoint of therefore needs the two modes at once -- compact stencils, exact\n", + "transpose on the accumulation. `Eq` takes the override:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "c38", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-11T20:10:38.660543Z", + "iopub.status.busy": "2026-08-11T20:10:38.660438Z", + "iopub.status.idle": "2026-08-11T20:10:38.767322Z", + "shell.execute_reply": "2026-08-11T20:10:38.767079Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "symmetric placement on a direct Operator: True\n" + ] + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "u = Function(name='u', grid=grid, space_order=4, staggered=(x, y))\n", + "w = Function(name='w', grid=grid, space_order=4, staggered=(x, y))\n", + "acc = Function(name='acc', grid=grid, space_order=4, staggered=NODE)\n", + "\n", + "# the Operator is 'direct'; this one equation is not\n", + "op = Operator(\n", + " [Eq(acc, acc + u * w, interp_mode='symmetric')],\n", + " sym_opt={'interp-mode': 'direct'},\n", + ")\n", + "print(\n", + " 'symmetric placement on a direct Operator:',\n", + " str(op) == str(\n", + " Operator(\n", + " [Eq(acc, acc + u * w)],\n", + " sym_opt={'interp-mode': 'symmetric'},\n", + " )\n", + " ),\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "id": "c39", + "metadata": {}, "source": [ "## When to use which\n", "\n", - "| Situation | Mode |\n", - "|------------------------------------------------------------------------|---------------|\n", - "| Acoustic / scalar-wave equations | `'direct'` |\n", - "| Elastic stress-strain or any $\\mathbf{I}\\,\\mathbf{A}\\,\\mathbf{I}^{\\!\\top}$ operator | `'symmetric'` |\n", - "| Adjoint-state inversion needing exact discrete adjoint | `'symmetric'` |\n", - "| Any equation where one factor already matches the target staggering | either |\n", + "| Situation | Mode |\n", + "|---------------------------------------------------------------------|---------------|\n", + "| Acoustic / scalar-wave equations | `'direct'` |\n", + "| Elastic stress-strain or any $\\mathbf{IAI}^{\\!\\top}$ operator | `'symmetric'` |\n", + "| Adjoint-state inversion needing exact discrete adjoint | `'symmetric'` |\n", + "| A gradient accumulated from staggered fields onto a nodal parameter | `'symmetric'` |\n", + "| Any equation where one factor already matches the target staggering | either |\n", "\n", "`'direct'` is the default because it is the cheaper and smaller stencil; pick\n", - "`'symmetric'` deliberately when you need the adjoint structure preserved." + "`'symmetric'` deliberately when you need the adjoint structure preserved --\n", + "per Operator with `sym_opt`, or per equation with `Eq(..., interp_mode=...)`\n", + "when only the accumulation needs it." ] } ], @@ -862,7 +1202,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.13" + "version": "3.13.14" } }, "nbformat": 4, From ec12fb13b7feb0d0d368deb072a12bef86d93344 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 12 Aug 2026 08:50:17 -0400 Subject: [PATCH 6/6] dsl: mute sympy matrix warning --- devito/types/basic.py | 48 +++++++++++++++++++++++++++++++------------ tests/test_tensors.py | 38 ++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 13 deletions(-) diff --git a/devito/types/basic.py b/devito/types/basic.py index 5456c6220b..81d20bb670 100644 --- a/devito/types/basic.py +++ b/devito/types/basic.py @@ -1,7 +1,7 @@ import abc import inspect import warnings -from contextlib import suppress +from contextlib import contextmanager, suppress from ctypes import POINTER, Structure, _Pointer, c_char, c_char_p from functools import cached_property, reduce from operator import mul @@ -1471,6 +1471,22 @@ def __getnewargs_ex__(self): return args, kwargs +@contextmanager +def ignore_non_expr_deprecation(): + """ + Suppress sympy's deprecation of non-`Expr` entries in a Matrix. + + A Devito tensor is a Matrix of Devito objects, some of which are legitimately + not `Expr` (for example a serialization string), so the deprecation does not + apply to us. Sympy emits it from `_dod_to_DomainMatrix`, reached whenever a + tensor is built from components, and from `_unify_element_sympy`, reached on + scalar multiplication, so both are wrapped in this. + """ + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=SymPyDeprecationWarning) + yield + + class AbstractTensor(sympy.ImmutableDenseMatrix, Basic, Pickable, Evaluable): """ @@ -1563,23 +1579,29 @@ def __subfunc_setup__(cls, *args, **kwargs): return [] @classmethod - def _sympify(self, arg): + 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 + # 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) try: - # Pure sympy object - return arg._sympy_() - except AttributeError: + return sympy.sympify(arg, strict=True) + except sympy.SympifyError: return arg @classmethod - def _eval_from_dok(cls, rows, cols, dok): - with warnings.catch_warnings(): - warnings.filterwarnings( - "ignore", - category=SymPyDeprecationWarning - ) - return super()._eval_from_dok(rows, cols, dok) + def _dod_to_DomainMatrix(cls, rows, cols, dod, types): + # Entry point of every matrix construction, from either a flat list + # (`_new`) or a dok (`_eval_from_dok`) + with ignore_non_expr_deprecation(): + return super()._dod_to_DomainMatrix(rows, cols, dod, types) + + @classmethod + def _unify_element_sympy(cls, rep, element): + # Entry point of scalar multiplication, e.g. `lam * tau` + with ignore_non_expr_deprecation(): + return super()._unify_element_sympy(rep, element) @property def grid(self): diff --git a/tests/test_tensors.py b/tests/test_tensors.py index 512c48cf55..53725a493e 100644 --- a/tests/test_tensors.py +++ b/tests/test_tensors.py @@ -1,7 +1,10 @@ +import warnings + import numpy as np import pytest import sympy from sympy import Matrix, Rational +from sympy.utilities.exceptions import SymPyDeprecationWarning from devito import ( Dimension, Eq, Function, Grid, TensorFunction, TensorTimeFunction, TimeFunction, @@ -510,6 +513,41 @@ def test_diag(func1): assert all(f2[i, i] == f1 for i in range(3)) +@pytest.mark.parametrize('func1', [Function, TimeFunction, VectorFunction, + VectorTimeFunction]) +def test_diag_sympified_zeros(func1): + """The zeros filling a diagonal tensor must be `Expr`, not plain `int`.""" + grid = Grid(tuple([5]*3)) + f1 = func1(name="f1", grid=grid) + + with warnings.catch_warnings(): + warnings.simplefilter("error", SymPyDeprecationWarning) + f2 = diag(div(f1) if isinstance(f1, (VectorFunction, + VectorTimeFunction)) else f1) + + assert all(isinstance(c, sympy.Expr) for c in f2.flat()) + + +def test_non_expr_components(): + """ + A tensor may legitimately hold non-`Expr` components, which sympy deprecates + for plain Matrices. None of the construction or arithmetic paths must warn. + """ + grid = Grid(tuple([5]*3)) + f1 = Function(name="f1", grid=grid) + # `S.true` is `Basic` but not `Expr`, as e.g. a serialization string would be + comps = [[f1 if i == j else sympy.S.true for i in range(3)] for j in range(3)] + + with warnings.catch_warnings(): + warnings.simplefilter("error", SymPyDeprecationWarning) + # Construction, from a flat list of components... + f2 = TensorFunction(name="f2", grid=grid, components=comps) + # ...and from a dok, on rebuild + assert isinstance(f2.T, TensorFunction) + # Scalar multiplication, which unifies the scalar with the components + assert isinstance(f1 * TensorFunction(name="f3", grid=grid), TensorFunction) + + @pytest.mark.parametrize('func1', [TensorFunction, VectorFunction]) def test_kwargs(func1): orders = Matrix([[1, 2], [3, 4]]) if func1 is TensorFunction else Matrix([1, 2])