diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py
index b77d2c158d..88e3cb214e 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)
+        # 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).indices_ref
+        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)
@@ -1119,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
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/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
 
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)))   # <m*S, q>\n",
+    "    rhs = float(np.sum(np.array(m.data) * np.array(g.data)))   # <m, (q*S)^T>\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,
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'<e1, C*t2> = {inner_e!r}, <C*e1, t2> = {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
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])