diff --git a/tests/firedrake_adjoint/test_assignment.py b/tests/firedrake_adjoint/test_assignment.py
index a9c04449..eb414836 100644
--- a/tests/firedrake_adjoint/test_assignment.py
+++ b/tests/firedrake_adjoint/test_assignment.py
@@ -77,22 +77,23 @@ def test_assign_tlm():
 def test_assign_tlm_with_constant():
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "CG", 1)
+    R = FunctionSpace(mesh, "R", 0)
 
     x = SpatialCoordinate(mesh)
     f = assemble(interpolate(x[0], V))
     g = assemble(interpolate(sin(x[0]), V))
-    c = Constant(5.0, domain=mesh)
+    c = Function(R, val=5.0)
 
     u = Function(V)
     u.interpolate(c * f**2)
 
-    c.block_variable.tlm_value = Constant(0.3, domain=mesh)
+    c.block_variable.tlm_value = Function(R, val=0.3)
     tape = get_working_tape()
     tape.evaluate_tlm()
     assert_allclose(u.block_variable.tlm_value.dat.data, 0.3 * f.dat.data ** 2)
 
     tape.reset_tlm_values()
-    c.block_variable.tlm_value = Constant(0.4, domain=mesh)
+    c.block_variable.tlm_value = Function(R, val=0.4)
     f.block_variable.tlm_value = g
     tape.evaluate_tlm()
     assert_allclose(u.block_variable.tlm_value.dat.data, 0.4 * f.dat.data ** 2 + 10. * f.dat.data * g.dat.data)
@@ -143,11 +144,11 @@ def test_assign_nonlincom():
 def test_assign_with_constant():
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "CG", 1)
-
+    R = FunctionSpace(mesh, "R", 0)
     x = SpatialCoordinate(mesh)
     f = assemble(interpolate(x[0], V))
-    c = Constant(3.0, domain=mesh)
-    d = Constant(2.0, domain=mesh)
+    c = Function(R, val=3.0)
+    d = Function(R, val=2.0)
     u = Function(V)
 
     u.assign(c*f+d**3)
@@ -198,9 +199,9 @@ def test_assign_nonlin_changing():
 def test_assign_constant_scale():
     mesh = UnitSquareMesh(10, 10)
     V = VectorFunctionSpace(mesh, "CG", 1)
-
+    R = FunctionSpace(mesh, "R", 0)
     f = Function(V)
-    c = Constant(2.0, domain=mesh)
+    c = Function(R, val=2.0)
     x, y = SpatialCoordinate(mesh)
     g = assemble(interpolate(as_vector([sin(y)+x, cos(x)*y]), V))
 
@@ -209,8 +210,7 @@ def test_assign_constant_scale():
     J = assemble(inner(f, f) ** 2  * dx)
 
     rf = ReducedFunctional(J, Control(c))
-    h = Constant(0.1)
-    r = taylor_to_dict(rf, c, h)
+    r = taylor_to_dict(rf, c, Constant(0.1))
 
     assert min(r["R0"]["Rate"]) > 0.9
     assert min(r["R1"]["Rate"]) > 1.9
diff --git a/tests/firedrake_adjoint/test_disk_checkpointing.py b/tests/firedrake_adjoint/test_disk_checkpointing.py
index e7daaf50..76119d54 100644
--- a/tests/firedrake_adjoint/test_disk_checkpointing.py
+++ b/tests/firedrake_adjoint/test_disk_checkpointing.py
@@ -19,7 +19,7 @@ def adjoint_example(fine, coarse):
     # AssembleBlock
     m = assemble(interpolate(sin(4*pi*x)*cos(4*pi*y), cg_space))
 
-    u, v = w.split()
+    u, v = w.subfunctions
     # FunctionAssignBlock, FunctionMergeBlock
     v.assign(m)
     # FunctionSplitBlock, GenericSolveBlock
diff --git a/tests/firedrake_adjoint/test_hessian.py b/tests/firedrake_adjoint/test_hessian.py
index a7043e8a..876823f2 100644
--- a/tests/firedrake_adjoint/test_hessian.py
+++ b/tests/firedrake_adjoint/test_hessian.py
@@ -107,8 +107,8 @@ def test_function():
 
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "Lagrange", 2)
-
-    c = Constant(4, domain=mesh)
+    R = FunctionSpace(mesh, "R", 0)
+    c = Function(R, val=4)
     control_c = Control(c)
     f = Function(V)
     f.vector()[:] = 3
@@ -116,7 +116,7 @@ def test_function():
 
     u = Function(V)
     v = TestFunction(V)
-    bc = DirichletBC(V, Constant(1, domain=mesh), "on_boundary")
+    bc = DirichletBC(V, Function(R, val=1), "on_boundary")
 
     F = inner(grad(u), grad(v)) * dx + u**2*v*dx - f ** 2 * v * dx - c**2*v*dx
     solve(F == 0, u, bc)
@@ -127,7 +127,7 @@ def test_function():
     dJdc, dJdf = compute_gradient(J, [control_c, control_f])
 
     # Step direction for derivatives and convergence test
-    h_c = Constant(1.0, domain=mesh)
+    h_c = Function(R, val=1.0) 
     h_f = Function(V)
     h_f.vector()[:] = 10*rng.random(V.dim())
 
@@ -148,13 +148,13 @@ def test_nonlinear():
 
     mesh = UnitSquareMesh(10, 10)
     V = FunctionSpace(mesh, "Lagrange", 1)
-
+    R = FunctionSpace(mesh, "R", 0)
     f = Function(V)
     f.vector()[:] = 5
 
     u = Function(V)
     v = TestFunction(V)
-    bc = DirichletBC(V, Constant(1, domain=mesh), "on_boundary")
+    bc = DirichletBC(V, Function(R, val=1), "on_boundary")
 
     F = inner(grad(u), grad(v)) * dx - u**2*v*dx - f * v * dx
     solve(F == 0, u, bc)
diff --git a/tests/firedrake_adjoint/test_optimisation.py b/tests/firedrake_adjoint/test_optimisation.py
index 9f10b5a5..40ad346e 100644
--- a/tests/firedrake_adjoint/test_optimisation.py
+++ b/tests/firedrake_adjoint/test_optimisation.py
@@ -9,9 +9,10 @@
 def test_optimisation_constant_control():
     """This tests a list of controls in a minimisation (through scipy L-BFGS-B)"""
     mesh = UnitSquareMesh(1, 1)
+    R = FunctionSpace(mesh, "R", 0)
 
     n = 3
-    x = [Constant(0., domain=mesh) for i in range(n)]
+    x = [Function(R) for i in range(n)]
     c = [Control(xi) for xi in x]
 
     # Rosenbrock function https://en.wikipedia.org/wiki/Rosenbrock_function
diff --git a/tests/firedrake_adjoint/test_projection.py b/tests/firedrake_adjoint/test_projection.py
index df19e897..60fe0492 100644
--- a/tests/firedrake_adjoint/test_projection.py
+++ b/tests/firedrake_adjoint/test_projection.py
@@ -135,36 +135,36 @@ def test_project_nonlin_changing():
 def test_self_project():
     mesh = UnitSquareMesh(1,1)
     V = FunctionSpace(mesh, "CG", 1)
+    R = FunctionSpace(mesh, "R", 0)
     u = Function(V)
-    c = Constant(1., domain=mesh)
+    c = Function(R, val=1.)
     u.project(u+c)
     J = assemble(u**2*dx)
     rf = ReducedFunctional(J, Control(c))
-    h = Constant(0.1)
-    assert taylor_test(rf, Constant(2., domain=mesh), h)
+    assert taylor_test(rf,  Function(R, val=2.), Constant(0.1))
 
 def test_self_project_function():
     mesh = UnitSquareMesh(1,1)
     V = FunctionSpace(mesh, "CG", 1)
+    R = FunctionSpace(mesh, "R", 0)
     u = Function(V)
-    c = Constant(1., domain=mesh)
+    c = Function(R, val=1.)
     project(u+c, u)
     project(u+c*u**2, u)
     J = assemble(u**2*dx)
     rf = ReducedFunctional(J, Control(c))
-    h = Constant(0.1)
-    assert taylor_test(rf, Constant(3., domain=mesh), h)
+    assert taylor_test(rf, Function(R, val=3.), Constant(0.1))
 
 def test_project_to_function_space():
     mesh = UnitSquareMesh(1,1)
     V = FunctionSpace(mesh, "CG", 1)
     W = FunctionSpace(mesh, "DG", 1)
+    R = FunctionSpace(mesh, "R", 0)
     u = Function(V)
     x = SpatialCoordinate(mesh)
     u.interpolate(x[0])
-    c = Constant(1., domain=mesh)
+    c = Function(R, val=1.)
     w = project((u+c)*u, W)
     J = assemble(w**2*dx)
     rf = ReducedFunctional(J, Control(c))
-    h = Constant(0.1)
-    assert taylor_test(rf, Constant(1., domain=mesh), h)
+    assert taylor_test(rf, Function(R, val=1.), Constant(0.1)) > 1.9
diff --git a/tests/firedrake_adjoint/test_reduced_functional.py b/tests/firedrake_adjoint/test_reduced_functional.py
index ae918d96..e9c827e0 100644
--- a/tests/firedrake_adjoint/test_reduced_functional.py
+++ b/tests/firedrake_adjoint/test_reduced_functional.py
@@ -9,21 +9,22 @@
 def test_constant():
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "Lagrange", 1)
+    R = FunctionSpace(mesh, "R", 0)
 
-    c = Constant(1, domain=mesh)
+    c = Function(R, val=1)
     f = Function(V)
     f.vector()[:] = 1
 
     u = Function(V)
     v = TestFunction(V)
-    bc = DirichletBC(V, Constant(1, domain=mesh), "on_boundary")
+    bc = DirichletBC(V, Function(R, val=1), "on_boundary")
 
     F = inner(grad(u), grad(v))*dx - f**2*v*dx
     solve(F == 0, u, bc)
 
     J = assemble(c**2*u*dx)
     Jhat = ReducedFunctional(J, Control(c))
-    assert taylor_test(Jhat, c, Constant(1, domain=mesh)) > 1.9
+    assert taylor_test(Jhat, c, Function(R, val=1)) > 1.9
 
 
 def test_function():
@@ -54,6 +55,7 @@ def test_wrt_function_dirichlet_boundary(control):
     mesh = UnitSquareMesh(10,10)
 
     V = FunctionSpace(mesh,"CG",1)
+    R = FunctionSpace(mesh,"R",0)
     u = TrialFunction(V)
     u_ = Function(V)
     v = TestFunction(V)
@@ -64,8 +66,8 @@ def test_wrt_function_dirichlet_boundary(control):
     bc2 = DirichletBC(V, 2, 2)
     bc = [bc1,bc2]
 
-    g1 = Constant(2, domain=mesh)
-    g2 = Constant(1, domain=mesh)
+    g1 = Function(R, val=2)
+    g2 = Function(R, val=1)
     f = Function(V)
     f.vector()[:] = 10
 
@@ -166,10 +168,11 @@ def test_mixed_boundary():
 def test_assemble_recompute():
     mesh = UnitSquareMesh(10, 10)
     V = FunctionSpace(mesh, "CG", 1)
+    R = FunctionSpace(mesh, "R", 0)
 
     f = Function(V)
     f.vector()[:] = 2
-    expr = Constant(assemble(f**2*dx), domain=mesh)
+    expr = Function(R).assign(assemble(f**2*dx))
     J = assemble(expr**2*dx(domain=mesh))
     Jhat = ReducedFunctional(J, Control(f))
 
diff --git a/tests/firedrake_adjoint/test_solving.py b/tests/firedrake_adjoint/test_solving.py
index 45edb07e..28b011ee 100644
--- a/tests/firedrake_adjoint/test_solving.py
+++ b/tests/firedrake_adjoint/test_solving.py
@@ -9,14 +9,14 @@
 def test_linear_problem():
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "Lagrange", 1)
-
+    R = FunctionSpace(mesh, "R", 0)
     f = Function(V)
     f.vector()[:] = 1
 
     u = TrialFunction(V)
     u_ = Function(V)
     v = TestFunction(V)
-    bc = DirichletBC(V, Constant(1, domain=mesh), "on_boundary")
+    bc = DirichletBC(V, Function(R, val=1), "on_boundary")
 
     def J(f):
         a = inner(grad(u), grad(v))*dx
@@ -60,13 +60,13 @@ def test_nonlinear_problem():
     """This tests whether nullspace and solver_parameters are passed on in adjoint solves"""
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "Lagrange", 1)
-
+    R = FunctionSpace(mesh, "R", 0)
     f = Function(V)
     f.vector()[:] = 1
 
     u = Function(V)
     v = TestFunction(V)
-    bc = DirichletBC(V, Constant(1, domain=mesh), "on_boundary")
+    bc = DirichletBC(V, Function(R, val=1), "on_boundary")
 
     def J(f):
         a = f*inner(grad(u), grad(v))*dx + u**2*v*dx - f*v*dx
@@ -160,6 +160,7 @@ def test_wrt_function_neumann_boundary():
     mesh = UnitSquareMesh(10,10)
 
     V = FunctionSpace(mesh,"CG",1)
+    R = FunctionSpace(mesh,"R",0)
     u = TrialFunction(V)
     u_ = Function(V)
     v = TestFunction(V)
@@ -169,8 +170,8 @@ def test_wrt_function_neumann_boundary():
     bc2 = DirichletBC(V, 2, 2)
     bc = [bc1,bc2]
 
-    g1 = Constant(2, domain=mesh)
-    g2 = Constant(1, domain=mesh)
+    g1 = Function(R, val=2)
+    g2 = Function(R, val=1)
     f = Function(V)
     f.vector()[:] = 10
 
@@ -188,13 +189,14 @@ def J(g1):
 def test_wrt_constant():
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "Lagrange", 1)
+    R = FunctionSpace(mesh, "R", 0)
 
-    c = Constant(1, domain=mesh)
+    c = Function(R, val=1)
 
     u = TrialFunction(V)
     u_ = Function(V)
     v = TestFunction(V)
-    bc = DirichletBC(V, Constant(1, domain=mesh), "on_boundary")
+    bc = DirichletBC(V, Function(R, val=1), "on_boundary")
 
     def J(c):
         a = inner(grad(u), grad(v))*dx
@@ -209,6 +211,7 @@ def test_wrt_constant_neumann_boundary():
     mesh = UnitSquareMesh(10,10)
 
     V = FunctionSpace(mesh,"CG",1)
+    R = FunctionSpace(mesh,"R",0)
     u = TrialFunction(V)
     u_ = Function(V)
     v = TestFunction(V)
@@ -218,8 +221,8 @@ def test_wrt_constant_neumann_boundary():
     bc2 = DirichletBC(V, 2, 2)
     bc = [bc1,bc2]
 
-    g1 = Constant(2, domain=mesh)
-    g2 = Constant(1, domain=mesh)
+    g1 = Function(R, val=2)
+    g2 = Function(R, val=1)
     f = Function(V)
     f.vector()[:] = 10
 
@@ -240,6 +243,7 @@ def test_time_dependent():
 
     # Defining function space, test and trial functions
     V = FunctionSpace(mesh,"CG",1)
+    R = FunctionSpace(mesh,"R",0)
     u = TrialFunction(V)
     u_ = Function(V)
     v = TestFunction(V)
@@ -252,7 +256,7 @@ def test_time_dependent():
     # Some variables
     T = 0.2
     dt = 0.1
-    f = Constant(1, domain=mesh)
+    f = Function(R, val=1)
 
     def J(f):
         u_1 = Function(V)
@@ -277,11 +281,12 @@ def test_two_nonlinear_solves():
     # regression test for firedrake issue #1841
     mesh = UnitSquareMesh(1,1)
     V = FunctionSpace(mesh, "CG", 1)
+    R = FunctionSpace(mesh, "R", 0)
     v = TestFunction(V)
     u0 = Function(V)
     u1 = Function(V)
 
-    ui = Constant(2.0, domain=mesh)
+    ui = Function(R, val=2.0)
     c = Control(ui)
     u0.assign(ui)
     F = dot(v, (u1-u0))*dx - dot(v, u0*u1)*dx
diff --git a/tests/firedrake_adjoint/test_split.py b/tests/firedrake_adjoint/test_split.py
index 2c3e27cc..cfe5692e 100644
--- a/tests/firedrake_adjoint/test_split.py
+++ b/tests/firedrake_adjoint/test_split.py
@@ -30,7 +30,7 @@ def main(ic, fnsplit=True):
     v = TestFunction(V2)
 
     if fnsplit:
-        ic_u = ic.split()[0]
+        ic_u = ic.sub(0)
     else:
         ic_u = split(ic)[0]
 
@@ -88,8 +88,9 @@ def test_fn_split_no_annotate(Z):
     w = TrialFunction(V2)
     v = TestFunction(V2)
 
-    ic_u = ic.split(annotate=True)[0]
-    ic_uv = ic.split(annotate=False)[0]
+    ic_u = ic.sub(0)
+    with stop_annotating():
+        ic_uv = ic.sub(0)
 
     mass = inner(w, v) * dx
     rhs = inner(ic_u, v) * dx
@@ -108,23 +109,24 @@ def test_fn_split_no_annotate(Z):
 
 def test_split_subvariables_update(Z):
     z = Function(Z)
-    u = z.split()[0]
+    u = z.sub(0)
     u.project(Constant(1.))
     assert np.allclose(z.sub(0).vector().dat.data, u.vector().dat.data)
 
 def test_merge_blocks():
     mesh = UnitSquareMesh(1,1)
     V = FunctionSpace(mesh, 'CG', 1)
+    R = FunctionSpace(mesh, 'R', 0)
     W = V * V
     w = Function(W)
-    w1, w2 = w.split()
-    w1_const = Constant(0.1, domain=mesh)
-    w2_const = Constant(0.2, domain=mesh)
+    w1, w2 = w.subfunctions
+    w1_const = Function(R, val=0.1)
+    w2_const = Function(R, val=0.2)
     w1.project(w1_const)
     w2.project(w2_const)
     J = assemble(w1*w1*dx)
     c = Control(w1_const)
     rf = ReducedFunctional(J, c)
-    assert taylor_test(rf, Constant(0.3, domain=mesh), Constant(0.01, domain=mesh)) > 1.95
+    assert taylor_test(rf, Function(R, val=0.3), Function(R, val=0.01)) > 1.95
 
 
diff --git a/tests/firedrake_adjoint/test_tlm.py b/tests/firedrake_adjoint/test_tlm.py
index 03321398..e80ecd5f 100644
--- a/tests/firedrake_adjoint/test_tlm.py
+++ b/tests/firedrake_adjoint/test_tlm.py
@@ -46,8 +46,8 @@ def test_tlm_bc():
     set_working_tape(tape)
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "Lagrange", 1)
-
-    c = Constant(1, domain=mesh)
+    R = FunctionSpace(mesh, "R", 0)
+    c = Function(R, val=1)
     f = Function(V)
     f.vector()[:] = 1
 
@@ -63,7 +63,7 @@ def test_tlm_bc():
 
     # Need to specify the domain for the constant as `ufl.action`, which requires `ufl.Constant`
     # to have a function space, will be applied on the tlm value.
-    c.block_variable.tlm_value = Constant(1, domain=mesh)
+    c.block_variable.tlm_value = Function(R, val=1)
     tape.evaluate_tlm()
 
     assert (taylor_test(Jhat, Constant(c), Constant(1), dJdm=J.block_variable.tlm_value) > 1.9)
@@ -212,9 +212,10 @@ def test_projection():
     set_working_tape(tape)
     mesh = UnitSquareMesh(10, 10)
     V = FunctionSpace(mesh, "CG", 1)
+    R = FunctionSpace(mesh, "R", 0)
 
     bc = DirichletBC(V, Constant(1), "on_boundary")
-    k = Constant(2.0, domain=mesh)
+    k = Function(R, val=2)
     x, y = SpatialCoordinate(mesh)
     expr = sin(k*x)
     f = project(expr, V)