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-latest
           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 func is self.expr:
+            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/tools.py

@@ -284,12 +284,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 expr.is_Staggered 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/elastic/elastic_example.py

@@ -37,7 +37,7 @@ def run(shape=(50, 50), spacing=(20.0, 20.0), tn=1000.0,
 def test_elastic(dtype):
     _, _, _, [rec1, rec2, v, tau] = run(dtype=dtype)
     assert np.isclose(norm(rec1), 19.25636, atol=1e-3, rtol=0)
-    assert np.isclose(norm(rec2), 0.6276, atol=1e-3, rtol=0)
+    assert np.isclose(norm(rec2), 0.6444, atol=1e-3, rtol=0)
 
 
 @pytest.mark.parametrize('shape', [(51, 51), (16, 16, 16)])

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.20283400000000001, 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.005394999999999986, 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.006213000000000028, 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.0058910000000000256, 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.005846000000000029, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 14,
@@ -518,15 +518,15 @@
      "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.20509600000000006, 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.006368999999999996, 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.007619000000000035, 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.006970000000000019, 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.007279000000000034, 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.25 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.21811100000000017, 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.007676999999999997, 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.008889000000000003, 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.008588000000000002, 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.29 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.25438299999999997, 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.009171000000000014, 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.009659999999999995, 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.009119000000000006, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])"
       ]
      },
      "execution_count": 26,
@@ -921,12 +921,12 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\frac{f(x, y)}{h_{x}} - \\frac{f(x - h_x, y)}{h_{x}}$"
+       "$\\displaystyle - \\frac{f(x, y)}{h_{x}} + \\frac{f(x + h_x, y)}{h_{x}}$"
       ],
       "text/plain": [
-       "f(x, y)   f(x - hₓ, y)\n",
-       "─────── - ────────────\n",
-       "   hₓ          hₓ     "
+       "  f(x, y)   f(x + hₓ, y)\n",
+       "- ─────── + ────────────\n",
+       "     hₓ          hₓ     "
       ]
      },
      "execution_count": 31,
@@ -946,16 +946,12 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\left(- \\frac{f(x - h_x, y - h_y)}{4 h_{x}} + \\frac{f(x + h_x, y + h_y)}{4 h_{x}}\\right) + \\left(- \\frac{f(x - h_x, y + h_y)}{4 h_{x}} + \\frac{f(x + h_x, y - h_y)}{4 h_{x}}\\right)$"
+       "$\\displaystyle \\left(- \\frac{f(x, y)}{2 h_{x}} + \\frac{f(x + h_x, y - h_y)}{2 h_{x}}\\right) + \\left(- \\frac{f(x, y)}{2 h_{x}} + \\frac{f(x + h_x, y + h_y)}{2 h_{x}}\\right)$"
       ],
       "text/plain": [
-       "  f(x - hₓ, y - h_y)   f(x + hₓ, y + h_y)     f(x - hₓ, y + h_y)   f(x + hₓ, y\n",
-       "- ────────────────── + ────────────────── + - ────────────────── + ───────────\n",
-       "         4⋅hₓ                 4⋅hₓ                   4⋅hₓ                 4⋅hₓ\n",
-       "\n",
-       " - h_y)\n",
-       "───────\n",
-       "       "
+       "  f(x, y)   f(x + hₓ, y - h_y)     f(x, y)   f(x + hₓ, y + h_y)\n",
+       "- ─────── + ────────────────── + - ─────── + ──────────────────\n",
+       "    2⋅hₓ           2⋅hₓ              2⋅hₓ           2⋅hₓ       "
       ]
      },
      "execution_count": 32,
@@ -1023,15 +1019,15 @@
      "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.524889000000001, 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.010354000000000039, 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.011986999999999972, 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.011274000000000013, 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.011272999999999995, 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]