diff --git a/goalie/error_estimation.py b/goalie/error_estimation.py
index d7a18b1..76155d4 100644
--- a/goalie/error_estimation.py
+++ b/goalie/error_estimation.py
@@ -98,7 +98,7 @@ def get_dwr_indicator(F, adjoint_error, test_space=None):
     elif not isinstance(adjoint_error, dict):
         raise TypeError(
             "Expected 'adjoint_error' to be a Function or dict, not"
-            " '{type(adjoint_error)}'."
+            f" '{type(adjoint_error)}'."
         )
 
     # Process input for test_space as a dictionary
@@ -111,7 +111,7 @@ def get_dwr_indicator(F, adjoint_error, test_space=None):
     elif not isinstance(test_space, dict):
         raise TypeError(
             "Expected 'test_space' to be a FunctionSpace or dict, not"
-            " '{type(test_space)}'."
+            f" '{type(test_space)}'."
         )
 
     # Construct the mapping for each component
@@ -122,7 +122,7 @@ def get_dwr_indicator(F, adjoint_error, test_space=None):
         if not isinstance(fs, WithGeometry):
             raise TypeError(
                 f"Expected 'test_space['{key}']' to be a FunctionSpace, not"
-                " '{type(fs)}'."
+                f" '{type(fs)}'."
             )
         if F.ufl_domain() != err.function_space().mesh():
             raise ValueError(
diff --git a/pyproject.toml b/pyproject.toml
index 3b03d8d..1a4fd52 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -52,7 +52,6 @@ select = [
 ]
 ignore = [
   "C901",  # function is too complex (TODO #165: enable this)
-  "F403",  # unable to detect undefined names
   "F405",  # name may be undefined, or defined from star imports
 ]
 [tool.ruff.lint.per-file-ignores]
@@ -61,6 +60,7 @@ ignore = [
 ]
 "demos/*.py" = [
   "E402",  # module level import not at top of file
+  "F403",  # unable to detect undefined names
 ]
 
 [tool.pytest.ini_options]
diff --git a/test/adjoint/examples/burgers.py b/test/adjoint/examples/burgers.py
index 6e1042e..390e9f0 100644
--- a/test/adjoint/examples/burgers.py
+++ b/test/adjoint/examples/burgers.py
@@ -7,8 +7,14 @@
     https://firedrakeproject.org/demos/burgers.py.html
 """
 
-from firedrake import *
+import ufl
 from firedrake.__future__ import interpolate
+from firedrake.assemble import assemble
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace, VectorFunctionSpace
+from firedrake.solving import solve
+from firedrake.ufl_expr import TestFunction
+from firedrake.utility_meshes import UnitSquareMesh
 
 # Problem setup
 n = 32
@@ -22,9 +28,7 @@
 
 
 def get_function_spaces(mesh):
-    r"""
-    :math:`\mathbb P2` space.
-    """
+    r""":math:`\mathbb P2` space."""
     return {"uv_2d": VectorFunctionSpace(mesh, "CG", 2)}
 
 
@@ -47,9 +51,9 @@ def solver(i):
         nu = Function(R).assign(0.0001)
         v = TestFunction(fs)
         F = (
-            inner((u - u_) / dtc, v) * dx
-            + inner(dot(u, nabla_grad(u)), v) * dx
-            + nu * inner(grad(u), grad(v)) * dx
+            ufl.inner((u - u_) / dtc, v) * ufl.dx
+            + ufl.inner(ufl.dot(u, ufl.nabla_grad(u)), v) * ufl.dx
+            + nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
         )
 
         # Time integrate from t_start to t_end
@@ -72,8 +76,10 @@ def get_initial_condition(self):
     Initial condition which is sinusoidal in the x-direction.
     """
     init_fs = self.function_spaces["uv_2d"][0]
-    x, y = SpatialCoordinate(self.meshes[0])
-    return {"uv_2d": assemble(interpolate(as_vector([sin(pi * x), 0]), init_fs))}
+    x, y = ufl.SpatialCoordinate(self.meshes[0])
+    return {
+        "uv_2d": assemble(interpolate(ufl.as_vector([ufl.sin(ufl.pi * x), 0]), init_fs))
+    }
 
 
 def get_qoi(self, i):
@@ -86,7 +92,7 @@ def get_qoi(self, i):
 
     def time_integrated_qoi(t):
         u = self.fields["uv_2d"][0]
-        return dtc * inner(u, u) * ds(2)
+        return dtc * ufl.inner(u, u) * ufl.ds(2)
 
     def end_time_qoi():
         return time_integrated_qoi(end_time)
diff --git a/test/adjoint/examples/point_discharge2d.py b/test/adjoint/examples/point_discharge2d.py
index e0f18ce..05af42a 100644
--- a/test/adjoint/examples/point_discharge2d.py
+++ b/test/adjoint/examples/point_discharge2d.py
@@ -1,23 +1,24 @@
 """
-Problem specification for a simple
-advection-diffusion test case with a
-point source, from [Riadh et al. 2014].
-
-This test case is notable for Goalie
-because it has an analytical solution,
-meaning the effectivity index can be
-computed.
-
-[Riadh et al. 2014] A. Riadh, G.
-    Cedric, M. Jean, "TELEMAC modeling
-    system: 2D hydrodynamics TELEMAC-2D
-    software release 7.0 user manual."
-    Paris: R&D, Electricite de France,
-    p. 134 (2014).
+Problem specification for a simple advection-diffusion test case with a point source,
+from [Riadh et al. 2014].
+
+This test case is notable for Goalie because it has an analytical solution, meaning the
+effectivity index can be computed.
+
+[Riadh et al. 2014] A. Riadh, G. Cedric, M. Jean, "TELEMAC modeling system: 2D
+    hydrodynamics TELEMAC-2D software release 7.0 user manual." Paris: R&D, Electricite
+    de France, p. 134 (2014).
 """
 
 import numpy as np
-from firedrake import *
+import ufl
+from firedrake.assemble import assemble
+from firedrake.bcs import DirichletBC
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace
+from firedrake.solving import solve
+from firedrake.ufl_expr import TestFunction
+from firedrake.utility_meshes import RectangleMesh
 
 from goalie.math import bessk0
 
@@ -43,19 +44,15 @@ def get_function_spaces(mesh):
 
 def source(mesh):
     """
-    Gaussian approximation to a point source
-    at (2, 5) with discharge rate 100 on a
+    Gaussian approximation to a point source at (2, 5) with discharge rate 100 on a
     given mesh.
     """
-    x, y = SpatialCoordinate(mesh)
-    return 100.0 * exp(-((x - src_x) ** 2 + (y - src_y) ** 2) / src_r**2)
+    x, y = ufl.SpatialCoordinate(mesh)
+    return 100.0 * ufl.exp(-((x - src_x) ** 2 + (y - src_y) ** 2) / src_r**2)
 
 
 def get_solver(self):
-    """
-    Advection-diffusion equation
-    solved using a direct method.
-    """
+    """Advection-diffusion equation solved using a direct method."""
 
     def solver(i):
         fs = self.function_spaces["tracer_2d"][i]
@@ -67,22 +64,22 @@ def solver(i):
         D = Function(R).assign(0.1)
         u_x = Function(R).assign(1.0)
         u_y = Function(R).assign(0.0)
-        u = as_vector([u_x, u_y])
-        h = CellSize(self[i])
+        u = ufl.as_vector([u_x, u_y])
+        h = ufl.CellSize(self[i])
         S = source(self[i])
 
         # Stabilisation parameter
-        unorm = sqrt(dot(u, u))
+        unorm = ufl.sqrt(ufl.dot(u, u))
         tau = 0.5 * h / unorm
-        tau = min_value(tau, unorm * h / (6 * D))
+        tau = ufl.min_value(tau, unorm * h / (6 * D))
 
         # Setup variational problem
         psi = TestFunction(fs)
-        psi = psi + tau * dot(u, grad(psi))
+        psi = psi + tau * ufl.dot(u, ufl.grad(psi))
         F = (
-            S * psi * dx
-            - dot(u, grad(c)) * psi * dx
-            - inner(D * grad(c), grad(psi)) * dx
+            S * psi * ufl.dx
+            - ufl.dot(u, ufl.grad(c)) * psi * ufl.dx
+            - ufl.inner(D * ufl.grad(c), ufl.grad(psi)) * ufl.dx
         )
 
         # Zero Dirichlet condition on the left-hand (inlet) boundary
@@ -104,32 +101,30 @@ def solver(i):
 
 def get_qoi(self, i):
     """
-    Quantity of interest which integrates
-    the tracer concentration over an offset
+    Quantity of interest which integrates the tracer concentration over an offset
     receiver region.
     """
 
     def steady_qoi():
         c = self.fields["tracer_2d"]
-        x, y = SpatialCoordinate(self[i])
-        kernel = conditional((x - rec_x) ** 2 + (y - rec_y) ** 2 < rec_r**2, 1, 0)
-        area = assemble(kernel * dx)
-        area_analytical = pi * rec_r**2
+        x, y = ufl.SpatialCoordinate(self[i])
+        kernel = ufl.conditional((x - rec_x) ** 2 + (y - rec_y) ** 2 < rec_r**2, 1, 0)
+        area = assemble(kernel * ufl.dx)
+        area_analytical = ufl.pi * rec_r**2
         scaling = 1.0 if np.allclose(area, 0.0) else area_analytical / area
-        return scaling * kernel * c * dx
+        return scaling * kernel * c * ufl.dx
 
     return steady_qoi
 
 
 def analytical_solution(mesh):
     """
-    Analytical solution as represented on
-    a given mesh. See [Riadh et al. 2014].
+    Analytical solution as represented on a given mesh. See [Riadh et al. 2014].
     """
-    x, y = SpatialCoordinate(mesh)
+    x, y = ufl.SpatialCoordinate(mesh)
     R = FunctionSpace(mesh, "R", 0)
     u = Function(R).assign(1.0)
     D = Function(R).assign(0.1)
     Pe = 0.5 * u / D
-    r = max_value(sqrt((x - src_x) ** 2 + (y - src_y) ** 2), src_r)
-    return 0.5 / (pi * D) * exp(Pe * (x - src_x)) * bessk0(Pe * r)
+    r = ufl.max_value(ufl.sqrt((x - src_x) ** 2 + (y - src_y) ** 2), src_r)
+    return 0.5 / (ufl.pi * D) * ufl.exp(Pe * (x - src_x)) * bessk0(Pe * r)
diff --git a/test/adjoint/examples/point_discharge3d.py b/test/adjoint/examples/point_discharge3d.py
index 29a4d36..2fb0779 100644
--- a/test/adjoint/examples/point_discharge3d.py
+++ b/test/adjoint/examples/point_discharge3d.py
@@ -1,32 +1,28 @@
 """
-Problem specification for a simple
-advection-diffusion test case with a
-point source. Extended from
-[Riadh et al. 2014] as in
-[Wallwork et al. 2020].
-
-This test case is notable for Goalie
-because it is in 3D and has an
-analytical solution, meaning the
-effectivity index can be computed.
-
-[Riadh et al. 2014] A. Riadh, G.
-    Cedric, M. Jean, "TELEMAC modeling
-    system: 2D hydrodynamics TELEMAC-2D
-    software release 7.0 user manual."
-    Paris: R&D, Electricite de France,
-    p. 134 (2014).
-
-[Wallwork et al. 2020] J.G. Wallwork,
-    N. Barral, D.A. Ham, M.D. Piggott,
-    "Anisotropic Goal-Oriented Mesh
-    Adaptation in Firedrake". In:
-    Proceedings of the 28th International
-    Meshing Roundtable (2020).
+Problem specification for a simple advection-diffusion test case with a point source.
+Extended from [Riadh et al. 2014] as in [Wallwork et al. 2020].
+
+This test case is notable for Goalie because it is in 3D and has an analytical solution,
+meaning the effectivity index can be computed.
+
+[Riadh et al. 2014] A. Riadh, G. Cedric, M. Jean, "TELEMAC modeling system: 2D
+    hydrodynamics TELEMAC-2D software release 7.0 user manual." Paris: R&D, Electricite
+    de France, p. 134 (2014).
+
+[Wallwork et al. 2020] J.G. Wallwork, N. Barral, D.A. Ham, M.D. Piggott, "Anisotropic
+    Goal-Oriented Mesh Adaptation in Firedrake". In: Proceedings of the 28th
+    International Meshing Roundtable (2020).
 """
 
 import numpy as np
-from firedrake import *
+import ufl
+from firedrake.assemble import assemble
+from firedrake.bcs import DirichletBC
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace
+from firedrake.solving import solve
+from firedrake.ufl_expr import TestFunction
+from firedrake.utility_meshes import BoxMesh
 
 from goalie.math import bessk0
 
@@ -44,29 +40,23 @@
 
 
 def get_function_spaces(mesh):
-    r"""
-    :math:`\mathbb P1` space.
-    """
+    r""":math:`\mathbb P1` space."""
     return {"tracer_3d": FunctionSpace(mesh, "CG", 1)}
 
 
 def source(mesh):
     """
-    Gaussian approximation to a point source
-    at (2, 5, 5) with discharge rate 100 on a
+    Gaussian approximation to a point source at (2, 5, 5) with discharge rate 100 on a
     given mesh.
     """
-    x, y, z = SpatialCoordinate(mesh)
-    return 100.0 * exp(
+    x, y, z = ufl.SpatialCoordinate(mesh)
+    return 100.0 * ufl.exp(
         -((x - src_x) ** 2 + (y - src_y) ** 2 + (z - src_z) ** 2) / src_r**2
     )
 
 
 def get_solver(self):
-    """
-    Advection-diffusion equation
-    solved using a direct method.
-    """
+    """Advection-diffusion equation solved using a direct method."""
 
     def solver(i):
         fs = self.function_spaces["tracer_3d"][i]
@@ -79,22 +69,22 @@ def solver(i):
         u_x = Function(R).assign(1.0)
         u_y = Function(R).assign(0.0)
         u_z = Function(R).assign(0.0)
-        u = as_vector([u_x, u_y, u_z])
-        h = CellSize(self[i])
+        u = ufl.as_vector([u_x, u_y, u_z])
+        h = ufl.CellSize(self[i])
         S = source(self[i])
 
         # Stabilisation parameter
-        unorm = sqrt(dot(u, u))
+        unorm = ufl.sqrt(ufl.dot(u, u))
         tau = 0.5 * h / unorm
-        tau = min_value(tau, unorm * h / (6 * D))
+        tau = ufl.min_value(tau, unorm * h / (6 * D))
 
         # Setup variational problem
         psi = TestFunction(fs)
-        psi = psi + tau * dot(u, grad(psi))
+        psi = psi + tau * ufl.dot(u, ufl.grad(psi))
         F = (
-            S * psi * dx
-            - dot(u, grad(c)) * psi * dx
-            - inner(D * grad(c), grad(psi)) * dx
+            S * psi * ufl.dx
+            - ufl.dot(u, ufl.grad(c)) * psi * ufl.dx
+            - ufl.inner(D * ufl.grad(c), ufl.grad(psi)) * ufl.dx
         )
 
         # Zero Dirichlet condition on the left-hand (inlet) boundary
@@ -116,34 +106,32 @@ def solver(i):
 
 def get_qoi(self, i):
     """
-    Quantity of interest which integrates
-    the tracer concentration over an offset
+    Quantity of interest which integrates the tracer concentration over an offset
     receiver region.
     """
 
     def steady_qoi():
         c = self.fields["tracer_3d"]
-        x, y, z = SpatialCoordinate(self[i])
-        kernel = conditional(
+        x, y, z = ufl.SpatialCoordinate(self[i])
+        kernel = ufl.conditional(
             (x - rec_x) ** 2 + (y - rec_y) ** 2 + (z - rec_z) ** 2 < rec_r**2, 1, 0
         )
-        area = assemble(kernel * dx)
-        area_analytical = pi * rec_r**2
+        area = assemble(kernel * ufl.dx)
+        area_analytical = ufl.pi * rec_r**2
         scaling = 1.0 if np.allclose(area, 0.0) else area_analytical / area
-        return scaling * kernel * c * dx
+        return scaling * kernel * c * ufl.dx
 
     return steady_qoi
 
 
 def analytical_solution(mesh):
-    """
-    Analytical solution as represented on
-    a given mesh.
-    """
-    x, y, z = SpatialCoordinate(mesh)
+    """Analytical solution as represented on a given mesh."""
+    x, y, z = ufl.SpatialCoordinate(mesh)
     R = FunctionSpace(mesh, "R", 0)
     u = Function(R).assign(1.0)
     D = Function(R).assign(0.1)
     Pe = 0.5 * u / D
-    r = max_value(sqrt((x - src_x) ** 2 + (y - src_y) ** 2 + (z - src_z) ** 2), src_r)
-    return 0.5 / (pi * D) * exp(Pe * (x - src_x)) * bessk0(Pe * r)
+    r = ufl.max_value(
+        ufl.sqrt((x - src_x) ** 2 + (y - src_y) ** 2 + (z - src_z) ** 2), src_r
+    )
+    return 0.5 / (ufl.pi * D) * ufl.exp(Pe * (x - src_x)) * bessk0(Pe * r)
diff --git a/test/adjoint/examples/steady_flow_past_cyl.py b/test/adjoint/examples/steady_flow_past_cyl.py
index 811456e..7a4bbf2 100644
--- a/test/adjoint/examples/steady_flow_past_cyl.py
+++ b/test/adjoint/examples/steady_flow_past_cyl.py
@@ -1,11 +1,8 @@
 """
-Problem specification for a simple steady
-state flow-past-a-cylinder test case which
+Problem specification for a simple steady state flow-past-a-cylinder test case which
 solves a Stokes problem.
 
-The test case is notable for Goalie
-because the prognostic equation is
-nonlinear.
+The test case is notable for Goalie because the prognostic equation is nonlinear.
 
 Code here is based on that found at
     https://nbviewer.jupyter.org/github/firedrakeproject/firedrake/blob/master/docs/notebooks/06-pde-constrained-optimisation.ipynb
@@ -13,8 +10,15 @@
 
 import os
 
-from firedrake import *
+import ufl
 from firedrake.__future__ import interpolate
+from firedrake.assemble import assemble
+from firedrake.bcs import DirichletBC
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace, VectorFunctionSpace
+from firedrake.mesh import Mesh
+from firedrake.solving import solve
+from firedrake.ufl_expr import TestFunctions
 
 mesh = Mesh(os.path.join(os.path.dirname(__file__), "mesh-with-hole.msh"))
 fields = ["up"]
@@ -26,17 +30,12 @@
 
 
 def get_function_spaces(mesh):
-    r"""
-    Taylor-Hood :math:`\mathbb P2-\mathbb P1` space.
-    """
+    r"""Taylor-Hood :math:`\mathbb P2-\mathbb P1` space."""
     return {"up": VectorFunctionSpace(mesh, "CG", 2) * FunctionSpace(mesh, "CG", 1)}
 
 
 def get_solver(self):
-    """
-    Stokes problem solved using a
-    direct method.
-    """
+    """Stokes problem solved using a direct method."""
 
     def solver(i):
         W = self.function_spaces["up"][i]
@@ -45,18 +44,18 @@ def solver(i):
         # Define variational problem
         R = FunctionSpace(self[i], "R", 0)
         nu = Function(R).assign(1.0)
-        u, p = split(up)
+        u, p = ufl.split(up)
         v, q = TestFunctions(W)
         F = (
-            inner(dot(u, nabla_grad(u)), v) * dx
-            + nu * inner(grad(u), grad(v)) * dx
-            - inner(p, div(v)) * dx
-            - inner(q, div(u)) * dx
+            ufl.inner(ufl.dot(u, ufl.nabla_grad(u)), v) * ufl.dx
+            + nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+            - ufl.inner(p, ufl.div(v)) * ufl.dx
+            - ufl.inner(q, ufl.div(u)) * ufl.dx
         )
 
         # Define inflow and no-slip boundary conditions
-        y = SpatialCoordinate(self[i])[1]
-        u_inflow = as_vector([y * (10 - y) / 25.0, 0])
+        y = ufl.SpatialCoordinate(self[i])[1]
+        u_inflow = ufl.as_vector([y * (10 - y) / 25.0, 0])
         noslip = DirichletBC(W.sub(0), (0, 0), (3, 5))
         inflow = DirichletBC(W.sub(0), assemble(interpolate(u_inflow, W.sub(0))), 1)
         bcs = [inflow, noslip, DirichletBC(W.sub(0), 0, 4)]
@@ -84,12 +83,11 @@ def solver(i):
 
 def get_initial_condition(self):
     """
-    Dummy initial condition function which
-    acts merely to pass over the
+    Dummy initial condition function which acts merely to pass over the
     :class:`FunctionSpace`.
     """
-    x, y = SpatialCoordinate(self[0])
-    u_inflow = as_vector([y * (10 - y) / 25.0, 0])
+    x, y = ufl.SpatialCoordinate(self[0])
+    u_inflow = ufl.as_vector([y * (10 - y) / 25.0, 0])
     up = Function(self.function_spaces["up"][0])
     u, p = up.subfunctions
     u.interpolate(u_inflow)
@@ -97,13 +95,10 @@ def get_initial_condition(self):
 
 
 def get_qoi(self, i):
-    """
-    Quantity of interest which integrates
-    pressure over the boundary of the hole.
-    """
+    """Quantity of interest which integrates pressure over the boundary of the hole."""
 
     def steady_qoi():
-        u, p = split(self.fields["up"])
-        return p * ds(4)
+        u, p = ufl.split(self.fields["up"])
+        return p * ufl.ds(4)
 
     return steady_qoi
diff --git a/test/adjoint/test_adjoint.py b/test/adjoint/test_adjoint.py
index b917260..40377eb 100644
--- a/test/adjoint/test_adjoint.py
+++ b/test/adjoint/test_adjoint.py
@@ -7,12 +7,17 @@
 import sys
 import unittest
 
+import numpy as np
 import pyadjoint
 import pytest
 from animate.utility import errornorm, norm
-from firedrake import Cofunction, UnitTriangleMesh
+from firedrake.cofunction import Cofunction
+from firedrake.utility_meshes import UnitTriangleMesh
 
-from goalie_adjoint import *
+from goalie.adjoint import AdjointMeshSeq
+from goalie.log import DEBUG, pyrint, set_log_level
+from goalie.time_partition import TimeInterval, TimePartition
+from goalie.utility import AttrDict, VTKFile
 
 sys.path.append(os.path.join(os.path.dirname(__file__), "examples"))
 
diff --git a/test/adjoint/test_adjoint_mesh_seq.py b/test/adjoint/test_adjoint_mesh_seq.py
index 9ad2365..c5e230a 100644
--- a/test/adjoint/test_adjoint_mesh_seq.py
+++ b/test/adjoint/test_adjoint_mesh_seq.py
@@ -24,10 +24,10 @@
 from parameterized import parameterized
 from pyadjoint.block_variable import BlockVariable
 
+from goalie.adjoint import AdjointMeshSeq
 from goalie.go_mesh_seq import GoalOrientedMeshSeq
-from goalie.log import *
-from goalie.time_partition import TimeInterval
-from goalie_adjoint import *
+from goalie.log import WARNING
+from goalie.time_partition import TimeInterval, TimePartition
 
 
 class TestBlockLogic(unittest.TestCase):
diff --git a/test/adjoint/test_demos.py b/test/adjoint/test_demos.py
index 9d03686..333c1ec 100644
--- a/test/adjoint/test_demos.py
+++ b/test/adjoint/test_demos.py
@@ -10,7 +10,7 @@
 
 import pytest
 
-from goalie.log import *
+from goalie.log import WARNING, set_log_level
 
 cwd = os.path.abspath(os.path.dirname(__file__))
 demo_dir = os.path.abspath(os.path.join(cwd, "..", "..", "demos"))
diff --git a/test/adjoint/test_fp_iteration.py b/test/adjoint/test_fp_iteration.py
index 08d8763..93a873f 100644
--- a/test/adjoint/test_fp_iteration.py
+++ b/test/adjoint/test_fp_iteration.py
@@ -2,17 +2,25 @@
 import unittest
 from unittest.mock import MagicMock
 
-from firedrake import *
+import numpy as np
+import ufl
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace
+from firedrake.utility_meshes import UnitSquareMesh, UnitTriangleMesh
 from parameterized import parameterized
 
-from goalie_adjoint import *
+from goalie.adjoint import AdjointMeshSeq
+from goalie.go_mesh_seq import GoalOrientedAdaptParameters, GoalOrientedMeshSeq
+from goalie.mesh_seq import MeshSeq
+from goalie.options import AdaptParameters
+from goalie.time_partition import TimeInstant, TimePartition
 
 
 def constant_qoi(mesh_seq, index):
     R = FunctionSpace(mesh_seq[index], "R", 0)
 
     def qoi():
-        return Function(R).assign(1) * dx
+        return Function(R).assign(1) * ufl.dx
 
     return qoi
 
@@ -21,7 +29,7 @@ def oscillating_qoi(mesh_seq, index):
     R = FunctionSpace(mesh_seq[index], "R", 0)
 
     def qoi():
-        return Function(R).assign(1 if mesh_seq.fp_iteration % 2 == 0 else 2) * dx
+        return Function(R).assign(1 if mesh_seq.fp_iteration % 2 == 0 else 2) * ufl.dx
 
     return qoi
 
diff --git a/test/adjoint/test_utils.py b/test/adjoint/test_utils.py
index 91d54ad..ee3b3ba 100644
--- a/test/adjoint/test_utils.py
+++ b/test/adjoint/test_utils.py
@@ -1,10 +1,14 @@
 import unittest
 
 import numpy as np
-from firedrake import *
+import ufl
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace
+from firedrake.utility_meshes import UnitSquareMesh
 
-from goalie.adjoint import annotate_qoi
-from goalie_adjoint import *
+from goalie.adjoint import AdjointMeshSeq, annotate_qoi
+from goalie.go_mesh_seq import GoalOrientedAdaptParameters
+from goalie.time_partition import TimeInterval
 
 
 class TestAdjointUtils(unittest.TestCase):
@@ -27,7 +31,7 @@ def get_qoi(mesh_seq, i):
             R = FunctionSpace(mesh_seq[i], "R", 0)
 
             def qoi():
-                return Function(R).assign(1) * dx
+                return Function(R).assign(1) * ufl.dx
 
             return qoi
 
@@ -39,7 +43,7 @@ def get_qoi(mesh_seq, i):
             R = FunctionSpace(mesh_seq[i], "R", 0)
 
             def qoi(t):
-                return Function(R).assign(1) * dx
+                return Function(R).assign(1) * ufl.dx
 
             return qoi
 
@@ -51,7 +55,7 @@ def get_qoi(mesh_seq, i):
             R = FunctionSpace(mesh_seq[i], "R", 0)
 
             def qoi():
-                return Function(R).assign(1) * dx
+                return Function(R).assign(1) * ufl.dx
 
             return qoi
 
@@ -66,7 +70,7 @@ def get_qoi(mesh_seq, i):
             R = FunctionSpace(mesh_seq[i], "R", 0)
 
             def qoi(t):
-                return Function(R).assign(1) * dx
+                return Function(R).assign(1) * ufl.dx
 
             return qoi
 
@@ -81,7 +85,7 @@ def get_qoi(mesh_seq, i):
             R = FunctionSpace(mesh_seq[i], "R", 0)
 
             def qoi(t, r):
-                return Function(R).assign(1) * dx
+                return Function(R).assign(1) * ufl.dx
 
             return qoi
 
diff --git a/test/sensors.py b/test/sensors.py
index 0f957eb..ac84792 100644
--- a/test/sensors.py
+++ b/test/sensors.py
@@ -6,11 +6,9 @@
 moving geometries. Diss. 2011.
 """
 
-import firedrake
-from ufl import *
-from utility import uniform_mesh
+import ufl
 
-__all__ = ["bowl", "hyperbolic", "multiscale", "interweaved", "mesh_for_sensors"]
+__all__ = ["bowl", "hyperbolic", "multiscale", "interweaved"]
 
 
 def bowl(*coords):
@@ -18,21 +16,16 @@ def bowl(*coords):
 
 
 def hyperbolic(x, y):
-    return conditional(
-        abs(x * y) < 2 * pi / 50, 0.01 * sin(50 * x * y), sin(50 * x * y)
+    return ufl.conditional(
+        abs(x * y) < 2 * ufl.pi / 50, 0.01 * ufl.sin(50 * x * y), ufl.sin(50 * x * y)
     )
 
 
 def multiscale(x, y):
-    return 0.1 * sin(50 * x) + atan(0.1 / (sin(5 * y) - 2 * x))
+    return 0.1 * ufl.sin(50 * x) + ufl.atan(0.1 / (ufl.sin(5 * y) - 2 * x))
 
 
 def interweaved(x, y):
-    return atan(0.1 / (sin(5 * y) - 2 * x)) + atan(0.5 / (sin(3 * y) - 7 * x))
-
-
-def mesh_for_sensors(dim, n):
-    mesh = uniform_mesh(dim, n, length=2)
-    coords = firedrake.Function(mesh.coordinates)
-    coords.interpolate(coords - as_vector([1] * dim))
-    return firedrake.Mesh(coords)
+    return ufl.atan(0.1 / (ufl.sin(5 * y) - 2 * x)) + ufl.atan(
+        0.5 / (ufl.sin(3 * y) - 7 * x)
+    )
diff --git a/test/test_error_estimation.py b/test/test_error_estimation.py
index bc74fe6..1ea19a6 100644
--- a/test/test_error_estimation.py
+++ b/test/test_error_estimation.py
@@ -1,6 +1,11 @@
 import unittest
 
-from firedrake import *
+import numpy as np
+import ufl
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace
+from firedrake.ufl_expr import TestFunction, TrialFunction
+from firedrake.utility_meshes import UnitSquareMesh, UnitTriangleMesh
 from parameterized import parameterized
 
 from goalie.error_estimation import (
@@ -38,20 +43,20 @@ def test_form_type_error(self):
         self.assertEqual(str(cm.exception), msg)
 
     def test_exterior_facet_integral(self):
-        F = self.one * ds(1) - self.one * ds(2)
+        F = self.one * ufl.ds(1) - self.one * ufl.ds(2)
         indicator = form2indicator(F)
         self.assertAlmostEqual(indicator.dat.data[0], -1)
         self.assertAlmostEqual(indicator.dat.data[1], 1)
 
     def test_interior_facet_integral(self):
-        F = avg(self.one) * dS
+        F = ufl.avg(self.one) * ufl.dS
         indicator = form2indicator(F)
-        self.assertAlmostEqual(indicator.dat.data[0], sqrt(2))
-        self.assertAlmostEqual(indicator.dat.data[1], sqrt(2))
+        self.assertAlmostEqual(indicator.dat.data[0], np.sqrt(2))
+        self.assertAlmostEqual(indicator.dat.data[1], np.sqrt(2))
 
     def test_cell_integral(self):
-        x, y = SpatialCoordinate(self.mesh)
-        F = conditional(x + y < 1, 1, 0) * dx
+        x, y = ufl.SpatialCoordinate(self.mesh)
+        F = ufl.conditional(x + y < 1, 1, 0) * ufl.dx
         indicator = form2indicator(F)
         self.assertAlmostEqual(indicator.dat.data[0], 0)
         self.assertAlmostEqual(indicator.dat.data[1], 0.5)
@@ -88,14 +93,14 @@ def test_absolute_value_type_error(self):
 
     def test_unit_time_instant(self):
         mesh_seq = self.mesh_seq(time_partition=TimeInstant("field", time=1.0))
-        mesh_seq.indicators["field"][0][0].assign(form2indicator(self.one * dx))
+        mesh_seq.indicators["field"][0][0].assign(form2indicator(self.one * ufl.dx))
         estimator = mesh_seq.error_estimate()
         self.assertAlmostEqual(estimator, 1)  # 1 * (0.5 + 0.5)
 
     @parameterized.expand([[False], [True]])
     def test_unit_time_instant_abs(self, absolute_value):
         mesh_seq = self.mesh_seq(time_partition=TimeInstant("field", time=1.0))
-        mesh_seq.indicators["field"][0][0].assign(form2indicator(-self.one * dx))
+        mesh_seq.indicators["field"][0][0].assign(form2indicator(-self.one * ufl.dx))
         estimator = mesh_seq.error_estimate(absolute_value=absolute_value)
         self.assertAlmostEqual(
             estimator, 1 if absolute_value else -1
@@ -103,7 +108,7 @@ def test_unit_time_instant_abs(self, absolute_value):
 
     def test_half_time_instant(self):
         mesh_seq = self.mesh_seq(time_partition=TimeInstant("field", time=0.5))
-        mesh_seq.indicators["field"][0][0].assign(form2indicator(self.one * dx))
+        mesh_seq.indicators["field"][0][0].assign(form2indicator(self.one * ufl.dx))
         estimator = mesh_seq.error_estimate()
         self.assertAlmostEqual(estimator, 0.5)  # 0.5 * (0.5 + 0.5)
 
@@ -111,7 +116,7 @@ def test_time_partition_same_timestep(self):
         mesh_seq = self.mesh_seq(
             time_partition=TimePartition(1.0, 2, [0.5, 0.5], ["field"])
         )
-        mesh_seq.indicators["field"][0][0].assign(form2indicator(2 * self.one * dx))
+        mesh_seq.indicators["field"][0][0].assign(form2indicator(2 * self.one * ufl.dx))
         estimator = mesh_seq.error_estimate()
         self.assertAlmostEqual(estimator, 1)  # 2 * 0.5 * (0.5 + 0.5)
 
@@ -119,7 +124,7 @@ def test_time_partition_different_timesteps(self):
         mesh_seq = self.mesh_seq(
             time_partition=TimePartition(1.0, 2, [0.5, 0.25], ["field"])
         )
-        indicator = form2indicator(self.one * dx)
+        indicator = form2indicator(self.one * ufl.dx)
         mesh_seq.indicators["field"][0][0].assign(indicator)
         mesh_seq.indicators["field"][1][0].assign(indicator)
         mesh_seq.indicators["field"][1][1].assign(indicator)
@@ -132,7 +137,7 @@ def test_time_instant_multiple_fields(self):
         mesh_seq = self.mesh_seq(
             time_partition=TimeInstant(["field1", "field2"], time=1.0)
         )
-        indicator = form2indicator(self.one * dx)
+        indicator = form2indicator(self.one * ufl.dx)
         mesh_seq.indicators["field1"][0][0].assign(indicator)
         mesh_seq.indicators["field2"][0][0].assign(indicator)
         estimator = mesh_seq.error_estimate()
@@ -148,7 +153,7 @@ def setUp(self):
         super().setUp()
         self.two = Function(self.fs, name="Dos")
         self.two.assign(2)
-        self.F = self.one * self.test * dx
+        self.F = self.one * self.test * ufl.dx
 
     def test_form_type_error(self):
         with self.assertRaises(TypeError) as cm:
@@ -248,7 +253,3 @@ def test_convert_adjoint_error_mismatch(self):
             get_dwr_indicator(self.F, self.two, test_space=test_space)
         msg = "Key 'Dos' does not exist in the test space provided."
         self.assertEqual(str(cm.exception), msg)
-
-
-if __name__ == "__main__":
-    unittest.main()  # pragma: no cover
diff --git a/test/test_function_data.py b/test/test_function_data.py
index b88375f..87576c5 100644
--- a/test/test_function_data.py
+++ b/test/test_function_data.py
@@ -3,12 +3,18 @@
 """
 
 import abc
+import os
 import unittest
 from tempfile import TemporaryDirectory
 
-from firedrake import *
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace
+from firedrake.mesh import MeshHierarchy
+from firedrake.utility_meshes import UnitTriangleMesh
 
-from goalie import *
+from goalie.function_data import AdjointSolutionData, ForwardSolutionData, IndicatorData
+from goalie.time_partition import TimePartition
+from goalie.utility import AttrDict
 
 
 class BaseTestCases:
@@ -453,7 +459,3 @@ def test_transfer_prolong(self):
                             source_function.dat.data.all()
                             == target_function.dat.data.all()
                         )
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/test/test_math.py b/test/test_math.py
index 50f6ade..5571bb0 100644
--- a/test/test_math.py
+++ b/test/test_math.py
@@ -4,7 +4,7 @@
 import scipy as sp
 from firedrake import UnitTriangleMesh
 
-from goalie.math import *
+from goalie.math import bessi0, bessk0
 
 
 class TestBessel(unittest.TestCase):
diff --git a/test/test_metric.py b/test/test_metric.py
index e877596..6c2637e 100644
--- a/test/test_metric.py
+++ b/test/test_metric.py
@@ -6,14 +6,19 @@
 
 import numpy as np
 import pytest
+import sensors
 from animate.metric import RiemannianMetric
 from animate.utility import errornorm
-from firedrake import *
+from firedrake.functionspace import TensorFunctionSpace
 from parameterized import parameterized
-from sensors import *
-from utility import uniform_mesh
+from utility import mesh_for_sensors, uniform_mesh
 
-from goalie import *
+from goalie.metric import (
+    enforce_variable_constraints,
+    ramp_complexity,
+    space_time_normalise,
+)
+from goalie.time_partition import TimeInterval, TimePartition
 
 
 class BaseClasses:
@@ -127,18 +132,18 @@ def test_target_complexity_negative_error(self):
     @pytest.mark.slow
     @parameterized.expand(
         [
-            (bowl, 1),
-            (bowl, 2),
-            (bowl, np.inf),
-            (hyperbolic, 1),
-            (hyperbolic, 2),
-            (hyperbolic, np.inf),
-            (multiscale, 1),
-            (multiscale, 2),
-            (multiscale, np.inf),
-            (interweaved, 1),
-            (interweaved, 2),
-            (interweaved, np.inf),
+            (sensors.bowl, 1),
+            (sensors.bowl, 2),
+            (sensors.bowl, np.inf),
+            (sensors.hyperbolic, 1),
+            (sensors.hyperbolic, 2),
+            (sensors.hyperbolic, np.inf),
+            (sensors.multiscale, 1),
+            (sensors.multiscale, 2),
+            (sensors.multiscale, np.inf),
+            (sensors.interweaved, 1),
+            (sensors.interweaved, 2),
+            (sensors.interweaved, np.inf),
         ]
     )
     def test_consistency(self, sensor, degree):
diff --git a/test/test_options.py b/test/test_options.py
index 849ae18..89e54dd 100644
--- a/test/test_options.py
+++ b/test/test_options.py
@@ -1,6 +1,6 @@
 import unittest
 
-from goalie.options import *
+from goalie.options import AdaptParameters, GoalOrientedAdaptParameters
 
 
 class TestAdaptParameters(unittest.TestCase):
@@ -148,7 +148,3 @@ def test_estimator_rtol_type_error(self):
             " 'str'."
         )
         self.assertEqual(str(cm.exception), msg)
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/test/test_parallel.py b/test/test_parallel.py
index 544f6d5..917a233 100644
--- a/test/test_parallel.py
+++ b/test/test_parallel.py
@@ -1,5 +1,6 @@
 import pytest
-from firedrake import *
+from firedrake.utility_meshes import UnitCubeMesh, UnitSquareMesh
+from pyop2.mpi import COMM_WORLD
 
 from goalie.mesh_seq import MeshSeq
 from goalie.time_partition import TimeInterval
diff --git a/test/test_time_partition.py b/test/test_time_partition.py
index f4a2360..2f57a5f 100644
--- a/test/test_time_partition.py
+++ b/test/test_time_partition.py
@@ -339,7 +339,3 @@ def test_time_partition_slice(self):
         self.assertAlmostEqual(tp12.timesteps, timesteps[1:3])
         self.assertEqual(tp0.field_names, self.field_names)
         self.assertEqual(tp0.field_names, tp12.field_names)
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/test/test_utility.py b/test/test_utility.py
index 9713e26..29c2d76 100644
--- a/test/test_utility.py
+++ b/test/test_utility.py
@@ -7,10 +7,12 @@
 import shutil
 import unittest
 
-from firedrake import *
+from firedrake.function import Function
+from firedrake.functionspace import FunctionSpace
 from utility import uniform_mesh
 
-from goalie import *
+from goalie.error_estimation import effectivity_index
+from goalie.utility import create_directory
 
 pointwise_norm_types = [["l1"], ["l2"], ["linf"]]
 integral_scalar_norm_types = [["L1"], ["L2"], ["L4"], ["H1"], ["HCurl"]]
@@ -78,7 +80,3 @@ def test_create_directory():
         raise OSError(
             f"Can't remove {fpath} because it isn't empty. Contents: {ls}."
         ) from e
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/test/utility.py b/test/utility.py
index 87961e6..b4800b8 100644
--- a/test/utility.py
+++ b/test/utility.py
@@ -7,6 +7,8 @@
 
 from goalie.metric import RiemannianMetric
 
+__all__ = ["uniform_mesh", "uniform_metric", "mesh_for_sensors"]
+
 
 def uniform_mesh(dim, n, length=1, **kwargs):
     args = [n] * dim + [length]
@@ -18,3 +20,10 @@ def uniform_metric(function_space, scaling):
     metric = RiemannianMetric(function_space)
     metric.interpolate(scaling * ufl.Identity(dim))
     return metric
+
+
+def mesh_for_sensors(dim, n):
+    mesh = uniform_mesh(dim, n, length=2)
+    coords = firedrake.Function(mesh.coordinates)
+    coords.interpolate(coords - ufl.as_vector([1] * dim))
+    return firedrake.Mesh(coords)