diff --git a/.github/workflows/pytest-core-nompi.yml b/.github/workflows/pytest-core-nompi.yml
index a843d127d91..f57c3575ca4 100644
--- a/.github/workflows/pytest-core-nompi.yml
+++ b/.github/workflows/pytest-core-nompi.yml
@@ -58,7 +58,7 @@ jobs:
           sympy: "1.10"
 
         - name: pytest-ubuntu-py39-gcc7-omp
-          python-version: '3.9'
+          python-version: '3.7'
           os: ubuntu-20.04
           arch: "gcc-7"
           language: "openmp"
diff --git a/.github/workflows/tutorials.yml b/.github/workflows/tutorials.yml
index 72bb3c67c82..30ddb25d763 100644
--- a/.github/workflows/tutorials.yml
+++ b/.github/workflows/tutorials.yml
@@ -29,29 +29,29 @@ jobs:
       fail-fast: false
       matrix:
         name: [
-          tutos-ubuntu-gcc-py37,
-          tutos-osx-gcc-py37,
-          tutos-osx-clang-py37,
-          tutos-docker-gcc-py37
+          tutos-ubuntu-gcc-py38,
+          tutos-osx-gcc-py38,
+          tutos-osx-clang-py38,
+          tutos-docker-gcc-py38
           ]
 
         include:
-          - name: tutos-ubuntu-gcc-py37
+          - name: tutos-ubuntu-gcc-py38
             os: ubuntu-latest
             compiler: gcc
             language: "openmp"
 
-          - name: tutos-osx-gcc-py37
+          - name: tutos-osx-gcc-py38
             os: macos-latest
             compiler: gcc-11
             language: "openmp"
 
-          - name: tutos-osx-clang-py37
+          - name: tutos-osx-clang-py38
             os: macos-latest
             compiler: clang
             language: "C"
 
-          - name: tutos-docker-gcc-py37
+          - name: tutos-docker-gcc-py38
             os: ubuntu-latest
             compiler: gcc
             language: "openmp"
@@ -72,13 +72,13 @@ jobs:
         xcode-version: latest-stable
 
     - name: Build docker image
-      if: matrix.name == 'tutos-docker-gcc-py37'
+      if: matrix.name == 'tutos-docker-gcc-py38'
       run: |
           docker build . --file docker/Dockerfile.devito --tag devito_img
 
     - name: Set run prefix
       run: |
-          if [ "${{ matrix.name }}" == 'tutos-docker-gcc-py37' ]; then
+          if [ "${{ matrix.name }}" == 'tutos-docker-gcc-py38' ]; then
               echo "RUN_CMD=docker run --rm --name testrun devito_img" >> $GITHUB_ENV
           else
               echo "RUN_CMD=" >> $GITHUB_ENV
@@ -86,7 +86,7 @@ jobs:
       id: set-run
 
     - name: Install dependencies
-      if: matrix.name != 'tutos-docker-gcc-py37'
+      if: matrix.name != 'tutos-docker-gcc-py38'
       run: |
         python -m pip install --upgrade pip
         pip install -e .[tests]
diff --git a/examples/seismic/tutorials/05_staggered_acoustic.ipynb b/examples/seismic/tutorials/05_staggered_acoustic.ipynb
index 2b04a5f08cc..b5ff3e8c523 100644
--- a/examples/seismic/tutorials/05_staggered_acoustic.ipynb
+++ b/examples/seismic/tutorials/05_staggered_acoustic.ipynb
@@ -122,7 +122,7 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\left[\\begin{matrix}\\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},z \\right)}\\\\\\operatorname{v_{z}}{\\left(t + dt,x,z + \\frac{h_{z}}{2} \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(\\frac{\\partial}{\\partial x} p{\\left(t,x,z \\right)} + \\frac{\\operatorname{v_{x}}{\\left(t,x + \\frac{h_{x}}{2},z \\right)}}{dt}\\right)\\\\dt \\left(\\frac{\\partial}{\\partial z} p{\\left(t,x,z \\right)} + \\frac{\\operatorname{v_{z}}{\\left(t,x,z + \\frac{h_{z}}{2} \\right)}}{dt}\\right)\\end{matrix}\\right]$"
+       "$\\displaystyle \\left[\\begin{matrix}v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},z \\right)}\\\\v_{z}{\\left(t + dt,x,z + \\frac{h_{z}}{2} \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(\\frac{\\partial}{\\partial x} p{\\left(t,x,z \\right)} + \\frac{v_{x}{\\left(t,x + \\frac{h_{x}}{2},z \\right)}}{dt}\\right)\\\\dt \\left(\\frac{\\partial}{\\partial z} p{\\left(t,x,z \\right)} + \\frac{v_{z}{\\left(t,x,z + \\frac{h_{z}}{2} \\right)}}{dt}\\right)\\end{matrix}\\right]$"
       ],
       "text/plain": [
        "Eq(Vector(v_x(t + dt, x + h_x/2, z), v_z(t + dt, x, z + h_z/2)), Vector(dt*(Derivative(p(t, x, z), x) + v_x(t, x + h_x/2, z)/dt), dt*(Derivative(p(t, x, z), z) + v_z(t, x, z + h_z/2)/dt)))"
@@ -145,7 +145,7 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle p{\\left(t + dt,x,z \\right)} = dt \\left(16.0 \\frac{\\partial}{\\partial x} \\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},z \\right)} + 16.0 \\frac{\\partial}{\\partial z} \\operatorname{v_{z}}{\\left(t + dt,x,z + \\frac{h_{z}}{2} \\right)} + \\frac{p{\\left(t,x,z \\right)}}{dt}\\right)$"
+       "$\\displaystyle p{\\left(t + dt,x,z \\right)} = dt \\left(16.0 \\frac{\\partial}{\\partial x} v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},z \\right)} + 16.0 \\frac{\\partial}{\\partial z} v_{z}{\\left(t + dt,x,z + \\frac{h_{z}}{2} \\right)} + \\frac{p{\\left(t,x,z \\right)}}{dt}\\right)$"
       ],
       "text/plain": [
        "Eq(p(t + dt, x, z), dt*(16.0*Derivative(v_x(t + dt, x + h_x/2, z), x) + 16.0*Derivative(v_z(t + dt, x, z + h_z/2), z) + p(t, x, z)/dt))"
diff --git a/examples/seismic/tutorials/06_elastic.ipynb b/examples/seismic/tutorials/06_elastic.ipynb
index cec0767a43f..210f9936a34 100644
--- a/examples/seismic/tutorials/06_elastic.ipynb
+++ b/examples/seismic/tutorials/06_elastic.ipynb
@@ -153,7 +153,7 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\left[\\begin{matrix}\\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},z \\right)}\\\\\\operatorname{v_{z}}{\\left(t + dt,x,z + \\frac{h_{z}}{2} \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(0.555555555555556 \\frac{\\partial}{\\partial x} \\operatorname{t_{xx}}{\\left(t,x,z \\right)} + 0.555555555555556 \\frac{\\partial}{\\partial z} \\operatorname{t_{xz}}{\\left(t,x + \\frac{h_{x}}{2},z + \\frac{h_{z}}{2} \\right)} + \\frac{\\operatorname{v_{x}}{\\left(t,x + \\frac{h_{x}}{2},z \\right)}}{dt}\\right)\\\\dt \\left(0.555555555555556 \\frac{\\partial}{\\partial x} \\operatorname{t_{xz}}{\\left(t,x + \\frac{h_{x}}{2},z + \\frac{h_{z}}{2} \\right)} + 0.555555555555556 \\frac{\\partial}{\\partial z} \\operatorname{t_{zz}}{\\left(t,x,z \\right)} + \\frac{\\operatorname{v_{z}}{\\left(t,x,z + \\frac{h_{z}}{2} \\right)}}{dt}\\right)\\end{matrix}\\right]$"
+       "$\\displaystyle \\left[\\begin{matrix}v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},z \\right)}\\\\v_{z}{\\left(t + dt,x,z + \\frac{h_{z}}{2} \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(0.555555555555556 \\frac{\\partial}{\\partial x} t_{xx}{\\left(t,x,z \\right)} + 0.555555555555556 \\frac{\\partial}{\\partial z} t_{xz}{\\left(t,x + \\frac{h_{x}}{2},z + \\frac{h_{z}}{2} \\right)} + \\frac{v_{x}{\\left(t,x + \\frac{h_{x}}{2},z \\right)}}{dt}\\right)\\\\dt \\left(0.555555555555556 \\frac{\\partial}{\\partial x} t_{xz}{\\left(t,x + \\frac{h_{x}}{2},z + \\frac{h_{z}}{2} \\right)} + 0.555555555555556 \\frac{\\partial}{\\partial z} t_{zz}{\\left(t,x,z \\right)} + \\frac{v_{z}{\\left(t,x,z + \\frac{h_{z}}{2} \\right)}}{dt}\\right)\\end{matrix}\\right]$"
       ],
       "text/plain": [
        "Eq(Vector(v_x(t + dt, x + h_x/2, z), v_z(t + dt, x, z + h_z/2)), Vector(dt*(0.555555555555556*Derivative(t_xx(t, x, z), x) + 0.555555555555556*Derivative(t_xz(t, x + h_x/2, z + h_z/2), z) + v_x(t, x + h_x/2, z)/dt), dt*(0.555555555555556*Derivative(t_xz(t, x + h_x/2, z + h_z/2), x) + 0.555555555555556*Derivative(t_zz(t, x, z), z) + v_z(t, x, z + h_z/2)/dt)))"
diff --git a/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb b/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb
index f5c3c3f1062..d6d13331f0d 100644
--- a/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb
+++ b/examples/seismic/tutorials/06_elastic_varying_parameters.ipynb
@@ -353,7 +353,7 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\left[\\begin{matrix}\\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)}\\\\\\operatorname{v_{y}}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(\\left(\\frac{\\partial}{\\partial x} \\operatorname{t_{xx}}{\\left(t,x,y \\right)} + \\frac{\\partial}{\\partial y} \\operatorname{t_{xy}}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}\\right) b{\\left(x,y \\right)} + \\frac{\\operatorname{v_{x}}{\\left(t,x + \\frac{h_{x}}{2},y \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\\\dt \\left(\\left(\\frac{\\partial}{\\partial x} \\operatorname{t_{xy}}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)} + \\frac{\\partial}{\\partial y} \\operatorname{t_{yy}}{\\left(t,x,y \\right)}\\right) b{\\left(x,y \\right)} + \\frac{\\operatorname{v_{y}}{\\left(t,x,y + \\frac{h_{y}}{2} \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\end{matrix}\\right]$"
+       "$\\displaystyle \\left[\\begin{matrix}v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)}\\\\v_{y}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(\\left(\\frac{\\partial}{\\partial x} t_{xx}{\\left(t,x,y \\right)} + \\frac{\\partial}{\\partial y} t_{xy}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}\\right) b{\\left(x,y \\right)} + \\frac{v_{x}{\\left(t,x + \\frac{h_{x}}{2},y \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\\\dt \\left(\\left(\\frac{\\partial}{\\partial x} t_{xy}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)} + \\frac{\\partial}{\\partial y} t_{yy}{\\left(t,x,y \\right)}\\right) b{\\left(x,y \\right)} + \\frac{v_{y}{\\left(t,x,y + \\frac{h_{y}}{2} \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\end{matrix}\\right]$"
       ],
       "text/plain": [
        "Eq(Vector(v_x(t + dt, x + h_x/2, y), v_y(t + dt, x, y + h_y/2)), Vector(dt*((Derivative(t_xx(t, x, y), x) + Derivative(t_xy(t, x + h_x/2, y + h_y/2), y))*b(x, y) + v_x(t, x + h_x/2, y)/dt)*damp(x, y), dt*((Derivative(t_xy(t, x + h_x/2, y + h_y/2), x) + Derivative(t_yy(t, x, y), y))*b(x, y) + v_y(t, x, y + h_y/2)/dt)*damp(x, y)))"
@@ -376,7 +376,7 @@
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\left[\\begin{matrix}\\operatorname{t_{xx}}{\\left(t + dt,x,y \\right)} & \\operatorname{t_{xy}}{\\left(t + dt,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}\\\\\\operatorname{t_{xy}}{\\left(t + dt,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)} & \\operatorname{t_{yy}}{\\left(t + dt,x,y \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(\\left(\\frac{\\partial}{\\partial x} \\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial y} \\operatorname{v_{y}}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\operatorname{lam}{\\left(x,y \\right)} + 2 \\mu{\\left(x,y \\right)} \\frac{\\partial}{\\partial x} \\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\operatorname{t_{xx}}{\\left(t,x,y \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)} & dt \\left(\\left(\\frac{\\partial}{\\partial y} \\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial x} \\operatorname{v_{y}}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\mu{\\left(x,y \\right)} + \\frac{\\operatorname{t_{xy}}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\\\dt \\left(\\left(\\frac{\\partial}{\\partial y} \\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial x} \\operatorname{v_{y}}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\mu{\\left(x,y \\right)} + \\frac{\\operatorname{t_{xy}}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)} & dt \\left(\\left(\\frac{\\partial}{\\partial x} \\operatorname{v_{x}}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial y} \\operatorname{v_{y}}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\operatorname{lam}{\\left(x,y \\right)} + 2 \\mu{\\left(x,y \\right)} \\frac{\\partial}{\\partial y} \\operatorname{v_{y}}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)} + \\frac{\\operatorname{t_{yy}}{\\left(t,x,y \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\end{matrix}\\right]$"
+       "$\\displaystyle \\left[\\begin{matrix}t_{xx}{\\left(t + dt,x,y \\right)} & t_{xy}{\\left(t + dt,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}\\\\t_{xy}{\\left(t + dt,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)} & t_{yy}{\\left(t + dt,x,y \\right)}\\end{matrix}\\right] = \\left[\\begin{matrix}dt \\left(\\left(\\frac{\\partial}{\\partial x} v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial y} v_{y}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\operatorname{lam}{\\left(x,y \\right)} + 2 \\mu{\\left(x,y \\right)} \\frac{\\partial}{\\partial x} v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{t_{xx}{\\left(t,x,y \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)} & dt \\left(\\left(\\frac{\\partial}{\\partial y} v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial x} v_{y}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\mu{\\left(x,y \\right)} + \\frac{t_{xy}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\\\dt \\left(\\left(\\frac{\\partial}{\\partial y} v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial x} v_{y}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\mu{\\left(x,y \\right)} + \\frac{t_{xy}{\\left(t,x + \\frac{h_{x}}{2},y + \\frac{h_{y}}{2} \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)} & dt \\left(\\left(\\frac{\\partial}{\\partial x} v_{x}{\\left(t + dt,x + \\frac{h_{x}}{2},y \\right)} + \\frac{\\partial}{\\partial y} v_{y}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)}\\right) \\operatorname{lam}{\\left(x,y \\right)} + 2 \\mu{\\left(x,y \\right)} \\frac{\\partial}{\\partial y} v_{y}{\\left(t + dt,x,y + \\frac{h_{y}}{2} \\right)} + \\frac{t_{yy}{\\left(t,x,y \\right)}}{dt}\\right) \\operatorname{damp}{\\left(x,y \\right)}\\end{matrix}\\right]$"
       ],
       "text/plain": [
        "Eq(Matrix([[t_xx(t + dt, x, y), t_xy(t + dt, x + h_x/2, y + h_y/2)], [t_xy(t + dt, x + h_x/2, y + h_y/2), t_yy(t + dt, x, y)]]), Matrix([[dt*((Derivative(v_x(t + dt, x + h_x/2, y), x) + Derivative(v_y(t + dt, x, y + h_y/2), y))*lam(x, y) + 2*mu(x, y)*Derivative(v_x(t + dt, x + h_x/2, y), x) + t_xx(t, x, y)/dt)*damp(x, y), dt*((Derivative(v_x(t + dt, x + h_x/2, y), y) + Derivative(v_y(t + dt, x, y + h_y/2), x))*mu(x, y) + t_xy(t, x + h_x/2, y + h_y/2)/dt)*damp(x, y)], [dt*((Derivative(v_x(t + dt, x + h_x/2, y), y) + Derivative(v_y(t + dt, x, y + h_y/2), x))*mu(x, y) + t_xy(t, x + h_x/2, y + h_y/2)/dt)*damp(x, y), dt*((Derivative(v_x(t + dt, x + h_x/2, y), x) + Derivative(v_y(t + dt, x, y + h_y/2), y))*lam(x, y) + 2*mu(x, y)*Derivative(v_y(t + dt, x, y + h_y/2), y) + t_yy(t, x, y)/dt)*damp(x, y)]]))"