diff --git a/demos/boussinesq/boussinesq.py.rst b/demos/boussinesq/boussinesq.py.rst
index edfdf3c1a5..5cc6708cf0 100644
--- a/demos/boussinesq/boussinesq.py.rst
+++ b/demos/boussinesq/boussinesq.py.rst
@@ -184,7 +184,7 @@ implements a boundary condition that fixes a field at a single point. ::
 
            # Take the basis function with the largest abs value at bc_point
            v = TestFunction(V)
-           F = assemble(Interpolate(inner(v, v), Fvom))
+           F = assemble(interpolate(inner(v, v), Fvom))
            with F.dat.vec as Fvec:
                max_index, _ = Fvec.max()
            nodes = V.dof_dset.lgmap.applyInverse([max_index])
diff --git a/demos/multicomponent/multicomponent.py.rst b/demos/multicomponent/multicomponent.py.rst
index bf74e5d2e0..7a29d6ef1c 100644
--- a/demos/multicomponent/multicomponent.py.rst
+++ b/demos/multicomponent/multicomponent.py.rst
@@ -521,7 +521,7 @@ mathematically valid to do this)::
 
            # Take the basis function with the largest abs value at bc_point
            v = TestFunction(V)
-           F = assemble(Interpolate(inner(v, v), Fvom))
+           F = assemble(interpolate(inner(v, v), Fvom))
            with F.dat.vec as Fvec:
                max_index, _ = Fvec.max()
            nodes = V.dof_dset.lgmap.applyInverse([max_index])
diff --git a/firedrake/cofunction.py b/firedrake/cofunction.py
index f0fda10f63..9ee867c622 100644
--- a/firedrake/cofunction.py
+++ b/firedrake/cofunction.py
@@ -318,7 +318,7 @@ def interpolate(self,
         Parameters
         ----------
         expression
-            A dual UFL expression to interpolate.
+            A UFL BaseForm to adjoint interpolate.
         ad_block_tag
             An optional string for tagging the resulting assemble
             block on the Pyadjoint tape.
@@ -331,9 +331,9 @@ def interpolate(self,
         firedrake.cofunction.Cofunction
             Returns `self`
         """
-        from firedrake import interpolation, assemble
+        from firedrake import interpolate, assemble
         v, = self.arguments()
-        interp = interpolation.Interpolate(v, expression, **kwargs)
+        interp = interpolate(v, expression, **kwargs)
         return assemble(interp, tensor=self, ad_block_tag=ad_block_tag)
 
     @property
diff --git a/firedrake/external_operators/point_expr_operator.py b/firedrake/external_operators/point_expr_operator.py
index 3aa40e1d5b..4e7183e47f 100644
--- a/firedrake/external_operators/point_expr_operator.py
+++ b/firedrake/external_operators/point_expr_operator.py
@@ -5,7 +5,7 @@
 
 import firedrake.ufl_expr as ufl_expr
 from firedrake.assemble import assemble
-from firedrake.interpolation import Interpolate
+from firedrake.interpolation import interpolate
 from firedrake.external_operators import AbstractExternalOperator, assemble_method
 
 
@@ -58,7 +58,7 @@ def assemble_operator(self, *args, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
         if len(V) < 2:
-            interp = Interpolate(expr, self.function_space())
+            interp = interpolate(expr, self.function_space())
             return assemble(interp)
         # Interpolation of UFL expressions for mixed functions is not yet supported
         # -> `Function.assign` might be enough in some cases.
@@ -72,7 +72,7 @@ def assemble_operator(self, *args, **kwargs):
     def assemble_Jacobian_action(self, *args, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
-        interp = Interpolate(expr, V)
+        interp = interpolate(expr, V)
 
         u, = [e for i, e in enumerate(self.ufl_operands) if self.derivatives[i] == 1]
         w = self.argument_slots()[-1]
@@ -83,7 +83,7 @@ def assemble_Jacobian_action(self, *args, **kwargs):
     def assemble_Jacobian(self, *args, assembly_opts, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
-        interp = Interpolate(expr, V)
+        interp = interpolate(expr, V)
 
         u, = [e for i, e in enumerate(self.ufl_operands) if self.derivatives[i] == 1]
         jac = ufl_expr.derivative(interp, u)
@@ -99,7 +99,7 @@ def assemble_Jacobian_adjoint(self, *args, assembly_opts, **kwargs):
     def assemble_Jacobian_adjoint_action(self, *args, **kwargs):
         V = self.function_space()
         expr = as_ufl(self.expr(*self.ufl_operands))
-        interp = Interpolate(expr, V)
+        interp = interpolate(expr, V)
 
         u, = [e for i, e in enumerate(self.ufl_operands) if self.derivatives[i] == 1]
         ustar = self.argument_slots()[0]
diff --git a/firedrake/function.py b/firedrake/function.py
index a628ac6599..38d6b6e76b 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -382,9 +382,9 @@ def interpolate(self,
         firedrake.function.Function
             Returns `self`
         """
-        from firedrake import interpolation, assemble
+        from firedrake import interpolate, assemble
         V = self.function_space()
-        interp = interpolation.Interpolate(expression, V, **kwargs)
+        interp = interpolate(expression, V, **kwargs)
         return assemble(interp, tensor=self, ad_block_tag=ad_block_tag)
 
     def zero(self, subset=None):
@@ -697,7 +697,7 @@ def __init__(self, domain, point):
         self.point = point
 
     def __str__(self):
-        return "domain %s does not contain point %s" % (self.domain, self.point)
+        return f"Domain {self.domain} does not contain point {self.point}"
 
 
 class PointEvaluator:
@@ -712,7 +712,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             The mesh on which to embed the points.
         points : numpy.ndarray | list
             Array or list of points to evaluate at.
-        tolerance : float | None
+        tolerance : Optional[float]
             Tolerance to use when checking if a point is in a cell.
             If ``None`` (the default), the ``tolerance`` of the ``mesh`` is used.
         missing_points_behaviour : str
diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py
index 4fb288fb4d..fb60ed7da2 100644
--- a/firedrake/interpolation.py
+++ b/firedrake/interpolation.py
@@ -4,15 +4,15 @@
 import abc
 import warnings
 from collections.abc import Iterable
-from typing import Literal
 from functools import partial, singledispatch
-from typing import Hashable
+from typing import Hashable, Literal
 
 import FIAT
 import ufl
 import finat.ufl
-from ufl.algorithms import extract_arguments, extract_coefficients, replace
+from ufl.algorithms import extract_arguments, extract_coefficients
 from ufl.domain import as_domain, extract_unique_domain
+from ufl.duals import is_dual
 
 from pyop2 import op2
 from pyop2.caching import memory_and_disk_cache
@@ -25,13 +25,11 @@
 import finat
 
 import firedrake
-import firedrake.bcs
 from firedrake import tsfc_interface, utils, functionspaceimpl
 from firedrake.ufl_expr import Argument, Coargument, action, adjoint as expr_adjoint
 from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshMissingPointsError, VertexOnlyMeshTopology
 from firedrake.petsc import PETSc
 from firedrake.halo import _get_mtype as get_dat_mpi_type
-from firedrake.cofunction import Cofunction
 from mpi4py import MPI
 
 from pyadjoint import stop_annotating, no_annotations
@@ -48,7 +46,7 @@
 
 class Interpolate(ufl.Interpolate):
 
-    def __init__(self, expr, v,
+    def __init__(self, expr, V,
                  subset=None,
                  access=None,
                  allow_missing_dofs=False,
@@ -60,7 +58,7 @@ def __init__(self, expr, v,
         ----------
         expr : ufl.core.expr.Expr or ufl.BaseForm
                The UFL expression to interpolate.
-        v : firedrake.functionspaceimpl.WithGeometryBase or firedrake.ufl_expr.Coargument
+        V : firedrake.functionspaceimpl.WithGeometryBase or firedrake.ufl_expr.Coargument
             The function space to interpolate into or the coargument defined
             on the dual of the function space to interpolate into.
         subset : pyop2.types.set.Subset
@@ -95,20 +93,18 @@ def __init__(self, expr, v,
                 between a VOM and its input ordering. Defaults to ``True`` which uses SF broadcast
                 and reduce operations.
         """
-        # Check function space
         expr = ufl.as_ufl(expr)
-        if isinstance(v, functionspaceimpl.WithGeometry):
-            expr_args = extract_arguments(expr)
-            is_adjoint = len(expr_args) and expr_args[0].number() == 0
-            v = Argument(v.dual(), 1 if is_adjoint else 0)
+        if isinstance(V, functionspaceimpl.WithGeometry):
+            # Need to create a Firedrake Argument so that it has a .function_space() method
+            expr_arg_numbers = {arg.number() for arg in extract_arguments(expr) if not is_dual(arg)}
+            is_adjoint = len(expr_arg_numbers) and expr_arg_numbers == {0}
+            V = Argument(V.dual(), 1 if is_adjoint else 0)
 
-        V = v.arguments()[0].function_space()
-        if len(expr.ufl_shape) != len(V.value_shape):
-            raise RuntimeError(f'Rank mismatch: Expression rank {len(expr.ufl_shape)}, FunctionSpace rank {len(V.value_shape)}')
+        target_shape = V.arguments()[0].function_space().value_shape
+        if expr.ufl_shape != target_shape:
+            raise ValueError(f"Shape mismatch: Expression shape {expr.ufl_shape}, FunctionSpace shape {target_shape}.")
 
-        if expr.ufl_shape != V.value_shape:
-            raise RuntimeError('Shape mismatch: Expression shape {expr.ufl_shape}, FunctionSpace shape {V.value_shape}')
-        super().__init__(expr, v)
+        super().__init__(expr, V)
 
         # -- Interpolate data (e.g. `subset` or `access`) -- #
         self.interp_data = {"subset": subset,
@@ -174,32 +170,10 @@ def interpolate(expr, V, subset=None, access=None, allow_missing_dofs=False, def
        reduction (hence using MIN will compute the MIN between the
        existing values and any new values).
     """
-    if isinstance(V, (Cofunction, Coargument)):
-        dual_arg = V
-    elif isinstance(V, ufl.BaseForm):
-        rank = len(V.arguments())
-        if rank == 1:
-            dual_arg = V
-        else:
-            raise TypeError(f"Expected a one-form, provided form had {rank} arguments")
-    elif isinstance(V, functionspaceimpl.WithGeometry):
-        dual_arg = Coargument(V.dual(), 0)
-        expr_args = extract_arguments(ufl.as_ufl(expr))
-        if expr_args and expr_args[0].number() == 0:
-            warnings.warn("Passing argument numbered 0 in expression for forward interpolation is deprecated. "
-                          "Use a TrialFunction in the expression.")
-            v, = expr_args
-            expr = replace(expr, {v: v.reconstruct(number=1)})
-    else:
-        raise TypeError(f"V must be a FunctionSpace, Cofunction, Coargument or one-form, not a {type(V).__name__}")
-
-    interp = Interpolate(expr, dual_arg,
-                         subset=subset, access=access,
-                         allow_missing_dofs=allow_missing_dofs,
-                         default_missing_val=default_missing_val,
-                         matfree=matfree)
-
-    return interp
+    return Interpolate(
+        expr, V, subset=subset, access=access, allow_missing_dofs=allow_missing_dofs,
+        default_missing_val=default_missing_val, matfree=matfree
+    )
 
 
 class Interpolator(abc.ABC):
@@ -528,7 +502,7 @@ def __init__(
 
         from firedrake.assemble import assemble
         V_dest_vec = firedrake.VectorFunctionSpace(dest_mesh, ufl_scalar_element)
-        f_dest_node_coords = Interpolate(dest_mesh.coordinates, V_dest_vec)
+        f_dest_node_coords = interpolate(dest_mesh.coordinates, V_dest_vec)
         f_dest_node_coords = assemble(f_dest_node_coords)
         dest_node_coords = f_dest_node_coords.dat.data_ro.reshape(-1, dest_mesh_gdim)
         try:
@@ -553,7 +527,7 @@ def __init__(
         else:
             fs_type = partial(firedrake.TensorFunctionSpace, shape=shape)
         P0DG_vom = fs_type(self.vom_dest_node_coords_in_src_mesh, "DG", 0)
-        self.point_eval_interpolate = Interpolate(self.expr_renumbered, P0DG_vom)
+        self.point_eval_interpolate = interpolate(self.expr_renumbered, P0DG_vom)
         # The parallel decomposition of the nodes of V_dest in the DESTINATION
         # mesh (dest_mesh) is retrieved using the input_ordering attribute of the
         # VOM. This again is an interpolation operation, which, under the hood
@@ -561,7 +535,7 @@ def __init__(
         P0DG_vom_i_o = fs_type(
             self.vom_dest_node_coords_in_src_mesh.input_ordering, "DG", 0
         )
-        self.to_input_ordering_interpolate = Interpolate(
+        self.to_input_ordering_interpolate = interpolate(
             firedrake.TrialFunction(P0DG_vom), P0DG_vom_i_o
         )
         # The P0DG function outputted by the above interpolation has the
diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 87757a2e54..14ebf94d0e 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -4140,7 +4140,7 @@ def _parent_mesh_embedding(
     # nessesary, to other processes.
     P0DG = functionspace.FunctionSpace(parent_mesh, "DG", 0)
     with stop_annotating():
-        visible_ranks = interpolation.Interpolate(
+        visible_ranks = interpolation.interpolate(
             constant.Constant(parent_mesh.comm.rank), P0DG
         )
         visible_ranks = assemble(visible_ranks).dat.data_ro_with_halos.real
diff --git a/firedrake/mg/utils.py b/firedrake/mg/utils.py
index 37832b64dc..886cc7530c 100644
--- a/firedrake/mg/utils.py
+++ b/firedrake/mg/utils.py
@@ -143,7 +143,7 @@ def physical_node_locations(V):
         Vc = V.collapse().reconstruct(element=finat.ufl.VectorElement(element, dim=mesh.geometric_dimension()))
 
         # FIXME: This is unsafe for DG coordinates and CG target spaces.
-        locations = firedrake.assemble(firedrake.Interpolate(firedrake.SpatialCoordinate(mesh), Vc))
+        locations = firedrake.assemble(firedrake.interpolate(firedrake.SpatialCoordinate(mesh), Vc))
         return cache.setdefault(key, locations)
 
 
diff --git a/firedrake/preconditioners/gtmg.py b/firedrake/preconditioners/gtmg.py
index 6ce73cd6b4..2ac5df9a5d 100644
--- a/firedrake/preconditioners/gtmg.py
+++ b/firedrake/preconditioners/gtmg.py
@@ -4,7 +4,7 @@
 from firedrake.petsc import PETSc
 from firedrake.preconditioners.base import PCBase
 from firedrake.parameters import parameters
-from firedrake.interpolation import Interpolate
+from firedrake.interpolation import interpolate
 from firedrake.solving_utils import _SNESContext
 from firedrake.matrix_free.operators import ImplicitMatrixContext
 import firedrake.dmhooks as dmhooks
@@ -155,7 +155,7 @@ def initialize(self, pc):
             # Create interpolation matrix from coarse space to fine space
             fine_space = ctx.J.arguments()[0].function_space()
             coarse_test, coarse_trial = coarse_operator.arguments()
-            interp = assemble(Interpolate(coarse_trial, fine_space))
+            interp = assemble(interpolate(coarse_trial, fine_space))
             interp_petscmat = interp.petscmat
         restr_petscmat = appctx.get("restriction_matrix", None)
 
diff --git a/firedrake/preconditioners/hypre_ads.py b/firedrake/preconditioners/hypre_ads.py
index 89c10dc438..98443f2c75 100644
--- a/firedrake/preconditioners/hypre_ads.py
+++ b/firedrake/preconditioners/hypre_ads.py
@@ -1,7 +1,7 @@
 from firedrake.preconditioners.base import PCBase
 from firedrake.petsc import PETSc
 from firedrake.function import Function
-from firedrake.ufl_expr import TestFunction
+from firedrake.ufl_expr import TrialFunction
 from firedrake.dmhooks import get_function_space
 from firedrake.preconditioners.hypre_ams import chop
 from firedrake.interpolation import interpolate
@@ -31,12 +31,12 @@ def initialize(self, obj):
         NC1 = V.reconstruct(family="N1curl" if mesh.ufl_cell().is_simplex() else "NCE", degree=1)
         G_callback = appctx.get("get_gradient", None)
         if G_callback is None:
-            G = chop(assemble(interpolate(grad(TestFunction(P1)), NC1)).petscmat)
+            G = chop(assemble(interpolate(grad(TrialFunction(P1)), NC1)).petscmat)
         else:
             G = G_callback(P1, NC1)
         C_callback = appctx.get("get_curl", None)
         if C_callback is None:
-            C = chop(assemble(interpolate(curl(TestFunction(NC1)), V)).petscmat)
+            C = chop(assemble(interpolate(curl(TrialFunction(NC1)), V)).petscmat)
         else:
             C = C_callback(NC1, V)
 
diff --git a/firedrake/preconditioners/hypre_ams.py b/firedrake/preconditioners/hypre_ams.py
index 9a59702af4..594fe88590 100644
--- a/firedrake/preconditioners/hypre_ams.py
+++ b/firedrake/preconditioners/hypre_ams.py
@@ -2,7 +2,7 @@
 from firedrake.preconditioners.base import PCBase
 from firedrake.petsc import PETSc
 from firedrake.function import Function
-from firedrake.ufl_expr import TestFunction
+from firedrake.ufl_expr import TrialFunction
 from firedrake.dmhooks import get_function_space
 from firedrake.utils import complex_mode
 from firedrake.interpolation import interpolate
@@ -51,7 +51,7 @@ def initialize(self, obj):
         P1 = V.reconstruct(family="Lagrange", degree=1)
         G_callback = appctx.get("get_gradient", None)
         if G_callback is None:
-            G = chop(assemble(interpolate(grad(TestFunction(P1)), V)).petscmat)
+            G = chop(assemble(interpolate(grad(TrialFunction(P1)), V)).petscmat)
         else:
             G = G_callback(P1, V)
 
diff --git a/firedrake/preconditioners/patch.py b/firedrake/preconditioners/patch.py
index c139c07f7a..2fbe469e0f 100644
--- a/firedrake/preconditioners/patch.py
+++ b/firedrake/preconditioners/patch.py
@@ -4,7 +4,7 @@
 from firedrake.solving_utils import _SNESContext
 from firedrake.utils import cached_property, complex_mode, IntType
 from firedrake.dmhooks import get_appctx, push_appctx, pop_appctx
-from firedrake.interpolation import Interpolate
+from firedrake.interpolation import interpolate
 
 from collections import namedtuple
 import operator
@@ -660,7 +660,7 @@ def sort_entities(self, dm, axis, dir, ndiv=None, divisions=None):
             # with access descriptor MAX to define a consistent opinion
             # about where the vertices are.
             CGk = V.reconstruct(family="Lagrange")
-            coordinates = assemble(Interpolate(coordinates, CGk, access=op2.MAX))
+            coordinates = assemble(interpolate(coordinates, CGk, access=op2.MAX))
 
         select = partial(select_entity, dm=dm, exclude="pyop2_ghost")
         entities = [(p, self.coords(dm, p, coordinates)) for p in
diff --git a/firedrake/pyplot/mpl.py b/firedrake/pyplot/mpl.py
index 3cf010a1c9..d6a7aa5112 100644
--- a/firedrake/pyplot/mpl.py
+++ b/firedrake/pyplot/mpl.py
@@ -18,7 +18,7 @@
 import mpl_toolkits.mplot3d
 from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection
 from math import factorial
-from firedrake import (Interpolate, sqrt, inner, Function, SpatialCoordinate,
+from firedrake import (interpolate, sqrt, inner, Function, SpatialCoordinate,
                        FunctionSpace, VectorFunctionSpace, PointNotInDomainError,
                        Constant, assemble, dx)
 from firedrake.mesh import MeshGeometry
@@ -120,7 +120,7 @@ def triplot(mesh, axes=None, interior_kw={}, boundary_kw={}):
     if element.degree() != 1:
         # Interpolate to piecewise linear.
         V = VectorFunctionSpace(mesh, element.family(), 1)
-        coordinates = assemble(Interpolate(coordinates, V))
+        coordinates = assemble(interpolate(coordinates, V))
 
     coords = toreal(coordinates.dat.data_ro_with_halos, "real")
     result = []
@@ -215,7 +215,7 @@ def _plot_2d_field(method_name, function, *args, complex_component="real", **kwa
     if len(function.ufl_shape) == 1:
         element = function.ufl_element().sub_elements[0]
         Q = FunctionSpace(mesh, element)
-        function = assemble(Interpolate(sqrt(inner(function, function)), Q))
+        function = assemble(interpolate(sqrt(inner(function, function)), Q))
 
     num_sample_points = kwargs.pop("num_sample_points", 10)
     function_plotter = FunctionPlotter(mesh, num_sample_points)
@@ -326,7 +326,7 @@ def trisurf(function, *args, complex_component="real", **kwargs):
     if len(function.ufl_shape) == 1:
         element = function.ufl_element().sub_elements[0]
         Q = FunctionSpace(mesh, element)
-        function = assemble(Interpolate(sqrt(inner(function, function)), Q))
+        function = assemble(interpolate(sqrt(inner(function, function)), Q))
 
     num_sample_points = kwargs.pop("num_sample_points", 10)
     function_plotter = FunctionPlotter(mesh, num_sample_points)
@@ -355,7 +355,7 @@ def quiver(function, *, complex_component="real", **kwargs):
 
     coords = toreal(extract_unique_domain(function).coordinates.dat.data_ro, "real")
     V = extract_unique_domain(function).coordinates.function_space()
-    function_interp = assemble(Interpolate(function, V))
+    function_interp = assemble(interpolate(function, V))
     vals = toreal(function_interp.dat.data_ro, complex_component)
     C = np.linalg.norm(vals, axis=1)
     return axes.quiver(*(coords.T), *(vals.T), C, **kwargs)
@@ -816,7 +816,7 @@ def _bezier_plot(function, axes, complex_component="real", **kwargs):
     mesh = function.function_space().mesh()
     if deg == 0:
         V = FunctionSpace(mesh, "DG", 1)
-        interp = assemble(Interpolate(function, V))
+        interp = assemble(interpolate(function, V))
         return _bezier_plot(interp, axes, complex_component=complex_component,
                             **kwargs)
     y_vals = _bezier_calculate_points(function)
diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py
index 6540611752..11bab9fc0f 100644
--- a/firedrake/utility_meshes.py
+++ b/firedrake/utility_meshes.py
@@ -11,7 +11,7 @@
     Function,
     Constant,
     assemble,
-    Interpolate,
+    interpolate,
     FiniteElement,
     interval,
     tetrahedron,
@@ -2351,7 +2351,7 @@ def OctahedralSphereMesh(
     )
     if degree > 1:
         # use it to build a higher-order mesh
-        m = assemble(Interpolate(ufl.SpatialCoordinate(m), VectorFunctionSpace(m, "CG", degree)))
+        m = assemble(interpolate(ufl.SpatialCoordinate(m), VectorFunctionSpace(m, "CG", degree)))
         m = mesh.Mesh(
             m,
             name=name,
@@ -2386,11 +2386,11 @@ def OctahedralSphereMesh(
     # Make a copy of the coordinates so that we can blend two different
     # mappings near the pole
     Vc = m.coordinates.function_space()
-    Xlatitudinal = assemble(Interpolate(
+    Xlatitudinal = assemble(interpolate(
         Constant(radius) * ufl.as_vector([x * scale, y * scale, znew]), Vc
     ))
     Vlow = VectorFunctionSpace(m, "CG", 1)
-    Xlow = assemble(Interpolate(Xlatitudinal, Vlow))
+    Xlow = assemble(interpolate(Xlatitudinal, Vlow))
     r = ufl.sqrt(Xlow[0] ** 2 + Xlow[1] ** 2 + Xlow[2] ** 2)
     Xradial = Constant(radius) * Xlow / r
 
diff --git a/tests/firedrake/adjoint/test_reduced_functional.py b/tests/firedrake/adjoint/test_reduced_functional.py
index acd07d259f..d4f9a863d0 100644
--- a/tests/firedrake/adjoint/test_reduced_functional.py
+++ b/tests/firedrake/adjoint/test_reduced_functional.py
@@ -214,7 +214,7 @@ def test_interpolate():
 
     f = Function(V)
     f.dat.data[:] = 2
-    J = assemble(Interpolate(f**2, c))
+    J = assemble(interpolate(f**2, c))
     Jhat = ReducedFunctional(J, Control(f))
 
     h = Function(V)
@@ -244,7 +244,7 @@ def test_interpolate_mixed():
     f1, f2 = split(f)
     exprs = [f2 * div(f1)**2, grad(f2) * div(f1)]
     expr = as_vector([e[i] for e in exprs for i in np.ndindex(e.ufl_shape)])
-    J = assemble(Interpolate(expr, c))
+    J = assemble(interpolate(expr, c))
     Jhat = ReducedFunctional(J, Control(f))
 
     h = Function(V)
diff --git a/tests/firedrake/external_operators/test_external_operators.py b/tests/firedrake/external_operators/test_external_operators.py
index b6153f3d1f..a47953a566 100644
--- a/tests/firedrake/external_operators/test_external_operators.py
+++ b/tests/firedrake/external_operators/test_external_operators.py
@@ -104,7 +104,7 @@ def test_assemble(V, f):
     assert isinstance(jac, MatrixBase)
 
     # Assemble the exact Jacobian, i.e. the interpolation matrix: `Interpolate(dexpr(u,v,w)/du, V)`
-    jac_exact = assemble(Interpolate(derivative(expr(u, v, w), u), V))
+    jac_exact = assemble(interpolate(derivative(expr(u, v, w), u), V))
     np.allclose(jac.petscmat[:, :], jac_exact.petscmat[:, :], rtol=1e-14)
 
     # -- dNdu(u, v, w; δu, v*) (TLM) -- #
diff --git a/tests/firedrake/multigrid/test_poisson_gtmg.py b/tests/firedrake/multigrid/test_poisson_gtmg.py
index f70e5c6825..a4154dd392 100644
--- a/tests/firedrake/multigrid/test_poisson_gtmg.py
+++ b/tests/firedrake/multigrid/test_poisson_gtmg.py
@@ -60,7 +60,7 @@ def p1_callback():
     if custom_transfer:
         P1 = get_p1_space()
         V = FunctionSpace(mesh, "DGT", degree - 1)
-        I = assemble(Interpolate(TrialFunction(P1), V)).petscmat
+        I = assemble(interpolate(TrialFunction(P1), V)).petscmat
         R = PETSc.Mat().createTranspose(I)
         appctx['interpolation_matrix'] = I
         appctx['restriction_matrix'] = R
diff --git a/tests/firedrake/regression/test_adjoint_operators.py b/tests/firedrake/regression/test_adjoint_operators.py
index 03557bf435..d0212636c4 100644
--- a/tests/firedrake/regression/test_adjoint_operators.py
+++ b/tests/firedrake/regression/test_adjoint_operators.py
@@ -729,7 +729,7 @@ def test_copy_function():
     g = f.copy(deepcopy=True)
     J = assemble(g*dx)
     rf = ReducedFunctional(J, Control(f))
-    a = assemble(Interpolate(-one, V))
+    a = assemble(interpolate(-one, V))
     assert np.isclose(rf(a), -J)
 
 
diff --git a/tests/firedrake/regression/test_function.py b/tests/firedrake/regression/test_function.py
index 85c6cdc0a0..1f9ed7a429 100644
--- a/tests/firedrake/regression/test_function.py
+++ b/tests/firedrake/regression/test_function.py
@@ -81,22 +81,22 @@ def test_firedrake_tensor_function_nonstandard_shape(W_nonstandard_shape):
 
 def test_mismatching_rank_interpolation(V):
     f = Function(V)
-    with pytest.raises(RuntimeError):
+    with pytest.raises(ValueError):
         f.interpolate(Constant((1, 2)))
     VV = VectorFunctionSpace(V.mesh(), 'CG', 1)
     f = Function(VV)
-    with pytest.raises(RuntimeError):
+    with pytest.raises(ValueError):
         f.interpolate(Constant((1, 2)))
     VVV = TensorFunctionSpace(V.mesh(), 'CG', 1)
     f = Function(VVV)
-    with pytest.raises(RuntimeError):
+    with pytest.raises(ValueError):
         f.interpolate(Constant((1, 2)))
 
 
 def test_mismatching_shape_interpolation(V):
     VV = VectorFunctionSpace(V.mesh(), 'CG', 1)
     f = Function(VV)
-    with pytest.raises(RuntimeError):
+    with pytest.raises(ValueError):
         f.interpolate(Constant([1] * (VV.value_shape[0] + 1)))
 
 
diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py
index 50e29b05cb..ccd4de13f0 100644
--- a/tests/firedrake/regression/test_interp_dual.py
+++ b/tests/firedrake/regression/test_interp_dual.py
@@ -54,7 +54,7 @@ def test_assemble_interp_adjoint_tensor(mesh, V1, f1):
 
 def test_assemble_interp_operator(V2, f1):
     # Check type
-    If1 = Interpolate(f1, V2)
+    If1 = Interpolate(f1, Argument(V2.dual(), 0))
     assert isinstance(If1, ufl.Interpolate)
 
     # -- I(f1, V2) -- #
@@ -89,7 +89,7 @@ def test_assemble_interp_matrix(V1, V2, f1):
 def test_assemble_interp_tlm(V1, V2, f1):
     # -- Action(I(v1, V2), f1) -- #
     v1 = TrialFunction(V1)
-    Iv1 = Interpolate(v1, V2)
+    Iv1 = Interpolate(v1, Argument(V2.dual(), 0))
     b = assemble(interpolate(f1, V2))
 
     assembled_action_Iv1 = assemble(action(Iv1, f1))
@@ -99,7 +99,7 @@ def test_assemble_interp_tlm(V1, V2, f1):
 def test_assemble_interp_adjoint_matrix(V1, V2):
     # -- Adjoint(I(v1, V2)) -- #
     v1 = TrialFunction(V1)
-    Iv1 = Interpolate(v1, V2)
+    Iv1 = Interpolate(v1, Argument(V2.dual(), 0))
 
     v2 = TestFunction(V2)
     c2 = assemble(conj(v2) * dx)
@@ -120,7 +120,7 @@ def test_assemble_interp_adjoint_matrix(V1, V2):
 def test_assemble_interp_adjoint_model(V1, V2):
     # -- Action(Adjoint(I(v1, v2)), fstar) -- #
     v1 = TrialFunction(V1)
-    Iv1 = Interpolate(v1, V2)
+    Iv1 = Interpolate(v1, Argument(V2.dual(), 0))
 
     fstar = Cofunction(V2.dual())
     v = Argument(V1, 0)
@@ -167,9 +167,9 @@ def test_assemble_base_form_operator_expressions(mesh):
     f2 = Function(V1).interpolate(sin(2*pi*y))
     f3 = Function(V1).interpolate(cos(2*pi*x))
 
-    If1 = Interpolate(f1, V2)
-    If2 = Interpolate(f2, V2)
-    If3 = Interpolate(f3, V2)
+    If1 = Interpolate(f1, Argument(V2.dual(), 0))
+    If2 = Interpolate(f2, Argument(V2.dual(), 0))
+    If3 = Interpolate(f3, Argument(V2.dual(), 0))
 
     # Sum of BaseFormOperators (1-form)
     res = assemble(If1 + If2 + If3)
@@ -234,7 +234,7 @@ def test_solve_interp_f(mesh):
 
     # -- Solution where the source term is interpolated via `ufl.Interpolate`
     u2 = Function(V1)
-    If = Interpolate(f1, V2)
+    If = Interpolate(f1, Argument(V2.dual(), 0))
     # This requires assembling If
     F2 = inner(grad(u2), grad(w))*dx + inner(u2, w)*dx - inner(If, w)*dx
     solve(F2 == 0, u2)
@@ -267,7 +267,7 @@ def test_solve_interp_u(mesh):
     # -- Solution where u2 is interpolated via `ufl.Interpolate` (mat-free)
     u2 = Function(V1)
     # Iu is the identity
-    Iu = Interpolate(u2, V1)
+    Iu = Interpolate(u2, Argument(V1.dual(), 0))
     # This requires assembling the action the Jacobian of Iu
     F2 = inner(grad(u2), grad(w))*dx + inner(Iu, w)*dx - inner(f, w)*dx
     solve(F2 == 0, u2, solver_parameters={"mat_type": "matfree",
@@ -278,7 +278,7 @@ def test_solve_interp_u(mesh):
     # Same problem with grad(Iu) instead of grad(Iu)
     u2 = Function(V1)
     # Iu is the identity
-    Iu = Interpolate(u2, V1)
+    Iu = Interpolate(u2, Argument(V1.dual(), 0))
     # This requires assembling the action the Jacobian of Iu
     F2 = inner(grad(Iu), grad(w))*dx + inner(Iu, w)*dx - inner(f, w)*dx
     solve(F2 == 0, u2, solver_parameters={"mat_type": "matfree",
diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py
index b2054843cd..875d0fb0c1 100644
--- a/tests/firedrake/regression/test_interpolate_zany.py
+++ b/tests/firedrake/regression/test_interpolate_zany.py
@@ -117,7 +117,7 @@ def test_interpolate_zany_into_vom(V, mesh, which, expr_at_vom):
     P0 = expr_at_vom.function_space()
 
     # Interpolate a Function into P0(vom)
-    f_at_vom = assemble(Interpolate(fexpr, P0))
+    f_at_vom = assemble(interpolate(fexpr, P0))
     assert numpy.allclose(f_at_vom.dat.data_ro, expr_at_vom.dat.data_ro)
 
     # Construct a Cofunction on P0(vom)*
@@ -125,10 +125,10 @@ def test_interpolate_zany_into_vom(V, mesh, which, expr_at_vom):
     expected_action = assemble(action(Fvom, expr_at_vom))
 
     # Interpolate a Function into Fvom
-    f_at_vom = assemble(Interpolate(fexpr, Fvom))
+    f_at_vom = assemble(interpolate(fexpr, Fvom))
     assert numpy.allclose(f_at_vom, expected_action)
 
     # Interpolate a TestFunction into Fvom
-    expr_vom = assemble(Interpolate(vexpr, Fvom))
+    expr_vom = assemble(interpolate(vexpr, Fvom))
     f_at_vom = assemble(action(expr_vom, f))
     assert numpy.allclose(f_at_vom, expected_action)
diff --git a/tests/firedrake/submesh/test_submesh_interpolate.py b/tests/firedrake/submesh/test_submesh_interpolate.py
index a26c1acb08..92e422cf6a 100644
--- a/tests/firedrake/submesh/test_submesh_interpolate.py
+++ b/tests/firedrake/submesh/test_submesh_interpolate.py
@@ -51,7 +51,7 @@ def _test_submesh_interpolate_cell_cell(mesh, subdomain_cond, fe_fesub):
     f = Function(V_).interpolate(f)
     v0 = Coargument(V.dual(), 0)
     v1 = TrialFunction(Vsub)
-    interp = Interpolate(v1, v0, allow_missing_dofs=True)
+    interp = interpolate(v1, v0, allow_missing_dofs=True)
     A = assemble(interp)
     g = assemble(action(A, gsub))
     assert assemble(inner(g - f, g - f) * dx(label_value)).real < 1e-14
@@ -165,7 +165,7 @@ def test_submesh_interpolate_subcell_subcell_2_processes():
     f_l.dat.data_with_halos[:] = 3.0
     v0 = Coargument(V_r.dual(), 0)
     v1 = TrialFunction(V_l)
-    interp = Interpolate(v1, v0, allow_missing_dofs=True)
+    interp = interpolate(v1, v0, allow_missing_dofs=True)
     A = assemble(interp)
     f_r = assemble(action(A, f_l))
     g_r = Function(V_r).interpolate(conditional(x < 2.001, 3.0, 0.0))