api: fix corner case staggered fd for centered x0

.github/workflows/pytest-core-nompi.yml

@@ -72,7 +72,7 @@ jobs:
 
         - name: pytest-ubuntu-py312-gcc13-omp
           python-version: '3.12'
-          os: ubuntu-20.04
+          os: ubuntu-24.04
           arch: "gcc-13"
           language: "openmp"
           sympy: "1.11"

devito/finite_differences/derivative.py

@@ -356,8 +356,12 @@ def _eval_at(self, func):
         setup where one could have Eq(u(x + h_x/2), v(x).dx)) in which case v(x).dx
         has to be computed at x=x + h_x/2.
         """
-        # If an x0 already exists do not overwrite it
-        x0 = self.x0 or func.indices_ref._getters
+        # If an x0 already exists or evaluating at the same function (i.e u = u.dx)
+        # do not overwrite it
+        if self.x0 or self.side is not None or func.function is self.expr.function:
+            return self
+
+        x0 = func.indices_ref._getters
         if self.expr.is_Add:
             # If `expr` has both staggered and non-staggered terms such as
             # `(u(x + h_x/2) + v(x)).dx` then we exploit linearity of FD to split

devito/finite_differences/differentiable.py

@@ -121,6 +121,13 @@ def _fd(self):
     def _symbolic_functions(self):
         return frozenset([i for i in self._functions if i.coefficients == 'symbolic'])
 
+    @cached_property
+    def function(self):
+        if len(self._functions) == 1:
+            return self._functions.pop()
+        else:
+            return None
+
     @cached_property
     def _uses_symbolic_coefficients(self):
         return bool(self._symbolic_functions)

devito/finite_differences/tools.py

@@ -275,6 +275,7 @@ def generate_indices(expr, dim, order, side=None, matvec=None, x0=None):
     -------
     An IndexSet, representing an ordered list of indices.
     """
+    stagg_fd = expr.is_Staggered and x0
     # Evaluation point
     x0 = sympify(((x0 or {}).get(dim) or expr.indices_ref[dim]))
 
@@ -284,12 +285,17 @@ def generate_indices(expr, dim, order, side=None, matvec=None, x0=None):
         iexpr = x0 + d * dim.spacing
         return IndexSet(dim, expr=iexpr), x0
 
-    # Shift for side
-    side = side or centered
-
     # Evaluation point relative to the expression's grid
     mid = (x0 - expr.indices_ref[dim]).subs({dim: 0, dim.spacing: 1})
 
+    # Centered scheme for staggered field create artifacts if not shifted
+    if stagg_fd and mid == 0 and not dim.is_Time and side is not centered:
+        mid = 0.5 if x0 is dim else -0.5
+        x0 = x0 + mid * dim.spacing
+
+    # Shift for side
+    side = side or centered
+
     # Indices range
     o_min = int(np.ceil(mid - order/2)) + side.val
     o_max = int(np.floor(mid + order/2)) + side.val

examples/seismic/tti/operators.py

@@ -299,10 +299,10 @@ def kernel_staggered_2d(model, u, v, **kwargs):
     else:
         # Stencils
         phdx = ((costheta*epsilon*u).dx - (sintheta*epsilon*u).dyc +
-                (costheta*delta*v).dxc - (sintheta*delta*v).dy)
+                (costheta*delta*v).dx - (sintheta*delta*v).dyc)
         u_vx = Eq(vx.backward, dampl * vx + dampl * s * phdx)
 
-        pvdz = ((sintheta*delta*u).dx + (costheta*delta*u).dyc +
+        pvdz = ((sintheta*delta*u).dxc + (costheta*delta*u).dy +
                 (sintheta*v).dxc + (costheta*v).dy)
         u_vz = Eq(vz.backward, dampl * vz + dampl * s * pvdz)
 

examples/seismic/tutorials/06_elastic_varying_parameters.ipynb

@@ -392,22 +392,22 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.22 s\n"
+      "Operator `Kernel` ran in 0.23 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.19680499999999987, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.20279700000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.005080999999999991, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.005569999999999994, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.005789000000000039, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.006224000000000026, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.005510000000000015, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.006004000000000021, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section4', rank=None),\n",
-       "  PerfEntry(time=0.005339000000000027, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.00578600000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 14,
@@ -511,22 +511,22 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.24 s\n"
+      "Operator `Kernel` ran in 0.23 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.20575800000000005, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.20478399999999966, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.006976999999999987, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.005923000000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.007339000000000036, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.00643000000000003, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.0069080000000000166, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.006100000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section4', rank=None),\n",
-       "  PerfEntry(time=0.006813000000000026, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.006256000000000025, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 17,
@@ -699,20 +699,20 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.21 s\n"
+      "Operator `Kernel` ran in 0.26 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.1906600000000002, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.22333499999999995, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.005354999999999988, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.008422999999999998, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.005681000000000019, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.009065000000000017, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.005547000000000023, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.009086000000000014, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 24,
@@ -796,20 +796,20 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.27 s\n"
+      "Operator `Kernel` ran in 0.22 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.22429000000000004, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.19835300000000017, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.013601000000000016, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.005520999999999994, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.012073000000000028, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.005916000000000023, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.012027000000000043, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.006079000000000015, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 26,
@@ -1016,22 +1016,22 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Operator `Kernel` ran in 0.57 s\n"
+      "Operator `Kernel` ran in 0.71 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
        "PerformanceSummary([(PerfKey(name='section0', rank=None),\n",
-       "  PerfEntry(time=0.5206050000000007, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.6343909999999997, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section1', rank=None),\n",
-       "  PerfEntry(time=0.010398000000000041, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.017495000000000028, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section2', rank=None),\n",
-       "  PerfEntry(time=0.012117999999999964, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.018152000000000022, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section3', rank=None),\n",
-       "  PerfEntry(time=0.011993000000000005, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
+       "  PerfEntry(time=0.017839, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),\n",
        " (PerfKey(name='section4', rank=None),\n",
-       "  PerfEntry(time=0.011515999999999974, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
+       "  PerfEntry(time=0.017724000000000007, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 34,

tests/test_derivatives.py

@@ -440,6 +440,20 @@ def test_new_x0_eval_at(self):
         v = Function(name="v", grid=grid, space_order=2)
         assert u.dx(x0=x - x.spacing/2)._eval_at(v).x0 == {x: x - x.spacing/2}
 
+    @pytest.mark.parametrize('stagg', [True, False])
+    def test_eval_at_centered(self, stagg):
+        grid = Grid((10,))
+        x = grid.dimensions[0]
+        stagg = NODE if stagg else x
+        x0 = x if stagg else x + .5 * x.spacing
+
+        u = Function(name="u", grid=grid, space_order=2, staggered=stagg)
+        v = Function(name="v", grid=grid, space_order=2, staggered=stagg)
+
+        assert u.dx._eval_at(v).evaluate == u.dx(x0=x0).evaluate
+        assert v.dx._eval_at(u).evaluate == v.dx(x0=x0).evaluate
+        assert u.dx._eval_at(u).evaluate == u.dx.evaluate
+
     def test_fd_new_lo(self):
         grid = Grid((10,))
         x = grid.dimensions[0]