From 9f93edd2508d9265093ba318224d3e4d5ad063a6 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Fri, 19 Dec 2025 17:18:44 +0000
Subject: [PATCH 1/4] api: Extend side API to all derivative operator shortcuts

---
 devito/finite_differences/derivative.py     |  2 +-
 devito/finite_differences/differentiable.py | 29 ++++++---
 devito/finite_differences/operators.py      | 24 ++++++--
 devito/types/tensor.py                      | 57 ++++++++++++++----
 tests/test_derivatives.py                   | 65 ++++++++++++++++++++-
 5 files changed, 152 insertions(+), 25 deletions(-)

diff --git a/devito/finite_differences/derivative.py b/devito/finite_differences/derivative.py
index bcfff8b550..e3a03e3718 100644
--- a/devito/finite_differences/derivative.py
+++ b/devito/finite_differences/derivative.py
@@ -38,7 +38,7 @@ class Derivative(sympy.Derivative, Differentiable, Pickable):
         Derivative order.
     side : Side or tuple of Side, optional, default=centered
         Side of the finite difference location, centered (at x), left (at x - 1)
-        or right (at x +1).
+        or right (at x + 1).
     transpose : Transpose, optional, default=direct
         Forward (matvec=direct) or transpose (matvec=transpose) mode of the
         finite difference.
diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py
index c45b543398..c2b139badc 100644
--- a/devito/finite_differences/differentiable.py
+++ b/devito/finite_differences/differentiable.py
@@ -303,12 +303,12 @@ def shift(self, dim, shift):
         return self._subs(dim, dim + shift)
 
     @property
-    def laplace(self):
+    def laplace(self, **kwargs):
         """
         Generates a symbolic expression for the Laplacian, the second
         derivative w.r.t all spatial Dimensions.
         """
-        return self.laplacian()
+        return self.laplacian(**kwargs)
 
     def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         """
@@ -329,16 +329,20 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
-            Custom weights for the finite differences.
+            Custom weights for the finite difference coefficients.
         """
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         order = order or self.space_order
         space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
         derivs = tuple(f'd{d.name}2' for d in space_dims)
         return Add(*[getattr(self, d)(x0=shift_x0(shift, space_dims[i], None, i),
-                                      method=method, fd_order=order, w=w)
+                                      method=method, fd_order=order, side=side, w=w)
                      for i, d in enumerate(derivs)])
 
     def div(self, shift=None, order=None, method='FD', **kwargs):
@@ -357,15 +361,20 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
         order = order or self.space_order
         return Add(*[getattr(self, f'd{d.name}')(x0=shift_x0(shift, d, None, i),
-                                                 fd_order=order, method=method, w=w)
+                                                 fd_order=order, method=method, side=side,
+                                                 w=w)
                      for i, d in enumerate(space_dims)])
 
     def grad(self, shift=None, order=None, method='FD', **kwargs):
@@ -384,16 +393,22 @@ def grad(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
-            Custom weights for the finite
+            Custom weights for the finite difference coefficients.
         """
         from devito.types.tensor import VectorFunction, VectorTimeFunction
         space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
         order = order or self.space_order
+
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         comps = [getattr(self, f'd{d.name}')(x0=shift_x0(shift, d, None, i),
-                                             fd_order=order, method=method, w=w)
+                                             fd_order=order, method=method, side=side,
+                                             w=w)
                  for i, d in enumerate(space_dims)]
         vec_func = VectorTimeFunction if self.is_TimeDependent else VectorFunction
         return vec_func(name=f'grad_{self.name}', time_order=self.time_order,
diff --git a/devito/finite_differences/operators.py b/devito/finite_differences/operators.py
index 17452ae664..330454a7fd 100644
--- a/devito/finite_differences/operators.py
+++ b/devito/finite_differences/operators.py
@@ -14,12 +14,16 @@ def div(func, shift=None, order=None, method='FD', **kwargs):
     method: str, optional, default='FD'
         Discretization method. Options are 'FD' (default) and
         'RSFD' (rotated staggered grid finite-difference).
+    side : Side or tuple of Side, optional, default=centered
+        Side of the finite difference location, centered (at x), left (at x - 1)
+        or right (at x + 1).
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
+    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
-        return func.div(shift=shift, order=order, method=method, w=w)
+        return func.div(shift=shift, order=order, method=method, side=side, w=w)
     except AttributeError:
         return 0
 
@@ -57,12 +61,16 @@ def grad(func, shift=None, order=None, method='FD', **kwargs):
     method: str, optional, default='FD'
         Discretization method. Options are 'FD' (default) and
         'RSFD' (rotated staggered grid finite-difference).
+    side : Side or tuple of Side, optional, default=centered
+        Side of the finite difference location, centered (at x), left (at x - 1)
+        or right (at x + 1).
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
+    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
-        return func.grad(shift=shift, order=order, method=method, w=w)
+        return func.grad(shift=shift, order=order, method=method, side=side, w=w)
     except AttributeError:
         raise AttributeError("Gradient not supported for class %s" % func.__class__)
 
@@ -100,12 +108,16 @@ def curl(func, shift=None, order=None, method='FD', **kwargs):
     method: str, optional, default='FD'
         Discretization method. Options are 'FD' (default) and
         'RSFD' (rotated staggered grid finite-difference).
+    side : Side or tuple of Side, optional, default=centered
+        Side of the finite difference location, centered (at x), left (at x - 1)
+        or right (at x + 1).
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
+    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
-        return func.curl(shift=shift, order=order, method=method, w=w)
+        return func.curl(shift=shift, order=order, method=method, side=side, w=w)
     except AttributeError:
         raise AttributeError("Curl only supported for 3D VectorFunction")
 
@@ -143,12 +155,16 @@ def laplace(func, shift=None, order=None, method='FD', **kwargs):
         Uses `func.space_order` when not specified
     method: str, optional, default='FD'
         Discretization method. Options are 'FD' (default) and 'RSFD'
+    side : Side or tuple of Side, optional, default=centered
+        Side of the finite difference location, centered (at x), left (at x - 1)
+        or right (at x + 1).
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
+    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
-        return func.laplacian(shift=shift, order=order, method=method, w=w)
+        return func.laplacian(shift=shift, order=order, method=method, side=side, w=w)
     except AttributeError:
         return 0
 
diff --git a/devito/types/tensor.py b/devito/types/tensor.py
index 30d5e9980b..be2ea60cd7 100644
--- a/devito/types/tensor.py
+++ b/devito/types/tensor.py
@@ -234,9 +234,13 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite differences.
         """
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         comps = []
         func = vec_func(self)
@@ -247,7 +251,7 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         for i in range(len(self.space_dimensions)):
             comps.append(sum([getattr(self[j, i], 'd%s' % d.name)
                               (x0=shift_x0(shift, d, i, j), fd_order=order,
-                               method=method, w=w)
+                               method=method, side=side, w=w)
                               for j, d in enumerate(space_dims)]))
         return func._new(comps)
 
@@ -277,9 +281,13 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite
         """
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         comps = []
         func = vec_func(self)
@@ -290,7 +298,7 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         for j in range(ndim):
             comps.append(sum([getattr(self[j, i], 'd%s2' % d.name)
                               (x0=shift_x0(shift, d, j, i), fd_order=order,
-                               method=method, w=w)
+                               method=method, side=side, w=w)
                               for i, d in enumerate(space_dims)]))
         return func._new(comps)
 
@@ -372,15 +380,20 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         shift_x0 = make_shift_x0(shift, (len(self.space_dimensions),))
         order = order or self.space_order
         space_dims = self.root_dimensions
         return sum([getattr(self[i], 'd%s' % d.name)(x0=shift_x0(shift, d, None, i),
-                                                     fd_order=order, method=method, w=w)
+                                                     fd_order=order, method=method,
+                                                     side=side, w=w)
                     for i, d in enumerate(space_dims)])
 
     @property
@@ -404,16 +417,21 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite
         """
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         func = vec_func(self)
         shift_x0 = make_shift_x0(shift, (len(self.space_dimensions),))
         order = order or self.space_order
         space_dims = self.root_dimensions
         comps = [sum([getattr(s, 'd%s2' % d.name)(x0=shift_x0(shift, d, None, i),
-                                                  fd_order=order, w=w, method=method)
+                                                  fd_order=order, method=method,
+                                                  side=side, w=w)
                       for i, d in enumerate(space_dims)])
                  for s in self]
         return func._new(comps)
@@ -432,29 +450,39 @@ def curl(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
         if len(self.space_dimensions) != 3:
             raise AttributeError("Curl only supported for 3D VectorFunction")
         # The curl of a VectorFunction is a VectorFunction
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         dims = self.root_dimensions
         derivs = ['d%s' % d.name for d in dims]
         shift_x0 = make_shift_x0(shift, (len(dims), len(dims)))
         order = order or self.space_order
         comp1 = (getattr(self[2], derivs[1])(x0=shift_x0(shift, dims[1], 2, 1),
-                                             fd_order=order, method=method, w=w) -
+                                             fd_order=order, method=method,
+                                             side=side, w=w) -
                  getattr(self[1], derivs[2])(x0=shift_x0(shift, dims[2], 1, 2),
-                                             fd_order=order, method=method, w=w))
+                                             fd_order=order, method=method,
+                                             side=side, w=w))
         comp2 = (getattr(self[0], derivs[2])(x0=shift_x0(shift, dims[2], 0, 2),
-                                             fd_order=order, method=method, w=w) -
+                                             fd_order=order, method=method,
+                                             side=side, w=w) -
                  getattr(self[2], derivs[0])(x0=shift_x0(shift, dims[0], 2, 0),
-                                             fd_order=order, method=method, w=w))
+                                             fd_order=order, method=method,
+                                             side=side, w=w))
         comp3 = (getattr(self[1], derivs[0])(x0=shift_x0(shift, dims[0], 1, 0),
-                                             fd_order=order, method=method, w=w) -
+                                             fd_order=order, method=method,
+                                             side=side, w=w) -
                  getattr(self[0], derivs[1])(x0=shift_x0(shift, dims[1], 0, 1),
-                                             fd_order=order, method=method, w=w))
+                                             fd_order=order, method=method,
+                                             side=side, w=w))
         func = vec_func(self)
         return func._new(3, 1, [comp1, comp2, comp3])
 
@@ -472,17 +500,22 @@ def grad(self, shift=None, order=None, method='FD', **kwargs):
         method: str, optional, default='FD'
             Discretization method. Options are 'FD' (default) and
             'RSFD' (rotated staggered grid finite-difference).
+        side : Side or tuple of Side, optional, default=centered
+            Side of the finite difference location, centered (at x), left (at x - 1)
+            or right (at x + 1).
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
+        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         func = tens_func(self)
         ndim = len(self.space_dimensions)
         shift_x0 = make_shift_x0(shift, (ndim, ndim))
         order = order or self.space_order
         space_dims = self.root_dimensions
-        comps = [[getattr(f, 'd%s' % d.name)(x0=shift_x0(shift, d, i, j), w=w,
-                                             fd_order=order, method=method)
+        comps = [[getattr(f, 'd%s' % d.name)(x0=shift_x0(shift, d, i, j),
+                                             fd_order=order, method=method,
+                                             side=side, w=w)
                   for j, d in enumerate(space_dims)]
                  for i, f in enumerate(self)]
         return func._new(comps)
diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py
index 8872abb743..a9b0940d82 100644
--- a/tests/test_derivatives.py
+++ b/tests/test_derivatives.py
@@ -3,7 +3,8 @@
 from sympy import sympify, simplify, diff, Float, Symbol
 
 from devito import (Grid, Function, TimeFunction, Eq, Operator, NODE, cos, sin,
-                    ConditionalDimension, left, right, centered, div, grad)
+                    ConditionalDimension, left, right, centered, div, grad,
+                    curl, laplace, VectorFunction)
 from devito.finite_differences import Derivative, Differentiable, diffify
 from devito.finite_differences.differentiable import (Add, EvalDerivative, IndexSum,
                                                       IndexDerivative, Weights,
@@ -624,6 +625,68 @@ def test_shifted_grad(self, shift, ndim):
                 gk = getattr(f, 'd%s' % d.name)(x0=x0, fd_order=order).evaluate
                 assert gi == gk
 
+    # TODO: Repeat all these with vector and tensor functions where defined
+    @pytest.mark.parametrize('side', [left, right, centered])
+    def test_grad_w_side(self, side):
+        grid = Grid(shape=(11, 11))
+        f = Function(name='f', grid=grid, space_order=2)
+
+        # Want to check that it's the same as constructing by hand
+        comps = (f.dx(side=side), f.dy(side=side))
+        expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
+                               components=comps, grid=grid).evaluate
+
+        expr2 = f.grad(side=side).evaluate
+        expr3 = grad(f, side=side).evaluate
+
+        assert expr1 == expr2
+        assert expr1 == expr3
+
+    @pytest.mark.parametrize('side', [left, right, centered])
+    def test_div_w_side(self, side):
+        grid = Grid(shape=(11, 11))
+        f = VectorFunction(name='f', grid=grid, space_order=2, staggered=(None, None))
+
+        expr1 = (f[0].dx(side=side) + f[1].dy(side=side)).evaluate
+
+        expr2 = f.div(side=side).evaluate
+        expr3 = div(f, side=side).evaluate
+
+        assert expr1 == expr2
+        assert expr1 == expr3
+
+    @pytest.mark.parametrize('side', [left, right, centered])
+    def test_curl_w_side(self, side):
+        grid = Grid(shape=(11, 11, 11))
+        f = VectorFunction(name='f', grid=grid, space_order=2,
+                           staggered=(None, None, None))
+
+        comps = (f[2].dy(side=side) - f[1].dz(side=side),
+                 f[0].dz(side=side) - f[2].dx(side=side),
+                 f[1].dx(side=side) - f[0].dy(side=side))
+
+        expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
+                               components=comps, grid=grid).evaluate
+
+        expr2 = f.curl(side=side).evaluate
+        expr3 = curl(f, side=side).evaluate
+
+        assert expr1 == expr2
+        assert expr1 == expr3
+
+    @pytest.mark.parametrize('side', [left, right, centered])
+    def test_laplace_w_side(self, side):
+        grid = Grid(shape=(11, 11))
+        f = Function(name='f', grid=grid, space_order=2)
+
+        expr1 = (f.dx2(side=side) + f.dy2(side=side)).evaluate
+
+        expr2 = f.laplacian(side=side).evaluate
+        expr3 = laplace(f, side=side).evaluate
+
+        assert expr1 == expr2
+        assert expr1 == expr3
+
     def test_substitution(self):
         grid = Grid((11, 11))
         f = Function(name="f", grid=grid, space_order=4)

From 039ddd0ebcc6b452980d4cb48a318514d5556aaf Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Fri, 19 Dec 2025 17:58:45 +0000
Subject: [PATCH 2/4] tests: Add a vector grad test and simplify existing tests

---
 tests/test_derivatives.py | 42 ++++++++++++++++++++-------------------
 1 file changed, 22 insertions(+), 20 deletions(-)

diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py
index a9b0940d82..c76d4eb22d 100644
--- a/tests/test_derivatives.py
+++ b/tests/test_derivatives.py
@@ -4,7 +4,7 @@
 
 from devito import (Grid, Function, TimeFunction, Eq, Operator, NODE, cos, sin,
                     ConditionalDimension, left, right, centered, div, grad,
-                    curl, laplace, VectorFunction)
+                    curl, laplace, VectorFunction, TensorFunction)
 from devito.finite_differences import Derivative, Differentiable, diffify
 from devito.finite_differences.differentiable import (Add, EvalDerivative, IndexSum,
                                                       IndexDerivative, Weights,
@@ -636,11 +636,22 @@ def test_grad_w_side(self, side):
         expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
                                components=comps, grid=grid).evaluate
 
-        expr2 = f.grad(side=side).evaluate
-        expr3 = grad(f, side=side).evaluate
+        assert expr1 == f.grad(side=side).evaluate
+        assert expr1 == grad(f, side=side).evaluate
 
-        assert expr1 == expr2
-        assert expr1 == expr3
+    @pytest.mark.parametrize('side', [left, right, centered])
+    def test_vector_grad_w_side(self, side):
+        grid = Grid(shape=(11, 11))
+        f = VectorFunction(name='f', grid=grid, space_order=2, staggered=(None, None))
+
+        comps = ((f[0].dx(side=side), f[0].dy(side=side)),
+                 (f[1].dx(side=side), f[1].dy(side=side)))
+
+        expr1 = TensorFunction(name=f"{f.name}_tens", space_order=f.space_order,
+                               components=comps, grid=grid).evaluate
+
+        assert expr1 == f.grad(side=side).evaluate
+        assert expr1 == grad(f, side=side).evaluate
 
     @pytest.mark.parametrize('side', [left, right, centered])
     def test_div_w_side(self, side):
@@ -649,11 +660,8 @@ def test_div_w_side(self, side):
 
         expr1 = (f[0].dx(side=side) + f[1].dy(side=side)).evaluate
 
-        expr2 = f.div(side=side).evaluate
-        expr3 = div(f, side=side).evaluate
-
-        assert expr1 == expr2
-        assert expr1 == expr3
+        assert expr1 == f.div(side=side).evaluate
+        assert expr1 == div(f, side=side).evaluate
 
     @pytest.mark.parametrize('side', [left, right, centered])
     def test_curl_w_side(self, side):
@@ -668,11 +676,8 @@ def test_curl_w_side(self, side):
         expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
                                components=comps, grid=grid).evaluate
 
-        expr2 = f.curl(side=side).evaluate
-        expr3 = curl(f, side=side).evaluate
-
-        assert expr1 == expr2
-        assert expr1 == expr3
+        assert expr1 == f.curl(side=side).evaluate
+        assert expr1 == curl(f, side=side).evaluate
 
     @pytest.mark.parametrize('side', [left, right, centered])
     def test_laplace_w_side(self, side):
@@ -681,11 +686,8 @@ def test_laplace_w_side(self, side):
 
         expr1 = (f.dx2(side=side) + f.dy2(side=side)).evaluate
 
-        expr2 = f.laplacian(side=side).evaluate
-        expr3 = laplace(f, side=side).evaluate
-
-        assert expr1 == expr2
-        assert expr1 == expr3
+        assert expr1 == f.laplacian(side=side).evaluate
+        assert expr1 == laplace(f, side=side).evaluate
 
     def test_substitution(self):
         grid = Grid((11, 11))

From 5454f56d305d334f70d10097860ec2677c4e76e5 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Mon, 22 Dec 2025 10:50:11 +0000
Subject: [PATCH 3/4] tests: Enhanced tests for side

---
 devito/finite_differences/differentiable.py | 13 +++---
 devito/types/tensor.py                      | 18 +++-----
 tests/test_derivatives.py                   | 48 ++++++++++++++++++++-
 3 files changed, 57 insertions(+), 22 deletions(-)

diff --git a/devito/finite_differences/differentiable.py b/devito/finite_differences/differentiable.py
index c2b139badc..cae122dfed 100644
--- a/devito/finite_differences/differentiable.py
+++ b/devito/finite_differences/differentiable.py
@@ -303,14 +303,14 @@ def shift(self, dim, shift):
         return self._subs(dim, dim + shift)
 
     @property
-    def laplace(self, **kwargs):
+    def laplace(self):
         """
         Generates a symbolic expression for the Laplacian, the second
         derivative w.r.t all spatial Dimensions.
         """
-        return self.laplacian(**kwargs)
+        return self.laplacian()
 
-    def laplacian(self, shift=None, order=None, method='FD', **kwargs):
+    def laplacian(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Laplacian of the Differentiable with shifted derivatives and custom
         FD order.
@@ -335,7 +335,6 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         order = order or self.space_order
         space_dims = self.root_dimensions
@@ -345,7 +344,7 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
                                       method=method, fd_order=order, side=side, w=w)
                      for i, d in enumerate(derivs)])
 
-    def div(self, shift=None, order=None, method='FD', **kwargs):
+    def div(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Divergence of the input Function.
 
@@ -367,7 +366,6 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         space_dims = self.root_dimensions
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
@@ -377,7 +375,7 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
                                                  w=w)
                      for i, d in enumerate(space_dims)])
 
-    def grad(self, shift=None, order=None, method='FD', **kwargs):
+    def grad(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Gradient of the input Function.
 
@@ -404,7 +402,6 @@ def grad(self, shift=None, order=None, method='FD', **kwargs):
         shift_x0 = make_shift_x0(shift, (len(space_dims),))
         order = order or self.space_order
 
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         comps = [getattr(self, f'd{d.name}')(x0=shift_x0(shift, d, None, i),
                                              fd_order=order, method=method, side=side,
diff --git a/devito/types/tensor.py b/devito/types/tensor.py
index be2ea60cd7..e83679d11d 100644
--- a/devito/types/tensor.py
+++ b/devito/types/tensor.py
@@ -220,7 +220,7 @@ def values(self):
         else:
             return super().values()
 
-    def div(self, shift=None, order=None, method='FD', **kwargs):
+    def div(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Divergence of the TensorFunction (is a VectorFunction).
 
@@ -240,7 +240,6 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite differences.
         """
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         comps = []
         func = vec_func(self)
@@ -262,7 +261,7 @@ def laplace(self):
         """
         return self.laplacian()
 
-    def laplacian(self, shift=None, order=None, method='FD', **kwargs):
+    def laplacian(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Laplacian of the TensorFunction with shifted derivatives and custom
         FD order.
@@ -287,7 +286,6 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite
         """
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         comps = []
         func = vec_func(self)
@@ -366,7 +364,7 @@ def __str__(self):
 
     __repr__ = __str__
 
-    def div(self, shift=None, order=None, method='FD', **kwargs):
+    def div(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Divergence of the VectorFunction, creates the divergence Function.
 
@@ -386,7 +384,6 @@ def div(self, shift=None, order=None, method='FD', **kwargs):
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         shift_x0 = make_shift_x0(shift, (len(self.space_dimensions),))
         order = order or self.space_order
@@ -403,7 +400,7 @@ def laplace(self):
         """
         return self.laplacian()
 
-    def laplacian(self, shift=None, order=None, method='FD', **kwargs):
+    def laplacian(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Laplacian of the VectorFunction, creates the Laplacian VectorFunction.
 
@@ -423,7 +420,6 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite
         """
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         func = vec_func(self)
         shift_x0 = make_shift_x0(shift, (len(self.space_dimensions),))
@@ -436,7 +432,7 @@ def laplacian(self, shift=None, order=None, method='FD', **kwargs):
                  for s in self]
         return func._new(comps)
 
-    def curl(self, shift=None, order=None, method='FD', **kwargs):
+    def curl(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Gradient of the (3D) VectorFunction, creates the curl VectorFunction.
 
@@ -459,7 +455,6 @@ def curl(self, shift=None, order=None, method='FD', **kwargs):
         if len(self.space_dimensions) != 3:
             raise AttributeError("Curl only supported for 3D VectorFunction")
         # The curl of a VectorFunction is a VectorFunction
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         dims = self.root_dimensions
         derivs = ['d%s' % d.name for d in dims]
@@ -486,7 +481,7 @@ def curl(self, shift=None, order=None, method='FD', **kwargs):
         func = vec_func(self)
         return func._new(3, 1, [comp1, comp2, comp3])
 
-    def grad(self, shift=None, order=None, method='FD', **kwargs):
+    def grad(self, shift=None, order=None, method='FD', side=None, **kwargs):
         """
         Gradient of the VectorFunction, creates the gradient TensorFunction.
 
@@ -506,7 +501,6 @@ def grad(self, shift=None, order=None, method='FD', **kwargs):
         weights/w: list, tuple, or dict, optional, default=None
             Custom weights for the finite difference coefficients.
         """
-        side = kwargs.get("side")
         w = kwargs.get('weights', kwargs.get('w'))
         func = tens_func(self)
         ndim = len(self.space_dimensions)
diff --git a/tests/test_derivatives.py b/tests/test_derivatives.py
index c76d4eb22d..3d654d7a84 100644
--- a/tests/test_derivatives.py
+++ b/tests/test_derivatives.py
@@ -625,13 +625,13 @@ def test_shifted_grad(self, shift, ndim):
                 gk = getattr(f, 'd%s' % d.name)(x0=x0, fd_order=order).evaluate
                 assert gi == gk
 
-    # TODO: Repeat all these with vector and tensor functions where defined
     @pytest.mark.parametrize('side', [left, right, centered])
     def test_grad_w_side(self, side):
         grid = Grid(shape=(11, 11))
         f = Function(name='f', grid=grid, space_order=2)
 
-        # Want to check that it's the same as constructing by hand
+        # Want to check that it's the same as constructing by hand. Note that
+        # all the subsequent tests of this flavor work in the same way.
         comps = (f.dx(side=side), f.dy(side=side))
         expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
                                components=comps, grid=grid).evaluate
@@ -663,6 +663,21 @@ def test_div_w_side(self, side):
         assert expr1 == f.div(side=side).evaluate
         assert expr1 == div(f, side=side).evaluate
 
+    @pytest.mark.parametrize('side', [left, right, centered])
+    def test_tensor_div_w_side(self, side):
+        grid = Grid(shape=(11, 11))
+        f = TensorFunction(name='f', grid=grid, space_order=2,
+                           staggered=((None, None), (None, None)))
+
+        comps = (f[0, 0].dx(side=side) + f[0, 1].dy(side=side),
+                 f[0, 1].dx(side=side) + f[1, 1].dy(side=side))
+
+        expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
+                               components=comps, grid=grid).evaluate
+
+        assert expr1 == f.div(side=side).evaluate
+        assert expr1 == div(f, side=side).evaluate
+
     @pytest.mark.parametrize('side', [left, right, centered])
     def test_curl_w_side(self, side):
         grid = Grid(shape=(11, 11, 11))
@@ -689,6 +704,35 @@ def test_laplace_w_side(self, side):
         assert expr1 == f.laplacian(side=side).evaluate
         assert expr1 == laplace(f, side=side).evaluate
 
+    @pytest.mark.parametrize('side', [left, right, centered])
+    def test_vector_laplace_w_side(self, side):
+        grid = Grid(shape=(11, 11))
+        f = VectorFunction(name='f', grid=grid, space_order=2, staggered=(None, None))
+
+        comps = (f[0].dx2(side=side) + f[0].dy2(side=side),
+                 f[1].dx2(side=side) + f[1].dy2(side=side))
+
+        expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
+                               components=comps, grid=grid).evaluate
+
+        assert expr1 == f.laplacian(side=side).evaluate
+        assert expr1 == laplace(f, side=side).evaluate
+
+    @pytest.mark.parametrize('side', [centered])
+    def test_tensor_laplace_w_side(self, side):
+        grid = Grid(shape=(11, 11))
+        f = TensorFunction(name='f', grid=grid, space_order=2,
+                           staggered=((None, None), (None, None)))
+
+        comps = (f[0, 0].dx2(side=side) + f[0, 1].dy2(side=side),
+                 f[0, 1].dx2(side=side) + f[1, 1].dy2(side=side))
+
+        expr1 = VectorFunction(name=f"{f.name}_vec", space_order=f.space_order,
+                               components=comps, grid=grid).evaluate
+
+        assert expr1 == f.laplacian(side=side).evaluate
+        assert expr1 == laplace(f, side=side).evaluate
+
     def test_substitution(self):
         grid = Grid((11, 11))
         f = Function(name="f", grid=grid, space_order=4)

From 575d306b0a63e1dca5b1f3f8ad39205d4245c1d9 Mon Sep 17 00:00:00 2001
From: Edward Caunt <ed@devitocodes.com>
Date: Mon, 22 Dec 2025 10:52:19 +0000
Subject: [PATCH 4/4] api: Define side explicitly as a kwarg in operators

---
 devito/finite_differences/operators.py | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/devito/finite_differences/operators.py b/devito/finite_differences/operators.py
index 330454a7fd..f3f159f4ad 100644
--- a/devito/finite_differences/operators.py
+++ b/devito/finite_differences/operators.py
@@ -1,4 +1,4 @@
-def div(func, shift=None, order=None, method='FD', **kwargs):
+def div(func, shift=None, order=None, method='FD', side=None, **kwargs):
     """
     Divergence of the input Function.
 
@@ -20,7 +20,6 @@ def div(func, shift=None, order=None, method='FD', **kwargs):
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
-    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
         return func.div(shift=shift, order=order, method=method, side=side, w=w)
@@ -45,7 +44,7 @@ def div45(func, shift=None, order=None):
     return div(func, shift=shift, order=order, method='RSFD')
 
 
-def grad(func, shift=None, order=None, method='FD', **kwargs):
+def grad(func, shift=None, order=None, method='FD', side=None, **kwargs):
     """
     Gradient of the input Function.
 
@@ -67,7 +66,6 @@ def grad(func, shift=None, order=None, method='FD', **kwargs):
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
-    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
         return func.grad(shift=shift, order=order, method=method, side=side, w=w)
@@ -92,7 +90,7 @@ def grad45(func, shift=None, order=None):
     return grad(func, shift=shift, order=order, method='RSFD')
 
 
-def curl(func, shift=None, order=None, method='FD', **kwargs):
+def curl(func, shift=None, order=None, method='FD', side=None, **kwargs):
     """
     Curl of the input Function. Only supported for VectorFunction
 
@@ -114,7 +112,6 @@ def curl(func, shift=None, order=None, method='FD', **kwargs):
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
-    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
         return func.curl(shift=shift, order=order, method=method, side=side, w=w)
@@ -140,7 +137,7 @@ def curl45(func, shift=None, order=None):
     return curl(func, shift=shift, order=order, method='RSFD')
 
 
-def laplace(func, shift=None, order=None, method='FD', **kwargs):
+def laplace(func, shift=None, order=None, method='FD', side=None, **kwargs):
     """
     Laplacian of the input Function.
 
@@ -161,7 +158,6 @@ def laplace(func, shift=None, order=None, method='FD', **kwargs):
     weights/w: list, tuple, or dict, optional, default=None
         Custom weights for the finite difference coefficients.
     """
-    side = kwargs.get("side")
     w = kwargs.get('weights', kwargs.get('w'))
     try:
         return func.laplacian(shift=shift, order=order, method=method, side=side, w=w)