examples: Update tti fd for better stability

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

@@ -51,8 +51,8 @@ jobs:
     - name: Test examples with MPI
       run: |
         python3 scripts/clear_devito_cache.py
-        DEVITO_MPI=1 mpirun -n 2 python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml examples/seismic/acoustic
-        DEVITO_MPI=1 mpirun -n 2 python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml examples/seismic/tti
+        DEVITO_MPI=1 mpirun -n 2 python3 -m pytest examples/seismic/acoustic
+        DEVITO_MPI=1 mpirun -n 2 python3 -m pytest examples/seismic/tti
 
     - name: Upload coverage to Codecov
       uses: codecov/codecov-action@v4

examples/seismic/tti/operators.py

@@ -74,18 +74,22 @@ def Gzz_centered(model, field):
     """
     b = getattr(model, 'b', 1)
     costheta, sintheta, cosphi, sinphi = trig_func(model)
+
     order1 = field.space_order // 2
-    Gz = -(sintheta * cosphi * field.dx(fd_order=order1) +
-           sintheta * sinphi * field.dy(fd_order=order1) +
-           costheta * field.dz(fd_order=order1))
+    x, y, z = field.grid.dimensions
+    dx, dy, dz = x.spacing/2, y.spacing/2, z.spacing/2
+
+    Gz = (sintheta * cosphi * field.dx(fd_order=order1, x0=x+dx) +
+          sintheta * sinphi * field.dy(fd_order=order1, x0=y+dy) +
+          costheta * field.dz(fd_order=order1, x0=z+dz))
 
-    Gzz = (b * Gz * costheta).dz(fd_order=order1).T
+    Gzz = (b * Gz * costheta).dz(fd_order=order1, x0=z-dz)
     # Add rotated derivative if angles are not zero. If angles are
     # zeros then `0*Gz = 0` and doesn't have any `.dy` ....
     if sintheta != 0:
-        Gzz += (b * Gz * sintheta * cosphi).dx(fd_order=order1).T
+        Gzz += (b * Gz * sintheta * cosphi).dx(fd_order=order1, x0=x-dx)
     if sinphi != 0:
-        Gzz += (b * Gz * sintheta * sinphi).dy(fd_order=order1).T
+        Gzz += (b * Gz * sintheta * sinphi).dy(fd_order=order1, x0=y-dy)
 
     return Gzz
 
@@ -105,17 +109,21 @@ def Gzz_centered_2d(model, field):
     -------
     Rotated second order derivative w.r.t. z.
     """
+    b = getattr(model, 'b', 1)
     costheta, sintheta = trig_func(model)
+
     order1 = field.space_order // 2
-    b = getattr(model, 'b', 1)
-    Gz = -(sintheta * field.dx(fd_order=order1) +
-           costheta * field.dy(fd_order=order1))
-    Gzz = (b * Gz * costheta).dy(fd_order=order1).T
+    x, y = field.grid.dimensions
+    dx, dy = x.spacing/2, y.spacing/2
+
+    Gz = (sintheta * field.dx(fd_order=order1, x0=x+dx) +
+          costheta * field.dy(fd_order=order1, x0=y+dy))
+    Gzz = (b * Gz * costheta).dy(fd_order=order1, x0=y-dy)
 
     # Add rotated derivative if angles are not zero. If angles are
     # zeros then `0*Gz = 0` and doesn't have any `.dy` ....
     if sintheta != 0:
-        Gzz += (b * Gz * sintheta).dx(fd_order=order1).T
+        Gzz += (b * Gz * sintheta).dx(fd_order=order1, x0=x-dx)
     return Gzz
 
 

tests/test_dse.py

@@ -2826,7 +2826,7 @@ def test_opcounts(self, space_order, expected):
 
     @switchconfig(profiling='advanced')
     @pytest.mark.parametrize('space_order,exp_ops,exp_arrays', [
-        (4, 122, 6), (8, 235, 7)
+        (4, 122, 6), (8, 225, 7)
     ])
     def test_opcounts_adjoint(self, space_order, exp_ops, exp_arrays):
         wavesolver = self.tti_operator(space_order=space_order,

tests/test_gradient.py

@@ -244,7 +244,7 @@ def initializer(data):
         ('OT2', (70, 80), (15., 15.), iso_setup, 2),
         ('sls', (70, 80), (20., 20.), vsc_setup, 2),
         ('sls', (70, 80), (20., 20.), vsc_setup, 1),
-        ('centered', (70, 80), (15., 15.), tti_setup, 2),
+        ('centered', (70, 80), (20., 20.), tti_setup, 2),
     ])
     @pytest.mark.parametrize('space_order', [4])
     @pytest.mark.parametrize('dtype', [np.float32, np.float64])