From 4bf52f1ee879c3383b794b0a3daed6de7e9b0755 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 13:57:33 +0100
Subject: [PATCH 1/8] io tests

---
 tests/firedrake/output/test_io_function.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/firedrake/output/test_io_function.py b/tests/firedrake/output/test_io_function.py
index 9d440b5b31..bb1ae7616a 100644
--- a/tests/firedrake/output/test_io_function.py
+++ b/tests/firedrake/output/test_io_function.py
@@ -34,7 +34,7 @@ def _get_mesh(cell_type, comm):
         mesh = Mesh(mesh_file, name=mesh_name, comm=comm)
     elif cell_type == "tetrahedra":
         # TODO: Prepare more interesting mesh.
-        mesh = UnitCubeMesh(16, 16, 16, name=mesh_name, comm=comm)
+        mesh = UnitCubeMesh(4, 4, 4, name=mesh_name, comm=comm)
     elif cell_type == "tetrahedra_large":
         mesh = Mesh(join(os.environ["PETSC_DIR"], "share/petsc/datafiles/meshes/mesh-3d-box-innersphere.msh"),
                     name=mesh_name, comm=comm)
@@ -56,9 +56,9 @@ def _get_mesh(cell_type, comm):
         # Sanity check
         mesh = UnitSquareMesh(2, 2, quadrilateral=True, name=mesh_name)
     elif cell_type == "triangle_periodic":
-        mesh = PeriodicUnitSquareMesh(20, 20, name=mesh_name)
+        mesh = PeriodicUnitSquareMesh(10, 10, name=mesh_name)
     elif cell_type == "tetrahedra_periodic":
-        mesh = PeriodicUnitCubeMesh(10, 10, 10, name=mesh_name)
+        mesh = PeriodicUnitCubeMesh(4, 4, 4, name=mesh_name)
     elif cell_type == "triangle_3d":
         mesh = UnitIcosahedralSphereMesh(refinement_level=1, name=mesh_name)
         x = SpatialCoordinate(mesh)

From f369aae48745d7f6039d72556ceac5aa3f98a080 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 14:08:11 +0100
Subject: [PATCH 2/8] cross-mesh tests

---
 .../regression/test_interpolate_cross_mesh.py          | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/firedrake/regression/test_interpolate_cross_mesh.py b/tests/firedrake/regression/test_interpolate_cross_mesh.py
index a04e3edbc0..c7c56c9abf 100644
--- a/tests/firedrake/regression/test_interpolate_cross_mesh.py
+++ b/tests/firedrake/regression/test_interpolate_cross_mesh.py
@@ -199,8 +199,8 @@ def parameters(request):
         V_dest = TensorFunctionSpace(m_dest, "CG", 4)
         V_dest_2 = TensorFunctionSpace(m_dest, "DQ", 2)
     elif request.param == "spheresphere":
-        m_src = UnitCubedSphereMesh(5, name="src_sphere")
-        m_dest = UnitIcosahedralSphereMesh(5, name="dest_sphere")
+        m_src = UnitCubedSphereMesh(4, name="src_sphere")
+        m_dest = UnitIcosahedralSphereMesh(4, name="dest_sphere")
         coords = np.array(
             [
                 [0, 1, 0],
@@ -273,8 +273,8 @@ def test_interpolate_unitsquare_mixed():
     V_1 = FunctionSpace(m_src, "CG", 1)
     V_2 = FunctionSpace(m_src, "CG", 2)
     V_src = V_1 * V_2
-    V_3 = FunctionSpace(m_dest, "CG", 3)
-    V_4 = FunctionSpace(m_dest, "CG", 4)
+    V_3 = FunctionSpace(m_dest, "CG", 1)
+    V_4 = FunctionSpace(m_dest, "CG", 2)
     V_dest = V_3 * V_4
     f_src = Function(V_src)
     f_src.subfunctions[0].interpolate(expr_1)
@@ -510,7 +510,7 @@ def test_interpolate_cross_mesh(run_test, space, parameters):
     if m_src.name == "src_sphere" and m_dest.name == "dest_sphere":
         # Between immersed manifolds we will often be doing projection so we
         # need a higher tolerance for our tests
-        atol = 1e-3
+        atol = 5e-2
     else:
         atol = 1e-8  # default
     run_test(

From 8a27c6f10a26aaaaf20eca23e5a4f973c1fa9595 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 14:11:53 +0100
Subject: [PATCH 3/8] interp kmv

---
 tests/firedrake/regression/test_project_interp_KMV.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/firedrake/regression/test_project_interp_KMV.py b/tests/firedrake/regression/test_project_interp_KMV.py
index 7eef071841..096d222686 100644
--- a/tests/firedrake/regression/test_project_interp_KMV.py
+++ b/tests/firedrake/regression/test_project_interp_KMV.py
@@ -42,11 +42,11 @@ def run_interpolation(mesh, expr, p):
 def test_interpolation_KMV(mesh, max_degree, interpolation_expr):
     for p in range(1, max_degree):
         errors = [
-            run_interpolation(mesh(r), interpolation_expr, p) for r in range(3, 6)
+            run_interpolation(mesh(r), interpolation_expr, p) for r in range(2, 5)
         ]
         errors = np.asarray(errors)
         l2conv = np.log2(errors[:-1] / errors[1:])
-        assert (l2conv > p + 0.7).all()
+        assert (l2conv > p + 0.6).all()
 
 
 def run_projection(mesh, expr, p):

From b7b6b7a162068cb137233dfbb05edf256839b81c Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 14:32:25 +0100
Subject: [PATCH 4/8] test fdm

---
 tests/firedrake/regression/test_fdm.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/firedrake/regression/test_fdm.py b/tests/firedrake/regression/test_fdm.py
index 872b5eeb9b..60c77d09e0 100644
--- a/tests/firedrake/regression/test_fdm.py
+++ b/tests/firedrake/regression/test_fdm.py
@@ -119,7 +119,7 @@ def solve_riesz_map(problem, solver_parameters):
 @pytest.fixture(params=[2, 3],
                 ids=["Rectangle", "Box"])
 def mesh(request):
-    nx = 4
+    nx = 2
     distribution = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)}
     m = UnitSquareMesh(nx, nx, quadrilateral=True, distribution_parameters=distribution)
     if request.param == 3:

From c1171e9c89020e834746897810da960902ec709b Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 14:36:14 +0100
Subject: [PATCH 5/8] test external op

---
 .../external_operators/test_external_operators_adjoint.py     | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/firedrake/external_operators/test_external_operators_adjoint.py b/tests/firedrake/external_operators/test_external_operators_adjoint.py
index a563e55115..8b7ece1922 100644
--- a/tests/firedrake/external_operators/test_external_operators_adjoint.py
+++ b/tests/firedrake/external_operators/test_external_operators_adjoint.py
@@ -45,7 +45,7 @@ def assemble_Jacobian_adjoint_action(self, *args, **kwargs):
             y, _ = self.argument_slots()
             return y
 
-    mesh = UnitSquareMesh(50, 50)
+    mesh = UnitSquareMesh(10, 10)
     V = FunctionSpace(mesh, "CG", 1)
     x, y = SpatialCoordinate(mesh)
     u_exact = Function(V).interpolate(sin(pi * x) * sin(pi * y))
@@ -79,4 +79,4 @@ def J(f):
 
     f_opt = minimize(Jhat, tol=1e-4, method="BFGS")
 
-    assert assemble((f_exact - f_opt)**2 * dx) / assemble(f_exact**2 * dx) < 1e-5
+    assert assemble((f_exact - f_opt)**2 * dx) / assemble(f_exact**2 * dx) < 1e-4

From fe12e5c2b893d903cb458e88525a1ac61c5e2435 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 14:38:15 +0100
Subject: [PATCH 6/8] test poisson inverse

---
 tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py b/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
index 15b03e867c..96c25a3dfa 100644
--- a/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
+++ b/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
@@ -109,7 +109,7 @@ def test_poisson_inverse_conductivity(num_points):
     J_hat = ReducedFunctional(J, q_hat)
 
     # Estimate q using Newton-CG which evaluates the hessian action
-    minimize(J_hat, method='Newton-CG', options={'maxiter': 10, 'xtol': 1e-1, 'disp': True})
+    minimize(J_hat, method='Newton-CG', options={'maxiter': 2, 'xtol': 1e-1, 'disp': True})
 
 
 @pytest.mark.skipcomplex  # Taping for complex-valued 0-forms not yet done

From d165f0746709cce3a48dff274a4cac0f734d0734 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 14:48:24 +0100
Subject: [PATCH 7/8] vlasov poisson

---
 demos/vlasov_poisson_1d/vp1d.py.rst | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/demos/vlasov_poisson_1d/vp1d.py.rst b/demos/vlasov_poisson_1d/vp1d.py.rst
index be950f2ac3..f388375215 100644
--- a/demos/vlasov_poisson_1d/vp1d.py.rst
+++ b/demos/vlasov_poisson_1d/vp1d.py.rst
@@ -296,11 +296,21 @@ We then use this to build a solver. ::
 We are getting close to the time loop. We set up some timestepping
 parameters. ::
   
-  T = 50.0 # maximum timestep
-  t = 0. # model time
-  ndump = 100 # frequency of file dumps
-  dumpn = 0 # dump counter
-  nsteps = 5000
+  import os
+  if os.getenv("FIREDRAKE_CI") == "1":
+      # Setup for a faster test execution.
+      T = 0.5 # maximum timestep
+      t = 0. # model time
+      ndump = 50 # frequency of file dumps
+      dumpn = 0 # dump counter
+      nsteps = 50
+  else:
+      T = 50.0 # maximum timestep
+      t = 0. # model time
+      ndump = 100 # frequency of file dumps
+      dumpn = 0 # dump counter
+      nsteps = 5000
+
   dt = T/nsteps
   dtc.assign(dt)
 

From 599e391cb87db8dd42f3b2a4fcb6b01050a94984 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 6 Aug 2026 15:46:12 +0100
Subject: [PATCH 8/8] and the rest

---
 .../adaptive_multigrid.py.rst                 |  6 +++++
 demos/benney_luke/benney_luke.py.rst          |  6 +++++
 .../fast_diagonalisation_poisson.py.rst       | 25 +++++++++++++++---
 demos/ma-demo/ma-demo.py.rst                  |  6 +++++
 demos/matrix_free/rayleigh-benard.py.rst      |  5 ++++
 demos/netgen/netgen_mesh.py.rst               | 26 +++++++++++++++++--
 .../poisson_adaptive_eigenvalue.py.rst        |  6 +++++
 .../reynolds_robust_navier_stokes_hdiv.py.rst | 19 +++++++++++---
 .../saddle_point_systems.py.rst               | 26 +++++++++++++++++++
 docs/notebooks/07-geometric-multigrid.py      | 17 +++++++++---
 docs/notebooks/09-hybridisation.py            |  7 ++++-
 docs/notebooks/10-sum-factorisation.py        | 15 ++++++++++-
 docs/notebooks/12-HPC_demo.py                 |  6 +++++
 13 files changed, 156 insertions(+), 14 deletions(-)

diff --git a/demos/adaptive_multigrid/adaptive_multigrid.py.rst b/demos/adaptive_multigrid/adaptive_multigrid.py.rst
index 44ff4eac0f..126860835d 100644
--- a/demos/adaptive_multigrid/adaptive_multigrid.py.rst
+++ b/demos/adaptive_multigrid/adaptive_multigrid.py.rst
@@ -153,6 +153,12 @@ With these helper functions complete, we can solve the system iteratively. In th
 
    theta = 0.5
    refinements = 15
+
+   # Setup for faster test execution.
+   import os
+   if os.getenv("FIREDRAKE_CI") == "1":
+       refinements = 3
+
    est_errors = []
    sqrt_dofs = []
    mg_iterations = []
diff --git a/demos/benney_luke/benney_luke.py.rst b/demos/benney_luke/benney_luke.py.rst
index 0d035af52e..e08bc8fb46 100644
--- a/demos/benney_luke/benney_luke.py.rst
+++ b/demos/benney_luke/benney_luke.py.rst
@@ -100,6 +100,12 @@ Now we move on to defining parameters::
   mu = 0.01
   epsilon = 0.01
 
+  # Setup for faster test execution.
+  import os
+  if os.getenv("FIREDRAKE_CI") == "1":
+      T = 0.025
+      Nx = 20
+
   m = UnitIntervalMesh(Nx)
   mesh = ExtrudedMesh(m, layers=Ny)
   coords = mesh.coordinates
diff --git a/demos/fast_diagonalisation/fast_diagonalisation_poisson.py.rst b/demos/fast_diagonalisation/fast_diagonalisation_poisson.py.rst
index bc558ebc25..93fe1315e1 100644
--- a/demos/fast_diagonalisation/fast_diagonalisation_poisson.py.rst
+++ b/demos/fast_diagonalisation/fast_diagonalisation_poisson.py.rst
@@ -24,8 +24,14 @@ extruded hexahedral meshes, so we must create an :func:`~.ExtrudedMesh`. ::
 
   from firedrake import *
 
-  base = UnitSquareMesh(8, 8, quadrilateral=True)
-  mesh = ExtrudedMesh(base, 8)
+  # Setup for faster test execution.
+  import os
+  if os.getenv("FIREDRAKE_CI") == "1":
+      base = UnitSquareMesh(4, 4, quadrilateral=True)
+      mesh = ExtrudedMesh(base, 4)
+  else:
+      base = UnitSquareMesh(8, 8, quadrilateral=True)
+      mesh = ExtrudedMesh(base, 8)
 
 
 Defining the problem: the Poisson equation
@@ -88,7 +94,12 @@ using a sparse direct LU factorization. ::
   }
 
   fdm_lu_params = fdm_params(lu_params)
-  its = run_solve(5, fdm_lu_params)
+
+  if os.getenv("FIREDRAKE_CI") == "1":
+      its = run_solve(3, fdm_lu_params)
+  else:
+      its = run_solve(5, fdm_lu_params)
+
   print(f"LU iterations {its}")
 
 
@@ -123,6 +134,10 @@ implemented via :class:`~.ASMExtrudedStarPC` as we have an extruded mesh::
   print("FDM + ASM")
   print("Degree\tIterations")
   for degree in range(3, 6):
+
+      if os.getenv("FIREDRAKE_CI") == "1" and degree > 3:
+          break
+
       its = run_solve(degree, fdm_asm_params)
       print(f"{degree}\t{its}")
 
@@ -177,6 +192,10 @@ block, and the two-level additive Schwarz method on the facets. ::
   print('FDM + SC + ASM')
   print("Degree\tIterations")
   for degree in range(3, 6):
+
+      if os.getenv("FIREDRAKE_CI") == "1" and degree > 3:
+          break
+
       its = run_solve(degree, fdm_sc_asm_params)
       print(f"{degree}\t{its}")
 
diff --git a/demos/ma-demo/ma-demo.py.rst b/demos/ma-demo/ma-demo.py.rst
index 432042a20d..f9c0cb1d2e 100644
--- a/demos/ma-demo/ma-demo.py.rst
+++ b/demos/ma-demo/ma-demo.py.rst
@@ -76,6 +76,12 @@ mesh of quadrilaterals. ::
 
   from firedrake import *
   n = 100
+
+  # Setup for faster test execution.
+  import os
+  if os.getenv("FIREDRAKE_CI") == "1":
+      n = 20
+
   mesh = UnitSquareMesh(n, n, quadrilateral=True)
 
 We construct the quadratic function space for :math:`u`, ::
diff --git a/demos/matrix_free/rayleigh-benard.py.rst b/demos/matrix_free/rayleigh-benard.py.rst
index 8b0202d9fe..c33e83dd48 100644
--- a/demos/matrix_free/rayleigh-benard.py.rst
+++ b/demos/matrix_free/rayleigh-benard.py.rst
@@ -14,6 +14,11 @@ temperature. ::
 
   N = 128
 
+  # Setup for faster test execution.
+  import os
+  if os.getenv("FIREDRAKE_CI") == "1":
+      N = 16
+
   M = UnitSquareMesh(N, N)
 
   V = VectorFunctionSpace(M, "CG", 2)
diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst
index 6d0f488c83..cced2883b5 100755
--- a/demos/netgen/netgen_mesh.py.rst
+++ b/demos/netgen/netgen_mesh.py.rst
@@ -48,10 +48,19 @@ Now using the two of the predefined CSG geometries included in Netgen, a rectang
    geo.SetMaterial(2, "inner")
    geo.SetDomainMaxH(2, 0.02)
 
+   # Setup for faster test execution.
+   import os
+   if os.getenv("FIREDRAKE_CI") == "1":
+       geo.SetDomainMaxH(2, 0.1)
+
 Notice that the ``leftdomain`` and ``rightdomain`` attribute in the ``AddRectangle`` and ``AddCircle`` methods are used to set a domain index for the domain on the left and right side respectively of the rectangle and circle perimeter. It is worth mentioning that the perimeters are parametrised in a counterclockwise direction.
 We can now construct a mesh for the geometry we have defined and save it to a PVD file for visualisation. We will do so using the ``GenerateMesh`` method inside of the ``SplineGeometry`` class: ::
 
-   ngmsh = geo.GenerateMesh(maxh=0.1)
+   if os.getenv("FIREDRAKE_CI") == "1":
+       ngmsh = geo.GenerateMesh(maxh=0.3)
+   else:
+       ngmsh = geo.GenerateMesh(maxh=0.1)
+
    # Generating a Firedrake mesh from the NetGen mesh
    msh = Mesh(ngmsh)
    VTKFile("output/MeshExample1.pvd").write(msh)
@@ -141,6 +150,10 @@ We begin by defining some quantities of interest such as the desired tolerance,
 
    tolerance = 1e-16
    max_iterations = 10
+
+   if os.getenv("FIREDRAKE_CI") == "1":
+       max_iterations = 3
+
    exact = 3.375610652693620492628**2
 
 We create a function to solve the eigenvalue problem using SLEPc. We begin initialising the ``FunctionSpace``, the bilinear forms and linear functionals needed in the variational problem.
@@ -256,7 +269,12 @@ The `+,-,*` operators have respectively the meaning of union, set difference, an
    sphere = Sphere(Pnt(0.6, 0.6, 0.6), 0.5)
    geo = CSGeometry()
    geo.Add(cube-sphere)
-   ngmsh = geo.GenerateMesh(maxh=0.1)
+
+   if os.getenv("FIREDRAKE_CI") == "1":
+       ngmsh = geo.GenerateMesh(maxh=0.3)
+   else:
+       ngmsh = geo.GenerateMesh(maxh=0.1)
+
    msh = Mesh(ngmsh)
    VTKFile("output/MeshExample3.pvd").write(msh)
 
@@ -399,6 +417,10 @@ It is also possible to construct high-order meshes using the ``SplineGeometry``,
    if COMM_WORLD.rank == 0:
        geo = CSG2d()
        circle = Circle(center=(1,1), radius=0.1, bc="curve").Maxh(0.01)
+
+       if os.getenv("FIREDRAKE_CI") == "1":
+           circle = Circle(center=(1,1), radius=0.1, bc="curve").Maxh(0.05)
+
        rect = Rectangle(pmin=(0,1), pmax=(1,2),
                         bottom="b", left="l", top="t", right="r")
        geo.Add(rect-circle)
diff --git a/demos/poisson_adaptive_eigenvalue/poisson_adaptive_eigenvalue.py.rst b/demos/poisson_adaptive_eigenvalue/poisson_adaptive_eigenvalue.py.rst
index 884af071e1..b128aa32fa 100644
--- a/demos/poisson_adaptive_eigenvalue/poisson_adaptive_eigenvalue.py.rst
+++ b/demos/poisson_adaptive_eigenvalue/poisson_adaptive_eigenvalue.py.rst
@@ -130,6 +130,12 @@ We define a function to adapt the mesh by refining elements with large error ind
 Finally, we run the adaptive loop until the upper and lower bounds agree to within a tolerance. ::
 
   max_iterations = 20
+
+  # Setup for faster test execution.
+  import os
+  if os.getenv("FIREDRAKE_CI") == "1":
+      max_iterations = 2
+
   error_estimators = []
   dofs = []
   err = 1
diff --git a/demos/reynolds_robust_navier_stokes_hdiv/reynolds_robust_navier_stokes_hdiv.py.rst b/demos/reynolds_robust_navier_stokes_hdiv/reynolds_robust_navier_stokes_hdiv.py.rst
index d9bf47ad57..9033f2c594 100644
--- a/demos/reynolds_robust_navier_stokes_hdiv/reynolds_robust_navier_stokes_hdiv.py.rst
+++ b/demos/reynolds_robust_navier_stokes_hdiv/reynolds_robust_navier_stokes_hdiv.py.rst
@@ -199,9 +199,16 @@ which is necessary for vertex-star relaxation in parallel. ::
   print = PETSc.Sys.Print
 
   distribution_parameters = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)}
-  num_refinements = 2
-  base = UnitSquareMesh(16, 16, distribution_parameters=distribution_parameters)
-  mh = MeshHierarchy(base, num_refinements)
+
+  # Setup for faster test execution.
+  import os
+  if os.getenv("FIREDRAKE_CI") == "1":
+      base = UnitSquareMesh(8, 8, distribution_parameters=distribution_parameters)
+      mh = MeshHierarchy(base, 1)
+  else:
+      base = UnitSquareMesh(16, 16, distribution_parameters=distribution_parameters)
+      mh = MeshHierarchy(base, 2)
+
   mesh = mh[-1]
   n = FacetNormal(mesh)
   (x, y) = SpatialCoordinate(mesh)
@@ -396,7 +403,11 @@ precision. ::
   problem = NonlinearVariationalProblem(F, w, bcs, Jp=Jp)
   solver = NonlinearVariationalSolver(problem, solver_parameters=sp, pre_apply_bcs=False)
 
-  for Re_ in [1, 500]  + list(range(1000, 5001, 1000)):
+  for Re_ in [1, 500] + list(range(1000, 5001, 1000)):
+
+      if os.getenv("FIREDRAKE_CI") == "1" and Re_ > 500:
+          break
+
       Re.assign(Re_)
 
       # Solve
diff --git a/demos/saddle_point_pc/saddle_point_systems.py.rst b/demos/saddle_point_pc/saddle_point_systems.py.rst
index cc9117a9fd..0978b08bf0 100644
--- a/demos/saddle_point_pc/saddle_point_systems.py.rst
+++ b/demos/saddle_point_pc/saddle_point_systems.py.rst
@@ -173,6 +173,12 @@ solving it ::
 
     print("Naive preconditioning")
     for n in range(8):
+
+        # Setup for faster test execution.
+        import os
+        if os.getenv("FIREDRAKE_CI") == "1" and n == 3:
+            break
+
         solver, w = build_problem(n, parameters, block_matrix=False)
         solver.solve()
 
@@ -280,6 +286,10 @@ applying the action of blocks, so we can use a block matrix format. ::
 
     print("Exact full Schur complement")
     for n in range(8):
+
+        if os.getenv("FIREDRAKE_CI") == "1" and n == 3:
+            break
+
         solver, w = build_problem(n, parameters, block_matrix=True)
         solver.solve()
         print(w.function_space().mesh().unique().num_cells(), solver.snes.ksp.getIterationNumber())
@@ -365,6 +375,10 @@ Let's see what happens. ::
 
     print("Schur complement with S_p")
     for n in range(8):
+
+        if os.getenv("FIREDRAKE_CI") == "1" and n == 3:
+            break
+
         solver, w = build_problem(n, parameters, block_matrix=True)
         solver.solve()
         print(w.function_space().mesh().unique().num_cells(), solver.snes.ksp.getIterationNumber())
@@ -420,6 +434,10 @@ and so we no longer need a flexible Krylov method. ::
 
     print("Schur complement with S_p and inexact inner inverses")
     for n in range(8):
+
+        if os.getenv("FIREDRAKE_CI") == "1" and n == 3:
+            break
+
         solver, w = build_problem(n, parameters, block_matrix=True)
         solver.solve()
         print(w.function_space().mesh().unique().num_cells(), solver.snes.ksp.getIterationNumber())
@@ -485,6 +503,10 @@ variable. We can provide it as an :class:`~.AuxiliaryOperatorPC` via a python pr
 
     print("DG approximation for S_p")
     for n in range(8):
+
+        if os.getenv("FIREDRAKE_CI") == "1" and n == 3:
+            break
+
         solver, w = build_problem(n, parameters, aP=None, block_matrix=False)
         solver.solve()
         print(w.function_space().mesh().unique().num_cells(), solver.snes.ksp.getIterationNumber())
@@ -569,6 +591,10 @@ Let's see what the iteration count looks like now. ::
 
     print("Riesz-map preconditioner")
     for n in range(8):
+
+        if os.getenv("FIREDRAKE_CI") == "1" and n == 3:
+            break
+
         solver, w = build_problem(n, parameters, aP=riesz, block_matrix=True)
         solver.solve()
         print(w.function_space().mesh().unique().num_cells(), solver.snes.ksp.getIterationNumber())
diff --git a/docs/notebooks/07-geometric-multigrid.py b/docs/notebooks/07-geometric-multigrid.py
index 4f6fd0171e..d870f129d7 100644
--- a/docs/notebooks/07-geometric-multigrid.py
+++ b/docs/notebooks/07-geometric-multigrid.py
@@ -29,7 +29,12 @@
 # %%
 from firedrake import *
 
-coarse_mesh = RectangleMesh(15, 10, 1.5, 1)
+# Setup for faster test execution.
+import os
+if os.getenv("FIREDRAKE_CI") == "1":
+    coarse_mesh = RectangleMesh(5, 4, 1.5, 1)
+else:
+    coarse_mesh = RectangleMesh(15, 10, 1.5, 1)
 
 # %% [markdown]
 # Having made the coarse mesh, we create the hierarchy of meshes.  The second argument tells Firedrake how many levels of refinement to use.  Here we refine three times, so that in total we have four meshes.
@@ -116,7 +121,10 @@
 
 # %%
 def create_solver(parameters=None):
-    coarse_mesh = RectangleMesh(15, 10, 1.5, 1)
+    if os.getenv("FIREDRAKE_CI") == "1":
+        coarse_mesh = RectangleMesh(5, 4, 1.5, 1)
+    else:
+        coarse_mesh = RectangleMesh(15, 10, 1.5, 1)
     hierarchy = MeshHierarchy(coarse_mesh, 3)
     
     mesh = hierarchy[-1]
@@ -277,7 +285,10 @@ def form(self, pc, test, trial):
 
 # %%
 def create_solver(parameters=None):
-    coarse_mesh = RectangleMesh(15, 10, 1.5, 1)
+    if os.getenv("FIREDRAKE_CI") == "1":
+        coarse_mesh = RectangleMesh(5, 4, 1.5, 1)
+    else:
+        coarse_mesh = RectangleMesh(15, 10, 1.5, 1)
     hierarchy = MeshHierarchy(coarse_mesh, 3)
     
     mesh = hierarchy[-1]
diff --git a/docs/notebooks/09-hybridisation.py b/docs/notebooks/09-hybridisation.py
index 4cb120a7bd..a04e8e5f27 100644
--- a/docs/notebooks/09-hybridisation.py
+++ b/docs/notebooks/09-hybridisation.py
@@ -93,7 +93,12 @@
 # So now we need to specify the number of refinements (say 4). Let's also set the mesh `degree` to be cubic. With this choice of coordinate space, we can better resolve the actual curvature of the sphere using bendy quadrilateral elements:
 
 # %%
-mesh = CubedSphereMesh(radius=R0, refinement_level=4, degree=3)
+# Use a coarser mesh when running CI tests.
+import os
+if os.getenv("FIREDRAKE_CI") == "1":
+    mesh = CubedSphereMesh(radius=R0, refinement_level=2, degree=3)
+else:
+    mesh = CubedSphereMesh(radius=R0, refinement_level=4, degree=3)
 
 # %% [markdown]
 # And now we just initialize the global normals on this mesh:
diff --git a/docs/notebooks/10-sum-factorisation.py b/docs/notebooks/10-sum-factorisation.py
index c5f2e7a131..266aca1c49 100644
--- a/docs/notebooks/10-sum-factorisation.py
+++ b/docs/notebooks/10-sum-factorisation.py
@@ -64,7 +64,12 @@
 # We can create a hexahedral mesh by extruding a quadrilateral mesh.
 
 # %%
-mesh = ExtrudedMesh(UnitSquareMesh(10, 10, quadrilateral=True), 10)
+# Setup for faster test execution.
+import os
+if os.getenv("FIREDRAKE_CI") == "1":
+    mesh = ExtrudedMesh(UnitSquareMesh(2, 2, quadrilateral=True), 2)
+else:
+    mesh = ExtrudedMesh(UnitSquareMesh(10, 10, quadrilateral=True), 10)
 
 # %% [markdown]
 # Let's choose the continuous Lagrange element of degree 5 as our function space.
@@ -160,6 +165,10 @@ def gauss_lobatto_legendre_cube_rule(dimension, degree):
 # %%
 flops = defaultdict(list)
 ps = range(1, 33)  # polynomial degrees
+
+if os.getenv("FIREDRAKE_CI") == "1":
+    ps = [1, 2, 4, 8]
+
 modes = {
     'gll': {'mode': 'spectral', 'variant': 'spectral', 'rule': gauss_lobatto_legendre_cube_rule},
     'spectral': {'mode': 'spectral', 'variant': None, 'rule': lambda *args: None},
@@ -194,6 +203,10 @@ def gauss_lobatto_legendre_cube_rule(dimension, degree):
 # This might take some time to run
 flops_curl = defaultdict(list)
 ps_curl = range(1, 17)
+
+if os.getenv("FIREDRAKE_CI") == "1":
+    ps_curl = [1, 2, 4]
+
 for p in ps_curl:
     for mode in modes:
         element = FiniteElement('NCE', mesh.ufl_cell(), degree=p, variant=modes[mode]['variant'])
diff --git a/docs/notebooks/12-HPC_demo.py b/docs/notebooks/12-HPC_demo.py
index aabebc59c4..8abb38ee1a 100644
--- a/docs/notebooks/12-HPC_demo.py
+++ b/docs/notebooks/12-HPC_demo.py
@@ -35,6 +35,12 @@
 Nref = 2
 degree = 2
 
+# Setup for faster test execution.
+import os
+if os.getenv("FIREDRAKE_CI") == "1":
+    Nx = 2
+    Nref = 1
+
 # %% [markdown]
 # These three parameters determine the total number of degrees of freedom (DOFs) in our problem:
 # - `Nx` defines our coarse grid in the mesh hierarchy, it is used to construct a coarse cube mesh.