diff --git a/demos/full_waveform_inversion/full_waveform_inversion.py.rst b/demos/full_waveform_inversion/full_waveform_inversion.py.rst
index fcc2ffaadc..147cbc1860 100644
--- a/demos/full_waveform_inversion/full_waveform_inversion.py.rst
+++ b/demos/full_waveform_inversion/full_waveform_inversion.py.rst
@@ -109,10 +109,12 @@ built over the ``my_ensemble.comm`` (spatial) communicator.
         dt = 0.03  # time step in seconds
         final_time = 0.6  # final time in seconds
         nx, ny = 15, 15
+        ftol = 0.9  # optimisation tolerance
     else:
         dt = 0.002  # time step in seconds
         final_time = 1.0  # final time in seconds
         nx, ny = 80, 80
+        ftol = 1e-2  # optimisation tolerance
 
     mesh = UnitSquareMesh(nx, ny, comm=my_ensemble.comm)
 
@@ -278,7 +280,9 @@ To have the step 4, we need first to tape the forward problem. That is done by c
 
 We now instantiate :class:`~.EnsembleReducedFunctional`::
 
-    J_hat = EnsembleReducedFunctional(J_val, Control(c_guess), my_ensemble)
+    J_hat = EnsembleReducedFunctional(J_val,
+                                      Control(c_guess, riesz_map="l2"),
+                                      my_ensemble)
 
 which enables us to recompute :math:`J` and its gradient :math:`\nabla_{\mathtt{c\_guess}} J`,
 where the :math:`J_s` and its gradients :math:`\nabla_{\mathtt{c\_guess}} J_s` are computed in parallel
@@ -286,13 +290,18 @@ based on the ``my_ensemble`` configuration.
 
 
 **Steps 4-6**: The instance of the :class:`~.EnsembleReducedFunctional`, named ``J_hat``,
-is then passed as an argument to the ``minimize`` function::
+is then passed as an argument to the ``minimize`` function. The default ``minimize`` function
+uses ``scipy.minimize``, and wraps the ``ReducedFunctional`` in a ``ReducedFunctionalNumPy``
+that handles transferring data between Firedrake and numpy data structures. However, because
+we have a custom ``ReducedFunctional``, we need to do this ourselves::
 
-    c_optimised = minimize(J_hat, method="L-BFGS-B", options={"disp": True, "maxiter": 1},
-                            bounds=(1.5, 2.0), derivative_options={"riesz_representation": 'l2'}
-                            )
+    from pyadjoint.reduced_functional_numpy import ReducedFunctionalNumPy
+    Jnumpy = ReducedFunctionalNumPy(J_hat)
 
-The ``minimize`` function executes the optimisation algorithm until the stopping criterion (``maxiter``) is met.
+    c_optimised = minimize(Jnumpy, method="L-BFGS-B", options={"disp": True, "ftol": ftol},
+                           bounds=(1.5, 2.0))
+
+The ``minimize`` function executes the optimisation algorithm until the stopping criterion (``ftol``) is met.
 For 20 iterations, the predicted velocity model is shown in the following figure.
 
 .. image:: c_predicted.png
@@ -303,9 +312,7 @@ For 20 iterations, the predicted velocity model is shown in the following figure
 .. warning::
 
     The ``minimize`` function uses the SciPy library for optimisation. However, for scenarios that require higher
-    levels of spatial parallelism, you should assess whether SciPy is the most suitable option for your problem.
-    SciPy's optimisation algorithm is not inner-product-aware. Therefore, we configure the options with
-    ``derivative_options={"riesz_representation": 'l2'}`` to account for this requirement.
+    levels of spatial parallelism, you should assess whether SciPy is the most suitable option for your problem such as the pyadjoint's TAOSolver.
 
 .. note::
 
diff --git a/demos/shape_optimization/shape_optimization.py.rst b/demos/shape_optimization/shape_optimization.py.rst
index 69c0776451..ee30704900 100644
--- a/demos/shape_optimization/shape_optimization.py.rst
+++ b/demos/shape_optimization/shape_optimization.py.rst
@@ -70,9 +70,10 @@ and evaluate the objective function::
 
 We now turn the objective function into a reduced function so that pyadjoint
 (and UFL shape differentiation capability) can automatically compute shape
-gradients, that is, directions of steepest ascent::
+gradients, that is, directions of steepest ascent. We also set the relevant
+Riesz map for this problem::
 
-  Jred = ReducedFunctional(J, Control(dT))
+  Jred = ReducedFunctional(J, Control(dT, riesz_map="H1"))
   stop_annotating()
 
 We now have all the ingredients to implement a basic steepest descent shape
@@ -84,8 +85,7 @@ optimization algorithm with fixed step size.::
     File.write(mesh.coordinates)
 
     # compute the gradient (steepest ascent)
-    opts = {"riesz_representation": "H1"}
-    gradJ = Jred.derivative(options=opts)
+    gradJ = Jred.derivative(apply_riesz=True)
 
     # update domain
     dT -= 0.2*gradJ
diff --git a/firedrake/adjoint/__init__.py b/firedrake/adjoint/__init__.py
index c48b990420..d3d28e6129 100644
--- a/firedrake/adjoint/__init__.py
+++ b/firedrake/adjoint/__init__.py
@@ -29,7 +29,7 @@
 from firedrake.adjoint_utils import get_solve_blocks  # noqa F401
 
 from pyadjoint.verification import taylor_test, taylor_to_dict  # noqa F401
-from pyadjoint.drivers import compute_gradient, compute_hessian  # noqa F401
+from pyadjoint.drivers import compute_gradient, compute_derivative, compute_hessian  # noqa F401
 from pyadjoint.adjfloat import AdjFloat  # noqa F401
 from pyadjoint.control import Control  # noqa F401
 from pyadjoint import IPOPTSolver, ROLSolver, MinimizationProblem, \
diff --git a/firedrake/adjoint/ensemble_reduced_functional.py b/firedrake/adjoint/ensemble_reduced_functional.py
index 0abf6c67af..fc17f2b1ed 100644
--- a/firedrake/adjoint/ensemble_reduced_functional.py
+++ b/firedrake/adjoint/ensemble_reduced_functional.py
@@ -1,11 +1,12 @@
-from pyadjoint import ReducedFunctional
+from pyadjoint.reduced_functional import AbstractReducedFunctional, ReducedFunctional
 from pyadjoint.enlisting import Enlist
 from pyop2.mpi import MPI
 
-import firedrake
+from firedrake.function import Function
+from firedrake.cofunction import Cofunction
 
 
-class EnsembleReducedFunctional(ReducedFunctional):
+class EnsembleReducedFunctional(AbstractReducedFunctional):
     """Enable solving simultaneously reduced functionals in parallel.
 
     Consider a functional :math:`J` and its gradient :math:`\\dfrac{dJ}{dm}`,
@@ -34,7 +35,7 @@ class EnsembleReducedFunctional(ReducedFunctional):
 
     Parameters
     ----------
-    J : pyadjoint.OverloadedType
+    functional : pyadjoint.OverloadedType
         An instance of an OverloadedType, usually :class:`pyadjoint.AdjFloat`.
         This should be the functional that we want to reduce.
     control : pyadjoint.Control or list of pyadjoint.Control
@@ -86,28 +87,40 @@ class EnsembleReducedFunctional(ReducedFunctional):
     works, please refer to the `Firedrake manual
     <https://www.firedrakeproject.org/parallelism.html#ensemble-parallelism>`_.
     """
-    def __init__(self, J, control, ensemble, scatter_control=True,
-                 gather_functional=None, derivative_components=None,
-                 scale=1.0, tape=None, eval_cb_pre=lambda *args: None,
+    def __init__(self, functional, control, ensemble, scatter_control=True,
+                 gather_functional=None,
+                 derivative_components=None,
+                 scale=1.0, tape=None,
+                 eval_cb_pre=lambda *args: None,
                  eval_cb_post=lambda *args: None,
                  derivative_cb_pre=lambda controls: controls,
                  derivative_cb_post=lambda checkpoint, derivative_components, controls: derivative_components,
-                 hessian_cb_pre=lambda *args: None, hessian_cb_post=lambda *args: None):
-        super(EnsembleReducedFunctional, self).__init__(
-            J, control, derivative_components=derivative_components,
-            scale=scale, tape=tape, eval_cb_pre=eval_cb_pre,
-            eval_cb_post=eval_cb_post, derivative_cb_pre=derivative_cb_pre,
+                 hessian_cb_pre=lambda *args: None,
+                 hessian_cb_post=lambda *args: None):
+        self.local_reduced_functional = ReducedFunctional(
+            functional, control,
+            derivative_components=derivative_components,
+            scale=scale, tape=tape,
+            eval_cb_pre=eval_cb_pre,
+            eval_cb_post=eval_cb_post,
+            derivative_cb_pre=derivative_cb_pre,
             derivative_cb_post=derivative_cb_post,
-            hessian_cb_pre=hessian_cb_pre, hessian_cb_post=hessian_cb_post)
+            hessian_cb_pre=hessian_cb_pre,
+            hessian_cb_post=hessian_cb_post
+        )
 
         self.ensemble = ensemble
         self.scatter_control = scatter_control
         self.gather_functional = gather_functional
 
+    @property
+    def controls(self):
+        return self.local_reduced_functional.controls
+
     def _allgather_J(self, J):
         if isinstance(J, float):
             vals = self.ensemble.ensemble_comm.allgather(J)
-        elif isinstance(J, firedrake.Function):
+        elif isinstance(J, Function):
             #  allgather not implemented in ensemble.py
             vals = []
             for i in range(self.ensemble.ensemble_comm.size):
@@ -134,7 +147,7 @@ def __call__(self, values):
             The computed value. Typically of instance of :class:`pyadjoint.AdjFloat`.
 
         """
-        local_functional = super(EnsembleReducedFunctional, self).__call__(values)
+        local_functional = self.local_reduced_functional(values)
         ensemble_comm = self.ensemble.ensemble_comm
         if self.gather_functional:
             controls_g = self._allgather_J(local_functional)
@@ -142,22 +155,23 @@ def __call__(self, values):
         # if gather_functional is None then we do a sum
         elif isinstance(local_functional, float):
             total_functional = ensemble_comm.allreduce(sendobj=local_functional, op=MPI.SUM)
-        elif isinstance(local_functional, firedrake.Function):
+        elif isinstance(local_functional, Function):
             total_functional = type(local_functional)(local_functional.function_space())
             total_functional = self.ensemble.allreduce(local_functional, total_functional)
         else:
             raise NotImplementedError("This type of functional is not supported.")
         return total_functional
 
-    def derivative(self, adj_input=1.0, options=None):
+    def derivative(self, adj_input=1.0, apply_riesz=False):
         """Compute derivatives of a functional with respect to the control parameters.
 
         Parameters
         ----------
         adj_input : float
             The adjoint input.
-        options : dict
-            Additional options for the derivative computation.
+        apply_riesz: bool
+            If True, apply the Riesz map of each control in order to return
+            a primal gradient rather than a derivative in the dual space.
 
         Returns
         -------
@@ -171,29 +185,62 @@ def derivative(self, adj_input=1.0, options=None):
 
         if self.gather_functional:
             dJg_dmg = self.gather_functional.derivative(adj_input=adj_input,
-                                                        options=options)
+                                                        apply_riesz=False)
             i = self.ensemble.ensemble_comm.rank
             adj_input = dJg_dmg[i]
 
-        dJdm_local = super(EnsembleReducedFunctional, self).derivative(adj_input=adj_input, options=options)
+        dJdm_local = self.local_reduced_functional.derivative(adj_input=adj_input,
+                                                              apply_riesz=apply_riesz)
 
         if self.scatter_control:
             dJdm_local = Enlist(dJdm_local)
             dJdm_total = []
 
             for dJdm in dJdm_local:
-                if not isinstance(dJdm, (firedrake.Function, float)):
-                    raise NotImplementedError("This type of gradient is not supported.")
+                if not isinstance(dJdm, (Cofunction, Function, float)):
+                    raise NotImplementedError(
+                        f"Gradients of type {type(dJdm).__name__} are not supported.")
 
                 dJdm_total.append(
                     self.ensemble.allreduce(dJdm, type(dJdm)(dJdm.function_space()))
-                    if isinstance(dJdm, firedrake.Function)
+                    if isinstance(dJdm, (Cofunction, Function))
                     else self.ensemble.ensemble_comm.allreduce(sendobj=dJdm, op=MPI.SUM)
                 )
             return dJdm_local.delist(dJdm_total)
         return dJdm_local
 
-    def hessian(self, m_dot, options=None):
+    def tlm(self, m_dot):
+        """Return the action of the tangent linear model of the functional.
+
+        The tangent linear model is evaluated w.r.t. the control on a vector
+        m_dot, around the last supplied value of the control.
+
+        Parameters
+        ----------
+        m_dot : pyadjoint.OverloadedType
+            The direction in which to compute the action of the tangent linear model.
+
+        Returns
+        -------
+            pyadjoint.OverloadedType: The action of the tangent linear model in the
+            direction m_dot.  Should be an instance of the same type as the functional.
+        """
+        local_tlm = self.local_reduced_functional.tlm(m_dot)
+        ensemble_comm = self.ensemble.ensemble_comm
+        if self.gather_functional:
+            mdot_g = self._allgather_J(local_tlm)
+            total_tlm = self.gather_functional.tlm(mdot_g)
+        # if gather_functional is None then we do a sum
+        elif isinstance(local_tlm, float):
+            total_tlm = ensemble_comm.allreduce(sendobj=local_tlm, op=MPI.SUM)
+        elif isinstance(local_tlm, Function):
+            total_tlm = type(local_tlm)(local_tlm.function_space())
+            total_tlm = self.ensemble.allreduce(local_tlm, total_tlm)
+        else:
+            raise NotImplementedError("This type of functional is not supported.")
+        return total_tlm
+
+    def hessian(self, m_dot, hessian_input=None, evaluate_tlm=True, apply_riesz=False):
         """The Hessian is not yet implemented for ensemble reduced functional.
 
         Raises:
diff --git a/firedrake/adjoint_utils/blocks/solving.py b/firedrake/adjoint_utils/blocks/solving.py
index cda49fca7d..2c1fb4f876 100644
--- a/firedrake/adjoint_utils/blocks/solving.py
+++ b/firedrake/adjoint_utils/blocks/solving.py
@@ -86,8 +86,15 @@ def _init_solver_parameters(self, args, kwargs):
         self.assemble_kwargs = {}
 
     def __str__(self):
-        return "solve({} = {})".format(ufl2unicode(self.lhs),
-                                       ufl2unicode(self.rhs))
+        try:
+            lhs_string = ufl2unicode(self.lhs)
+        except AttributeError:
+            lhs_string = str(self.lhs)
+        try:
+            rhs_string = ufl2unicode(self.rhs)
+        except AttributeError:
+            rhs_string = str(self.rhs)
+        return "solve({} = {})".format(lhs_string, rhs_string)
 
     def _create_F_form(self):
         # Process the equation forms, replacing values with checkpoints,
@@ -742,7 +749,7 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
         c = block_variable.output
         c_rep = block_variable.saved_output
 
-        if isinstance(c, firedrake.Function):
+        if isinstance(c, (firedrake.Function, firedrake.Cofunction)):
             trial_function = firedrake.TrialFunction(c.function_space())
         elif isinstance(c, firedrake.Constant):
             mesh = F_form.ufl_domain()
@@ -779,7 +786,12 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
         replace_map[self.func] = self.get_outputs()[0].saved_output
         dFdm = replace(dFdm, replace_map)
 
-        dFdm = dFdm * adj_sol
+        if isinstance(dFdm, firedrake.Argument):
+            #  Corner case. Should be fixed more permanently upstream in UFL.
+            #  See: https://github.com/FEniCS/ufl/issues/395
+            dFdm = ufl.Action(dFdm, adj_sol)
+        else:
+            dFdm = dFdm * adj_sol
         dFdm = firedrake.assemble(dFdm, **self.assemble_kwargs)
 
         return dFdm
diff --git a/firedrake/adjoint_utils/constant.py b/firedrake/adjoint_utils/constant.py
index fed5978a6d..5f0c83f98b 100644
--- a/firedrake/adjoint_utils/constant.py
+++ b/firedrake/adjoint_utils/constant.py
@@ -2,7 +2,6 @@
 from pyadjoint.adjfloat import AdjFloat
 from pyadjoint.tape import get_working_tape, annotate_tape
 from pyadjoint.overloaded_type import OverloadedType, create_overloaded_object
-from pyadjoint.reduced_functional_numpy import gather
 
 from firedrake.functionspace import FunctionSpace
 from firedrake.adjoint_utils.blocks import ConstantAssignBlock
@@ -58,15 +57,8 @@ def wrapper(self, *args, **kwargs):
 
         return wrapper
 
-    def get_derivative(self, options={}):
-        return self._ad_convert_type(self.adj_value, options=options)
-
-    def _ad_convert_type(self, value, options={}):
-        if value is None:
-            # TODO: Should the default be 0 constant here or return just None?
-            return type(self)(numpy.zeros(self.ufl_shape))
-        value = gather(value)
-        return self._constant_from_values(value)
+    def _ad_init_zero(self, dual=False):
+        return type(self)(numpy.zeros(self.ufl_shape))
 
     def _ad_function_space(self, mesh):
         element = self.ufl_element()
diff --git a/firedrake/adjoint_utils/ensemble_function.py b/firedrake/adjoint_utils/ensemble_function.py
index ac5b9c3413..179d017285 100644
--- a/firedrake/adjoint_utils/ensemble_function.py
+++ b/firedrake/adjoint_utils/ensemble_function.py
@@ -55,6 +55,13 @@ def _ad_dot(self, other, options=None):
     def _ad_convert_riesz(self, value, options=None):
         raise NotImplementedError
 
+    def _ad_init_zero(self, dual=False):
+        from firedrake import EnsembleFunction, EnsembleCofunction
+        if dual:
+            return EnsembleCofunction(self.function_space().dual())
+        else:
+            return EnsembleFunction(self.function_space())
+
     def _ad_create_checkpoint(self):
         if disk_checkpointing():
             raise NotImplementedError(
diff --git a/firedrake/adjoint_utils/function.py b/firedrake/adjoint_utils/function.py
index 5b064545cb..db4732ed55 100644
--- a/firedrake/adjoint_utils/function.py
+++ b/firedrake/adjoint_utils/function.py
@@ -2,7 +2,7 @@
 import ufl
 from ufl.domain import extract_unique_domain
 from pyadjoint.overloaded_type import create_overloaded_object, FloatingType
-from pyadjoint.tape import annotate_tape, stop_annotating, get_working_tape, no_annotations
+from pyadjoint.tape import annotate_tape, stop_annotating, get_working_tape
 from firedrake.adjoint_utils.blocks import FunctionAssignBlock, ProjectBlock, SubfunctionBlock, FunctionMergeBlock, SupermeshProjectBlock
 import firedrake
 from .checkpointing import disk_checkpointing, CheckpointFunction, \
@@ -218,72 +218,15 @@ def _ad_create_checkpoint(self):
         else:
             return self.copy(deepcopy=True)
 
-    def _ad_convert_riesz(self, value, options=None):
-        from firedrake import Function, Cofunction
-
-        options = {} if options is None else options
-        riesz_representation = options.get("riesz_representation", "L2")
-        solver_options = options.get("solver_options", {})
-        V = options.get("function_space", self.function_space())
-        if value == 0.:
-            # In adjoint-based differentiation, value == 0. arises only when
-            # the functional is independent on the control variable.
-            return Function(V)
-
-        if not isinstance(value, (Cofunction, Function)):
-            raise TypeError("Expected a Cofunction or a Function")
-
-        if riesz_representation == "l2":
-            return Function(V, val=value.dat)
-
-        elif riesz_representation in ("L2", "H1"):
-            if not isinstance(value, Cofunction):
-                raise TypeError("Expected a Cofunction")
-
-            ret = Function(V)
-            a = self._define_riesz_map_form(riesz_representation, V)
-            firedrake.solve(a == value, ret, **solver_options)
-            return ret
-
-        elif callable(riesz_representation):
-            return riesz_representation(value)
-
-        else:
-            raise ValueError(
-                "Unknown Riesz representation %s" % riesz_representation)
+    def _ad_convert_riesz(self, value, riesz_map=None):
+        return value.riesz_representation(riesz_map=riesz_map or "L2")
 
-    def _define_riesz_map_form(self, riesz_representation, V):
-        from firedrake import TrialFunction, TestFunction
-
-        u = TrialFunction(V)
-        v = TestFunction(V)
-        if riesz_representation == "L2":
-            a = firedrake.inner(u, v)*firedrake.dx
-
-        elif riesz_representation == "H1":
-            a = firedrake.inner(u, v)*firedrake.dx \
-                + firedrake.inner(firedrake.grad(u), firedrake.grad(v))*firedrake.dx
-
-        else:
-            raise NotImplementedError(
-                "Unknown Riesz representation %s" % riesz_representation)
-        return a
-
-    @no_annotations
-    def _ad_convert_type(self, value, options=None):
-        # `_ad_convert_type` is not annotated, unlike `_ad_convert_riesz`
-        options = {} if options is None else options.copy()
-        options.setdefault("riesz_representation", "L2")
-        if options["riesz_representation"] is None:
-            if value == 0.:
-                # In adjoint-based differentiation, value == 0. arises only when
-                # the functional is independent on the control variable.
-                V = options.get("function_space", self.function_space())
-                return firedrake.Cofunction(V.dual())
-            else:
-                return value
+    def _ad_init_zero(self, dual=False):
+        from firedrake import Function, Cofunction
+        if dual:
+            return Cofunction(self.function_space().dual())
         else:
-            return self._ad_convert_riesz(value, options=options)
+            return Function(self.function_space())
 
     def _ad_restore_at_checkpoint(self, checkpoint):
         if isinstance(checkpoint, CheckpointBase):
@@ -292,9 +235,7 @@ def _ad_restore_at_checkpoint(self, checkpoint):
             return checkpoint
 
     def _ad_will_add_as_dependency(self):
-        """Method called when the object is added as a Block dependency.
-
-        """
+        """Method called when the object is added as a Block dependency."""
         with checkpoint_init_data():
             super()._ad_will_add_as_dependency()
 
@@ -302,7 +243,8 @@ def _ad_mul(self, other):
         from firedrake import Function
 
         r = Function(self.function_space())
-        # `self` can be a Cofunction in which case only left multiplication with a scalar is allowed.
+        # `self` can be a Cofunction in which case only left multiplication
+        # with a scalar is allowed.
         r.assign(other * self)
         return r
 
@@ -314,7 +256,10 @@ def _ad_add(self, other):
         return r
 
     def _ad_dot(self, other, options=None):
-        from firedrake import assemble
+        from firedrake import assemble, action, Cofunction
+
+        if isinstance(other, Cofunction):
+            return assemble(action(other, self))
 
         options = {} if options is None else options
         riesz_representation = options.get("riesz_representation", "L2")
@@ -404,3 +349,9 @@ def _ad_to_petsc(self, vec=None):
 
     def __deepcopy__(self, memodict={}):
         return self.copy(deepcopy=True)
+
+
+class CofunctionMixin(FunctionMixin):
+
+    def _ad_dot(self, other):
+        return firedrake.assemble(firedrake.action(self, other))
diff --git a/firedrake/assemble.py b/firedrake/assemble.py
index e5adb33a62..9ca4d17237 100644
--- a/firedrake/assemble.py
+++ b/firedrake/assemble.py
@@ -1249,7 +1249,7 @@ def _apply_bc(self, tensor, bc, u=None):
             else:
                 # The residual belongs to a mixed space that is dual on the boundary nodes
                 # and primal on the interior nodes. Therefore, this is a type-safe operation.
-                r = tensor.riesz_representation("l2")
+                r = firedrake.Function(tensor.function_space().dual(), val=tensor.dat)
                 bc.apply(r, u=u)
         elif isinstance(bc, EquationBCSplit):
             bc.zero(tensor)
diff --git a/firedrake/cofunction.py b/firedrake/cofunction.py
index 374f7ee695..a2dbb10321 100644
--- a/firedrake/cofunction.py
+++ b/firedrake/cofunction.py
@@ -1,3 +1,4 @@
+from functools import cached_property
 import numpy as np
 import ufl
 
@@ -9,16 +10,16 @@
 import firedrake.functionspaceimpl as functionspaceimpl
 from firedrake import utils, vector, ufl_expr
 from firedrake.utils import ScalarType
-from firedrake.adjoint_utils.function import FunctionMixin
+from firedrake.adjoint_utils.function import CofunctionMixin
 from firedrake.adjoint_utils.checkpointing import DelegatedFunctionCheckpoint
 from firedrake.adjoint_utils.blocks.function import CofunctionAssignBlock
 from firedrake.petsc import PETSc
 
 
-class Cofunction(ufl.Cofunction, FunctionMixin):
+class Cofunction(ufl.Cofunction, CofunctionMixin):
     r"""A :class:`Cofunction` represents a function on a dual space.
-    Like Functions, cofunctions are
-    represented as sums of basis functions:
+
+    Like Functions, cofunctions are represented as sums of basis functions:
 
     .. math::
 
@@ -34,7 +35,7 @@ class Cofunction(ufl.Cofunction, FunctionMixin):
     """
 
     @PETSc.Log.EventDecorator()
-    @FunctionMixin._ad_annotate_init
+    @CofunctionMixin._ad_annotate_init
     def __init__(self, function_space, val=None, name=None, dtype=ScalarType,
                  count=None):
         r"""
@@ -106,7 +107,7 @@ def _analyze_form_arguments(self):
         self._coefficients = (self,)
 
     @utils.cached_property
-    @FunctionMixin._ad_annotate_subfunctions
+    @CofunctionMixin._ad_annotate_subfunctions
     def subfunctions(self):
         r"""Extract any sub :class:`Cofunction`\s defined on the component spaces
         of this this :class:`Cofunction`'s :class:`.FunctionSpace`."""
@@ -228,39 +229,45 @@ def assign(self, expr, subset=None, expr_from_assemble=False):
             Assigner(self, expr, subset).assign()
         return self
 
-    def riesz_representation(self, riesz_map='L2', **solver_options):
-        """Return the Riesz representation of this :class:`Cofunction` with respect to the given Riesz map.
+    def riesz_representation(self, riesz_map='L2', *, bcs=None,
+                             solver_options=None,
+                             form_compiler_parameters=None):
+        """Return the Riesz representation of this :class:`Cofunction`.
 
-        Example: For a L2 Riesz map, the Riesz representation is obtained by solving
-        the linear system ``Mx = self``, where M is the L2 mass matrix, i.e. M = <u, v>
-        with u and v trial and test functions, respectively.
+        Example: For a L2 Riesz map, the Riesz representation is obtained by
+        solving the linear system ``Mx = self``, where M is the L2 mass matrix,
+        i.e. M = <u, v> with u and v trial and test functions, respectively.
 
         Parameters
         ----------
-        riesz_map : str or collections.abc.Callable
-                    The Riesz map to use (`l2`, `L2`, or `H1`). This can also be a callable.
-        solver_options : dict
-                         Solver options to pass to the linear solver:
-                            - solver_parameters: optional solver parameters.
-                            - nullspace: an optional :class:`.VectorSpaceBasis` (or :class:`.MixedVectorSpaceBasis`)
-                                         spanning the null space of the operator.
-                            - transpose_nullspace: as for the nullspace, but used to make the right hand side consistent.
-                            - near_nullspace: as for the nullspace, but used to add the near nullspace.
-                            - options_prefix: an optional prefix used to distinguish PETSc options.
-                                              If not provided a unique prefix will be created.
-                                              Use this option if you want to pass options to the solver from the command line
-                                              in addition to through the ``solver_parameters`` dict.
+        riesz_map : str or ufl.sobolevspace.SobolevSpace or
+        collections.abc.Callable
+            The Riesz map to use (`l2`, `L2`, or `H1`). This can also be a
+            callable.
+        bcs: DirichletBC or list of DirichletBC
+            Boundary conditions to apply to the Riesz map.
+        solver_options: dict
+            A dictionary of PETSc options to be passed to the solver.
+        form_compiler_parameters: dict
+            A dictionary of form compiler parameters to be passed to the
+            variational problem that solves for the Riesz map.
 
         Returns
         -------
         firedrake.function.Function
-            Riesz representation of this :class:`Cofunction` with respect to the given Riesz map.
+            Riesz representation of this :class:`Cofunction` with respect to
+            the given Riesz map.
         """
-        return self._ad_convert_riesz(self, options={"function_space": self.function_space().dual(),
-                                                     "riesz_representation": riesz_map,
-                                                     "solver_options": solver_options})
+        if not callable(riesz_map):
+            riesz_map = RieszMap(
+                self.function_space(), riesz_map, bcs=bcs,
+                solver_parameters=solver_options,
+                form_compiler_parameters=form_compiler_parameters
+            )
+
+        return riesz_map(self)
 
-    @FunctionMixin._ad_annotate_iadd
+    @CofunctionMixin._ad_annotate_iadd
     @utils.known_pyop2_safe
     def __iadd__(self, expr):
 
@@ -276,7 +283,7 @@ def __iadd__(self, expr):
         # Let Python hit `BaseForm.__add__` which relies on ufl.FormSum.
         return NotImplemented
 
-    @FunctionMixin._ad_annotate_isub
+    @CofunctionMixin._ad_annotate_isub
     @utils.known_pyop2_safe
     def __isub__(self, expr):
 
@@ -293,7 +300,7 @@ def __isub__(self, expr):
         # Let Python hit `BaseForm.__sub__` which relies on ufl.FormSum.
         return NotImplemented
 
-    @FunctionMixin._ad_annotate_imul
+    @CofunctionMixin._ad_annotate_imul
     def __imul__(self, expr):
 
         if np.isscalar(expr):
@@ -372,3 +379,133 @@ def __str__(self):
 
     def cell_node_map(self):
         return self.function_space().cell_node_map()
+
+
+class RieszMap:
+    """Return a map between dual and primal function spaces.
+
+    A `RieszMap` can be called on a `Cofunction` in the appropriate space to
+    yield the `Function` which is the Riesz representer under the given inner
+    product. Conversely, it can be called on a `Function` to apply the given
+    inner product and return a `Cofunction`.
+
+    Parameters
+    ----------
+    function_space_or_inner_product: FunctionSpace or ufl.Form
+        The space from which to map, or a bilinear form defining an inner
+        product.
+    sobolev_space: str or ufl.sobolevspace.SobolevSpace.
+        Used to determine the inner product.
+    bcs: DirichletBC or list of DirichletBC
+        Boundary conditions to apply to the Riesz map.
+    solver_parameters: dict
+        A dictionary of PETSc options to be passed to the solver.
+    form_compiler_parameters: dict
+        A dictionary of form compiler parameters to be passed to the
+        variational problem that solves for the Riesz map.
+    restrict: bool
+        If `True`, use restricted function spaces in the Riesz map solver.
+    """
+
+    def __init__(self, function_space_or_inner_product=None,
+                 sobolev_space=ufl.L2, *, bcs=None, solver_parameters=None,
+                 form_compiler_parameters=None, restrict=True,
+                 constant_jacobian=False):
+        if isinstance(function_space_or_inner_product, ufl.Form):
+            args = ufl.algorithms.extract_arguments(
+                function_space_or_inner_product
+            )
+            if len(args) != 2:
+                raise ValueError(f"inner_product has arity {len(args)}, "
+                                 "should be 2.")
+            function_space = args[0].function_space()
+            inner_product = function_space_or_inner_product
+        else:
+            function_space = function_space_or_inner_product
+            if hasattr(function_space, "function_space"):
+                function_space = function_space.function_space()
+            if ufl.duals.is_dual(function_space):
+                function_space = function_space.dual()
+
+            if str(sobolev_space) == "l2":
+                inner_product = "l2"
+            else:
+                from firedrake import TrialFunction, TestFunction
+                u = TrialFunction(function_space)
+                v = TestFunction(function_space)
+                inner_product = RieszMap._inner_product_form(
+                    sobolev_space, u, v
+                )
+
+        self._function_space = function_space
+        self._inner_product = inner_product
+        self._bcs = bcs
+        self._solver_parameters = solver_parameters or {}
+        self._form_compiler_parameters = form_compiler_parameters or {}
+        self._restrict = restrict
+        self._constant_jacobian = constant_jacobian
+
+    @staticmethod
+    def _inner_product_form(sobolev_space, u, v):
+        from firedrake import inner, dx, grad
+        inner_products = {
+            "L2": lambda u, v: inner(u, v)*dx,
+            "H1": lambda u, v: inner(u, v)*dx + inner(grad(u), grad(v))*dx
+        }
+        try:
+            return inner_products[str(sobolev_space)](u, v)
+        except KeyError:
+            raise ValueError("No inner product defined for Sobolev space "
+                             f"{sobolev_space}.")
+
+    @cached_property
+    def _solver(self):
+        from firedrake import (LinearVariationalSolver,
+                               LinearVariationalProblem, Function, Cofunction)
+        rhs = Cofunction(self._function_space.dual())
+        soln = Function(self._function_space)
+        lvp = LinearVariationalProblem(
+            self._inner_product, rhs, soln, bcs=self._bcs,
+            restrict=self._restrict,
+            constant_jacobian=self._constant_jacobian,
+            form_compiler_parameters=self._form_compiler_parameters)
+        solver = LinearVariationalSolver(
+            lvp, solver_parameters=self._solver_parameters
+        )
+        return solver.solve, rhs, soln
+
+    def __call__(self, value):
+        """Return the Riesz representer of a Function or Cofunction."""
+        from firedrake import Function, Cofunction
+
+        if ufl.duals.is_dual(value):
+            if value.function_space().dual() != self._function_space:
+                raise ValueError("Function space mismatch in RieszMap.")
+            output = Function(self._function_space)
+
+            if self._inner_product == "l2":
+                for o, c in zip(output.subfunctions, value.subfunctions):
+                    o.dat.data[:] = c.dat.data[:]
+            else:
+                solve, rhs, soln = self._solver
+                rhs.assign(value)
+                solve()
+                output = Function(self._function_space)
+                output.assign(soln)
+        elif ufl.duals.is_primal(value):
+            if value.function_space() != self._function_space:
+                raise ValueError("Function space mismatch in RieszMap.")
+
+            if self._inner_product == "l2":
+                output = Cofunction(self._function_space.dual())
+                for o, c in zip(output.subfunctions, value.subfunctions):
+                    o.dat.data[:] = c.dat.data[:]
+            else:
+                output = firedrake.assemble(
+                    firedrake.action(self._inner_product, value)
+                )
+        else:
+            raise ValueError(
+                f"Unable to ascertain if {value} is primal or dual."
+            )
+        return output
diff --git a/firedrake/function.py b/firedrake/function.py
index c4b4686065..a4a0fb1191 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -19,7 +19,7 @@
 from firedrake.utils import ScalarType, IntType, as_ctypes
 
 from firedrake import functionspaceimpl
-from firedrake.cofunction import Cofunction
+from firedrake.cofunction import Cofunction, RieszMap
 from firedrake import utils
 from firedrake import vector
 from firedrake.adjoint_utils import FunctionMixin
@@ -462,39 +462,29 @@ def assign(self, expr, subset=None):
         return self
 
     def riesz_representation(self, riesz_map='L2'):
-        """Return the Riesz representation of this :class:`Function` with respect to the given Riesz map.
+        """Return the Riesz representation of this :class:`Function`.
 
-        Example: For a L2 Riesz map, the Riesz representation is obtained by taking the action
-        of ``M`` on ``self``, where M is the L2 mass matrix, i.e. M = <u, v>
-        with u and v trial and test functions, respectively.
+        Example: For a L2 Riesz map, the Riesz representation is obtained by
+        taking the action of ``M`` on ``self``, where M is the L2 mass matrix,
+        i.e. M = <u, v> with u and v trial and test functions, respectively.
 
         Parameters
         ----------
-        riesz_map : str or collections.abc.Callable
-                    The Riesz map to use (`l2`, `L2`, or `H1`). This can also be a callable.
+        riesz_map : str or ufl.sobolevspace.SobolevSpace or
+        collections.abc.Callable
+            The Riesz map to use (`l2`, `L2`, or `H1`). This can also be a
+            callable which applies the Riesz map.
 
         Returns
         -------
         firedrake.cofunction.Cofunction
-            Riesz representation of this :class:`Function` with respect to the given Riesz map.
+            Riesz representation of this :class:`Function` with respect to the
+            given Riesz map.
         """
-        from firedrake.ufl_expr import action
-        from firedrake.assemble import assemble
+        if not callable(riesz_map):
+            riesz_map = RieszMap(self.function_space(), riesz_map)
 
-        V = self.function_space()
-        if riesz_map == "l2":
-            return Cofunction(V.dual(), val=self.dat)
-
-        elif riesz_map in ("L2", "H1"):
-            a = self._define_riesz_map_form(riesz_map, V)
-            return assemble(action(a, self))
-
-        elif callable(riesz_map):
-            return riesz_map(self)
-
-        else:
-            raise NotImplementedError(
-                "Unknown Riesz representation %s" % riesz_map)
+        return riesz_map(self)
 
     @FunctionMixin._ad_annotate_iadd
     def __iadd__(self, expr):
diff --git a/firedrake/ml/jax/fem_operator.py b/firedrake/ml/jax/fem_operator.py
index b9e1914846..24c14a96b3 100644
--- a/firedrake/ml/jax/fem_operator.py
+++ b/firedrake/ml/jax/fem_operator.py
@@ -92,7 +92,7 @@ def bwd(self, _, grad_output: "jax.Array") -> "jax.Array":
             adj_input = float(adj_input)
 
         # Compute adjoint model of `F`: delegated to pyadjoint.ReducedFunctional
-        adj_output = self.F.derivative(adj_input=adj_input, options={'riesz_representation': None})
+        adj_output = self.F.derivative(adj_input=adj_input, apply_riesz=False)
 
         # Tuplify adjoint output
         adj_output = (adj_output,) if not isinstance(adj_output, collections.abc.Sequence) else adj_output
diff --git a/firedrake/ml/pytorch/fem_operator.py b/firedrake/ml/pytorch/fem_operator.py
index f7956265d0..59a1482587 100644
--- a/firedrake/ml/pytorch/fem_operator.py
+++ b/firedrake/ml/pytorch/fem_operator.py
@@ -83,7 +83,7 @@ def backward(ctx, grad_output):
             adj_input = float(adj_input)
 
         # Compute adjoint model of `F`: delegated to pyadjoint.ReducedFunctional
-        adj_output = F.derivative(adj_input=adj_input, options={"riesz_representation": None})
+        adj_output = F.derivative(adj_input=adj_input)
 
         # Tuplify adjoint output
         adj_output = (adj_output,) if not isinstance(adj_output, collections.abc.Sequence) else adj_output
diff --git a/firedrake/ufl_expr.py b/firedrake/ufl_expr.py
index bbc8a0d933..e4fc89b175 100644
--- a/firedrake/ufl_expr.py
+++ b/firedrake/ufl_expr.py
@@ -45,6 +45,12 @@ def __init__(self, function_space, number, part=None):
                                        number, part=part)
         self._function_space = function_space
 
+    def arguments(self):
+        return (self,)
+
+    def coefficients(self):
+        return ()
+
     @utils.cached_property
     def cell_node_map(self):
         return self.function_space().cell_node_map
diff --git a/firedrake/variational_solver.py b/firedrake/variational_solver.py
index f6b13aacb3..0ece7d7c30 100644
--- a/firedrake/variational_solver.py
+++ b/firedrake/variational_solver.py
@@ -87,7 +87,7 @@ def __init__(self, F, u, bcs=None, J=None,
                 bcs = J.bcs
         if bcs and any(isinstance(bc, EquationBC) for bc in bcs):
             restrict = False
-        self.restrict = restrict
+        self.restrict = restrict and bcs
 
         if restrict and bcs:
             V_res = restricted_function_space(V, extract_subdomain_ids(bcs))
diff --git a/tests/firedrake/adjoint/test_assemble.py b/tests/firedrake/adjoint/test_assemble.py
index 3d99770901..40937ee9eb 100644
--- a/tests/firedrake/adjoint/test_assemble.py
+++ b/tests/firedrake/adjoint/test_assemble.py
@@ -38,7 +38,7 @@ def test_assemble_0_forms():
     # where stored as a Function instead of Vector()
     s = a1 + a2 + 2.0 * a3
     rf = ReducedFunctional(s, Control(u))
-    dJdm = rf.derivative()
+    dJdm = rf.derivative(apply_riesz=True)
     assert_allclose(dJdm.dat.data_ro, 1. + 2. * 4. + 6. * 16.)
 
 
@@ -58,7 +58,7 @@ def test_assemble_0_forms_mixed():
     s -= a3  # this is done deliberately to end up with an adj_input of 0.0 for the a3 AssembleBlock
     rf = ReducedFunctional(s, Control(u))
     # derivative is: (1+4*u)*dx - summing is equivalent to testing with 1
-    dJdm = rf.derivative()
+    dJdm = rf.derivative(apply_riesz=True)
     assert_allclose(dJdm.dat.data_ro[0], 1. + 4. * 7)
     assert_allclose(dJdm.dat.data_ro[1], 0.0)
 
diff --git a/tests/firedrake/adjoint/test_assignment.py b/tests/firedrake/adjoint/test_assignment.py
index 4bada8a383..4ce7c07f9c 100644
--- a/tests/firedrake/adjoint/test_assignment.py
+++ b/tests/firedrake/adjoint/test_assignment.py
@@ -179,11 +179,11 @@ def test_assign_with_constant():
     J = assemble(u ** 2 * dx)
 
     rf = ReducedFunctional(J, Control(c))
-    dJdc = rf.derivative()
+    dJdc = rf.derivative(apply_riesz=True)
     assert_approx_equal(float(dJdc), 10.)
 
     rf = ReducedFunctional(J, Control(d))
-    dJdd = rf.derivative()
+    dJdd = rf.derivative(apply_riesz=True)
     assert_approx_equal(float(dJdd), 228.)
 
 
diff --git a/tests/firedrake/adjoint/test_disk_checkpointing.py b/tests/firedrake/adjoint/test_disk_checkpointing.py
index f264e5cb11..24d8fbfbd9 100644
--- a/tests/firedrake/adjoint/test_disk_checkpointing.py
+++ b/tests/firedrake/adjoint/test_disk_checkpointing.py
@@ -72,7 +72,7 @@ def adjoint_example(fine, coarse=None):
 
     assert np.allclose(J, Jnew)
 
-    grad_Jnew = Jhat.derivative()
+    grad_Jnew = Jhat.derivative(apply_riesz=True)
 
     return Jnew, grad_Jnew
 
diff --git a/tests/firedrake/adjoint/test_ensemble_reduced_functional.py b/tests/firedrake/adjoint/test_ensemble_reduced_functional.py
index 6d5285a687..a9fe89d132 100644
--- a/tests/firedrake/adjoint/test_ensemble_reduced_functional.py
+++ b/tests/firedrake/adjoint/test_ensemble_reduced_functional.py
@@ -1,5 +1,6 @@
 from firedrake import *
 from firedrake.adjoint import *
+from pyadjoint.reduced_functional_numpy import ReducedFunctionalNumPy
 import pytest
 from numpy.testing import assert_allclose
 
@@ -32,18 +33,15 @@ def test_verification():
     J = assemble(x * x * dx(domain=mesh))
     rf = EnsembleReducedFunctional(J, Control(x), ensemble)
     ensemble_J = rf(x)
-    dJdm = rf.derivative()
     assert_allclose(ensemble_J, size, rtol=1e-12)
+    dJdm = rf.derivative()
     assert_allclose(dJdm.dat.data_ro, 2.0 * size, rtol=1e-12)
     assert taylor_test(rf, x, Function(R, val=0.1)) > 1.9
 
 
 @pytest.mark.parallel(nprocs=4)
 @pytest.mark.skipcomplex  # Taping for complex-valued 0-forms not yet done
-@pytest.mark.xfail(reason="Taylor's test fails because the inner product \
-                   between the perturbation and gradient is not allreduced \
-                   for `scatter_control=False`.")
-def test_verification_gather_functional_adjfloat():
+def test_verification_gather_functional_adjfloat_evaluation():
     ensemble = Ensemble(COMM_WORLD, 2)
     rank = ensemble.ensemble_comm.rank
     mesh = UnitSquareMesh(4, 4, comm=ensemble.comm)
@@ -61,7 +59,6 @@ def test_verification_gather_functional_adjfloat():
     dJdm = rf.derivative()
     assert_allclose(ensemble_J, 1.0**4+2.0**4, rtol=1e-12)
     assert_allclose(dJdm.dat.data_ro, 4*(rank+1)**3, rtol=1e-12)
-    assert taylor_test(rf, x, Function(R, val=0.1)) > 1.9
 
 
 @pytest.mark.parallel(nprocs=4)
@@ -69,7 +66,26 @@ def test_verification_gather_functional_adjfloat():
 @pytest.mark.xfail(reason="Taylor's test fails because the inner product \
                    between the perturbation and gradient is not allreduced \
                    for `scatter_control=False`.")
-def test_verification_gather_functional_Function():
+def test_verification_gather_functional_adjfloat_taylor():
+    ensemble = Ensemble(COMM_WORLD, 2)
+    rank = ensemble.ensemble_comm.rank
+    mesh = UnitSquareMesh(4, 4, comm=ensemble.comm)
+    R = FunctionSpace(mesh, "R", 0)
+    x = function.Function(R, val=rank+1)
+    J = assemble(x * x * dx(domain=mesh))
+    a = AdjFloat(1.0)
+    b = AdjFloat(1.0)
+    Jg_m = [Control(a), Control(b)]
+    Jg = ReducedFunctional(a**2 + b**2, Jg_m)
+    rf = EnsembleReducedFunctional(J, Control(x), ensemble,
+                                   scatter_control=False,
+                                   gather_functional=Jg)
+    assert taylor_test(rf, x, Function(R, val=0.1)) > 1.9
+
+
+@pytest.mark.parallel(nprocs=4)
+@pytest.mark.skipcomplex  # Taping for complex-valued 0-forms not yet done
+def test_verification_gather_functional_Function_evaluation():
     ensemble = Ensemble(COMM_WORLD, 2)
     rank = ensemble.ensemble_comm.rank
     mesh = UnitSquareMesh(4, 4, comm=ensemble.comm)
@@ -88,17 +104,39 @@ def test_verification_gather_functional_Function():
     dJdm = rf.derivative()
     assert_allclose(ensemble_J, 1.0**4+2.0**4, rtol=1e-12)
     assert_allclose(dJdm.dat.data_ro, 4*(rank+1)**3, rtol=1e-12)
+
+
+@pytest.mark.parallel(nprocs=4)
+@pytest.mark.skipcomplex  # Taping for complex-valued 0-forms not yet done
+@pytest.mark.xfail(reason="Taylor's test fails because the inner product \
+                   between the perturbation and gradient is not allreduced \
+                   for `scatter_control=False`.")
+def test_verification_gather_functional_Function_taylor():
+    ensemble = Ensemble(COMM_WORLD, 2)
+    rank = ensemble.ensemble_comm.rank
+    mesh = UnitSquareMesh(4, 4, comm=ensemble.comm)
+    R = FunctionSpace(mesh, "R", 0)
+    x = function.Function(R, val=rank+1)
+    J = Function(R).assign(x**2)
+    a = Function(R).assign(1.0)
+    b = Function(R).assign(1.0)
+    Jg_m = [Control(a), Control(b)]
+    Jg = assemble((a**2 + b**2)*dx)
+    Jghat = ReducedFunctional(Jg, Jg_m)
+    rf = EnsembleReducedFunctional(J, Control(x), ensemble,
+                                   scatter_control=False,
+                                   gather_functional=Jghat)
     assert taylor_test(rf, x, Function(R, val=0.1)) > 1.9
 
 
-@pytest.mark.parallel(nprocs=6)
+@pytest.mark.parallel(nprocs=3)
 @pytest.mark.skipcomplex  # Taping for complex-valued 0-forms not yet done
 def test_minimise():
     # Optimisation test using a list of controls.
     # This test is equivalent to the one found at:
-    # https://github.com/dolfin-adjoint/pyadjoint/blob/master/tests/firedrake_adjoint/test_optimisation.py#L9.
+    # https://github.com/firedrakeproject/firedrake/blob/master/tests/firedrake/adjoint/test_optimisation.py#L92
     # In this test, the functional is the result of an ensemble allreduce operation.
-    ensemble = Ensemble(COMM_WORLD, 2)
+    ensemble = Ensemble(COMM_WORLD, 1)
     mesh = UnitSquareMesh(4, 4, comm=ensemble.comm)
     R = FunctionSpace(mesh, "R", 0)
     n = 2
@@ -106,8 +144,10 @@ def test_minimise():
     c = [Control(xi) for xi in x]
     # Rosenbrock function https://en.wikipedia.org/wiki/Rosenbrock_function
     # with minimum at x = (1, 1, 1, ...)
-    f = 100*(x[1] - x[0]**2)**2 + (1 - x[0])**2
-    J = assemble(f * dx(domain=mesh))
-    rf = EnsembleReducedFunctional(J, c, ensemble)
-    result = minimize(rf)
+    with set_working_tape():
+        f = 100*(x[1] - x[0]**2)**2 + (1 - x[0])**2
+        J = assemble(f * dx(domain=mesh))
+        rf = EnsembleReducedFunctional(J, c, ensemble)
+    rf_np = ReducedFunctionalNumPy(rf)
+    result = minimize(rf_np)
     assert_allclose([float(xi) for xi in result], 1., rtol=1e-8)
diff --git a/tests/firedrake/adjoint/test_hessian.py b/tests/firedrake/adjoint/test_hessian.py
index 403cb8a86e..05c6a4bfad 100644
--- a/tests/firedrake/adjoint/test_hessian.py
+++ b/tests/firedrake/adjoint/test_hessian.py
@@ -60,8 +60,8 @@ def test_simple_solve():
     tape.evaluate_adj()
 
     m = f.copy(deepcopy=True)
-    dJdm = assemble(inner(Jhat.derivative(), h)*dx)
-    Hm = assemble(inner(Jhat.hessian(h), h)*dx)
+    dJdm = assemble(inner(Jhat.derivative(apply_riesz=True), h)*dx)
+    Hm = assemble(inner(Jhat.hessian(h, apply_riesz=True), h)*dx)
     assert taylor_test(Jhat, m, h, dJdm=dJdm, Hm=Hm) > 2.9
 
 
@@ -140,7 +140,7 @@ def test_function():
     J = assemble(c ** 2 * u ** 2 * dx)
 
     Jhat = ReducedFunctional(J, [control_c, control_f])
-    dJdc, dJdf = compute_gradient(J, [control_c, control_f])
+    dJdc, dJdf = compute_gradient(J, [control_c, control_f], apply_riesz=True)
 
     # Step direction for derivatives and convergence test
     h_c = Function(R, val=1.0)
@@ -148,11 +148,11 @@ def test_function():
     h_f.vector()[:] = 10*rng.random(V.dim())
 
     # Total derivative
-    dJdc, dJdf = compute_gradient(J, [control_c, control_f])
+    dJdc, dJdf = compute_gradient(J, [control_c, control_f], apply_riesz=True)
     dJdm = assemble(dJdc * h_c * dx + dJdf * h_f * dx)
 
     # Hessian
-    Hcc, Hff = compute_hessian(J, [control_c, control_f], [h_c, h_f])
+    Hcc, Hff = compute_hessian(J, [control_c, control_f], [h_c, h_f], apply_riesz=True)
     Hm = assemble(Hcc * h_c * dx + Hff * h_f * dx)
     assert taylor_test(Jhat, [c, f], [h_c, h_f], dJdm=dJdm, Hm=Hm) > 2.9
 
@@ -241,7 +241,8 @@ def test_dirichlet():
 
 
 @pytest.mark.skipcomplex
-def test_burgers():
+@pytest.mark.parametrize("solve_type", ["solve", "nlvs"])
+def test_burgers(solve_type):
     tape = Tape()
     set_working_tape(tape)
     n = 100
@@ -273,9 +274,24 @@ def Dt(u, u_, timestep):
     F = (Dt(u, ic, timestep)*v
          + u*u.dx(0)*v + nu*u.dx(0)*v.dx(0))*dx
     bc = DirichletBC(V, 0.0, "on_boundary")
-
     t = 0.0
-    solve(F == 0, u, bc, solver_parameters=params)
+
+    if solve_type == "nlvs":
+        use_nlvs = True
+    elif solve_type == "solve":
+        use_nlvs = False
+    else:
+        raise ValueError(f"Unrecognised solve type {solve_type}")
+
+    if use_nlvs:
+        solver = NonlinearVariationalSolver(
+            NonlinearVariationalProblem(F, u),
+            solver_parameters=params)
+
+    if use_nlvs:
+        solver.solve()
+    else:
+        solve(F == 0, u, bc, solver_parameters=params)
     u_.assign(u)
     t += float(timestep)
 
@@ -284,7 +300,10 @@ def Dt(u, u_, timestep):
 
     end = 0.2
     while (t <= end):
-        solve(F == 0, u, bc, solver_parameters=params)
+        if use_nlvs:
+            solver.solve()
+        else:
+            solve(F == 0, u, bc, solver_parameters=params)
         u_.assign(u)
 
         t += float(timestep)
diff --git a/tests/firedrake/adjoint/test_optimisation.py b/tests/firedrake/adjoint/test_optimisation.py
index c4521c0ee7..2189a7882b 100644
--- a/tests/firedrake/adjoint/test_optimisation.py
+++ b/tests/firedrake/adjoint/test_optimisation.py
@@ -66,27 +66,27 @@ def test_petsc_roundtrip_multiple():
     assert (u_2.dat.data_ro == u_2_test.dat.data_ro).all()
 
 
-def minimize_tao_lmvm(rf, *, convert_options=None):
+def minimize_tao_lmvm(rf):
     problem = MinimizationProblem(rf)
     solver = TAOSolver(problem, {"tao_type": "lmvm",
                                  "tao_monitor": None,
                                  "tao_converged_reason": None,
                                  "tao_gatol": 1.0e-5,
                                  "tao_grtol": 0.0,
-                                 "tao_gttol": 0.0},
-                       convert_options=convert_options)
+                                 "tao_gttol": 1.0e-6,
+                                 "tao_monitor": None})
     return solver.solve()
 
 
-def minimize_tao_nls(rf, *, convert_options=None):
+def minimize_tao_nls(rf):
     problem = MinimizationProblem(rf)
     solver = TAOSolver(problem, {"tao_type": "nls",
                                  "tao_monitor": None,
                                  "tao_converged_reason": None,
                                  "tao_gatol": 1.0e-5,
                                  "tao_grtol": 0.0,
-                                 "tao_gttol": 0.0},
-                       convert_options=convert_options)
+                                 "tao_gttol": 1.0e-6,
+                                 "tao_monitor": None})
     return solver.solve()
 
 
@@ -121,8 +121,13 @@ def _simple_helmholz_model(V, source):
     return u
 
 
+@pytest.mark.parametrize(
+    "riesz_representation",
+    [None,
+     "l2",
+     pytest.param("H1", marks=pytest.mark.xfail(reason="H1 is the wrong norm for this problem"))])
 @pytest.mark.skipcomplex
-def test_simple_inversion():
+def test_simple_inversion(riesz_representation):
     """Test inversion of source term in helmholze eqn."""
     mesh = UnitIntervalMesh(10)
     V = FunctionSpace(mesh, "CG", 1)
@@ -136,7 +141,7 @@ def test_simple_inversion():
 
     # now rerun annotated model with zero source
     source = Function(V)
-    c = Control(source)
+    c = Control(source, riesz_map=riesz_representation)
     u = _simple_helmholz_model(V, source)
 
     J = assemble(1e6 * (u - u_ref)**2*dx)
@@ -144,13 +149,6 @@ def test_simple_inversion():
 
     x = minimize(rf)
     assert_allclose(x.dat.data, source_ref.dat.data, rtol=1e-2)
-    rf(source)
-    x = minimize(rf, derivative_options={"riesz_representation": "l2"})
-    assert_allclose(x.dat.data, source_ref.dat.data, rtol=1e-2)
-    rf(source)
-    x = minimize(rf, derivative_options={"riesz_representation": "H1"})
-    # Assert that the optimisation does not converge for H1 representation
-    assert not np.allclose(x.dat.data, source_ref.dat.data, rtol=1e-2)
 
 
 @pytest.mark.parametrize("minimize", [minimize_tao_lmvm,
@@ -158,7 +156,7 @@ def test_simple_inversion():
 @pytest.mark.parametrize("riesz_representation", [None, "l2", "L2", "H1"])
 @pytest.mark.skipcomplex
 def test_tao_simple_inversion(minimize, riesz_representation):
-    """Test inversion of source term in helmholze eqn using TAO."""
+    """Test inversion of source term in helmholtz eqn using TAO."""
     mesh = UnitIntervalMesh(10)
     V = FunctionSpace(mesh, "CG", 1)
     source_ref = Function(V)
@@ -171,14 +169,13 @@ def test_tao_simple_inversion(minimize, riesz_representation):
 
     # now rerun annotated model with zero source
     source = Function(V)
-    c = Control(source)
+    c = Control(source, riesz_map=riesz_representation)
     u = _simple_helmholz_model(V, source)
 
     J = assemble(1e6 * (u - u_ref)**2*dx)
     rf = ReducedFunctional(J, c)
 
-    x = minimize(rf, convert_options=(None if riesz_representation is None
-                                      else {"riesz_representation": riesz_representation}))
+    x = minimize(rf)
     assert_allclose(x.dat.data, source_ref.dat.data, rtol=1e-2)
 
 
@@ -304,7 +301,7 @@ def test_simple_inversion_riesz_representation(tao_type):
         u_ref = _simple_helmholz_model(V, source_ref)
 
     def forward(source):
-        c = Control(source)
+        c = Control(source, riesz_map=riesz_representation)
         u = _simple_helmholz_model(V, source)
 
         J = assemble(1e6 * (u - u_ref)**2*dx)
@@ -315,9 +312,7 @@ def forward(source):
     source = Function(V)
     rf = forward(source)
     with stop_annotating():
-        solver = TAOSolver(
-            MinimizationProblem(rf), tao_parameters,
-            convert_options={"riesz_representation": riesz_representation})
+        solver = TAOSolver(MinimizationProblem(rf), tao_parameters)
         x = solver.solve()
         assert_allclose(x.dat.data, source_ref.dat.data, rtol=1e-2)
 
@@ -327,7 +322,7 @@ def forward(source):
                                      mfn_parameters=mfn_parameters)
 
     def forward_transform(source):
-        c = Control(source)
+        c = Control(source, riesz_map="l2")
         source = transform(source, TransformType.PRIMAL,
                            riesz_representation,
                            mfn_parameters=mfn_parameters)
@@ -340,8 +335,7 @@ def forward_transform(source):
 
     with stop_annotating():
         solver_transform = TAOSolver(
-            MinimizationProblem(rf_transform), tao_parameters,
-            convert_options={"riesz_representation": "l2"})
+            MinimizationProblem(rf_transform), tao_parameters)
         x_transform = transform(solver_transform.solve(), TransformType.PRIMAL,
                                 riesz_representation,
                                 mfn_parameters=mfn_parameters)
diff --git a/tests/firedrake/adjoint/test_reduced_functional.py b/tests/firedrake/adjoint/test_reduced_functional.py
index edd483c040..acd07d259f 100644
--- a/tests/firedrake/adjoint/test_reduced_functional.py
+++ b/tests/firedrake/adjoint/test_reduced_functional.py
@@ -31,7 +31,7 @@ def test_constant():
 
     c = Function(R, val=1)
     f = Function(V)
-    f.vector()[:] = 1
+    f.assign(1)
 
     u = Function(V)
     v = TestFunction(V)
@@ -52,7 +52,7 @@ def test_function():
 
     c = Constant(1)
     f = Function(V)
-    f.vector()[:] = 1
+    f.assign(1)
 
     u = Function(V)
     v = TestFunction(V)
@@ -65,19 +65,18 @@ def test_function():
     Jhat = ReducedFunctional(J, Control(f))
 
     h = Function(V)
-    h.vector()[:] = rand(V.dof_dset.size)
+    h.dat.data[:] = rand(V.dof_dset.size)
     assert taylor_test(Jhat, f, h) > 1.9
 
 
 @pytest.mark.skipcomplex
-@pytest.mark.parametrize("control", ["dirichlet", "neumann"])
+@pytest.mark.parametrize("control", ["dirichlet-pre", "dirichlet-post", "neumann"])
 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)
+    u = Function(V)
     v = TestFunction(V)
 
     x, y = SpatialCoordinate(mesh)
@@ -89,20 +88,21 @@ def test_wrt_function_dirichlet_boundary(control):
     g1 = Function(R, val=2)
     g2 = Function(R, val=1)
     f = Function(V)
-    f.vector()[:] = 10
+    f.assign(10)
 
-    a = inner(grad(u), grad(v))*dx
+    a = inner(grad(u), grad(v))*dx + u**2*v*dx
     L = inner(f, v)*dx + inner(g1, v)*ds(4) + inner(g2, v)*ds(3)
 
-    solve(a == L, u_, bc)
+    pre_apply_bcs = (control == "dirichlet-post")
+    solve(a - L == 0, u, bc, pre_apply_bcs=pre_apply_bcs)
 
-    J = assemble(u_**2*dx)
+    J = assemble(u**2*dx)
 
-    if control == "dirichlet":
+    if control.startswith("dirichlet"):
         Jhat = ReducedFunctional(J, Control(bc_func))
         g = bc_func
         h = Function(V)
-        h.vector()[:] = 1
+        h.assign(1)
     else:
         Jhat = ReducedFunctional(J, Control(g1))
         g = g1
@@ -131,10 +131,10 @@ def test_time_dependent():
     T = 0.5
     dt = 0.1
     f = Function(V)
-    f.vector()[:] = 1
+    f.assign(1)
 
     u_1 = Function(V)
-    u_1.vector()[:] = 1
+    u_1.assign(1)
     control = Control(u_1)
 
     a = u_1*u*v*dx + dt*f*inner(grad(u), grad(v))*dx
@@ -152,7 +152,7 @@ def test_time_dependent():
     Jhat = ReducedFunctional(J, control)
 
     h = Function(V)
-    h.vector()[:] = 1
+    h.assign(1)
     assert taylor_test(Jhat, control.tape_value(), h) > 1.9
 
 
@@ -172,7 +172,7 @@ def test_mixed_boundary():
     g1 = Constant(2)
     g2 = Constant(1)
     f = Function(V)
-    f.vector()[:] = 10
+    f.assign(10)
 
     a = f*inner(grad(u), grad(v))*dx
     L = inner(f, v)*dx + inner(g1, v)*ds(4) + inner(g2, v)*ds(3)
@@ -183,7 +183,7 @@ def test_mixed_boundary():
 
     Jhat = ReducedFunctional(J, Control(f))
     h = Function(V)
-    h.vector()[:] = 1
+    h.assign(1)
     assert taylor_test(Jhat, f, h) > 1.9
 
 
@@ -194,13 +194,13 @@ def test_assemble_recompute():
     R = FunctionSpace(mesh, "R", 0)
 
     f = Function(V)
-    f.vector()[:] = 2
+    f.assign(2)
     expr = Function(R).assign(assemble(f**2*dx))
     J = assemble(expr**2*dx(domain=mesh))
     Jhat = ReducedFunctional(J, Control(f))
 
     h = Function(V)
-    h.vector()[:] = 1
+    h.assign(1)
     assert taylor_test(Jhat, f, h) > 1.9
 
 
diff --git a/tests/firedrake/adjoint/test_shape_derivatives.py b/tests/firedrake/adjoint/test_shape_derivatives.py
index cb23522ce5..d5950b1586 100644
--- a/tests/firedrake/adjoint/test_shape_derivatives.py
+++ b/tests/firedrake/adjoint/test_shape_derivatives.py
@@ -39,8 +39,7 @@ def test_sin_weak_spatial():
     V = TestFunction(S)
     # Derivative (Cofunction)
     dJV = assemble(div(V)*sin(x[0])*dx + V[0]*cos(x[0])*dx)
-    # Apply L2 riesz representation to obtain the gradient.
-    actual = dJV.riesz_representation().vector().get_local()
+    actual = dJV.vector().get_local()
     assert np.allclose(computed, actual, rtol=1e-14)
 
 
@@ -95,8 +94,8 @@ def test_shape_hessian():
     h.interpolate(as_vector((cos(x[2]), A*cos(x[1]), A*x[1])))
 
     # Second order taylor
-    dJdm = assemble(inner(Jhat.derivative(), h)*dx)
-    Hm = assemble(inner(compute_hessian(J, c, h), h)*dx)
+    dJdm = assemble(inner(Jhat.derivative(apply_riesz=True), h)*dx)
+    Hm = assemble(inner(compute_hessian(J, c, h, apply_riesz=True), h)*dx)
     r2 = taylor_test(Jhat, s, h, dJdm=dJdm, Hm=Hm)
     print(r2)
     assert (r2 > 2.9)
@@ -150,8 +149,8 @@ def test_PDE_hessian_neumann():
     Jhat(s)
 
     # # Second order taylor
-    dJdm = assemble(inner(Jhat.derivative(), h)*dx)
-    Hm = assemble(inner(compute_hessian(J, c, h), h)*dx)
+    dJdm = assemble(inner(Jhat.derivative(apply_riesz=True), h)*dx)
+    Hm = assemble(inner(compute_hessian(J, c, h, apply_riesz=True), h)*dx)
     r2 = taylor_test(Jhat, s, h, dJdm=dJdm, Hm=Hm)
     assert (r2 > 2.95)
 
@@ -203,8 +202,8 @@ def test_PDE_hessian_dirichlet():
     Jhat(s)
 
     # # Second order taylor
-    dJdm = assemble(inner(Jhat.derivative(), h)*dx)
-    Hm = assemble(inner(compute_hessian(J, c, h), h)*dx)
+    dJdm = assemble(inner(Jhat.derivative(apply_riesz=True), h)*dx)
+    Hm = assemble(inner(compute_hessian(J, c, h, apply_riesz=True), h)*dx)
     r2 = taylor_test(Jhat, s, h, dJdm=dJdm, Hm=Hm)
     assert (r2 > 2.95)
 
diff --git a/tests/firedrake/adjoint/test_solving.py b/tests/firedrake/adjoint/test_solving.py
index d6e5a162c9..bfbf7b348c 100644
--- a/tests/firedrake/adjoint/test_solving.py
+++ b/tests/firedrake/adjoint/test_solving.py
@@ -76,7 +76,8 @@ def J(f):
 
 
 @pytest.mark.skipcomplex
-def test_nonlinear_problem():
+@pytest.mark.parametrize("pre_apply_bcs", (True, False))
+def test_nonlinear_problem(pre_apply_bcs):
     """This tests whether nullspace and solver_parameters are passed on in adjoint solves"""
     mesh = IntervalMesh(10, 0, 1)
     V = FunctionSpace(mesh, "Lagrange", 1)
@@ -91,7 +92,7 @@ def test_nonlinear_problem():
     def J(f):
         a = f*inner(grad(u), grad(v))*dx + u**2*v*dx - f*v*dx
         L = 0
-        solve(a == L, u, bc)
+        solve(a == L, u, bc, pre_apply_bcs=pre_apply_bcs)
         return assemble(u**2*dx)
 
     _test_adjoint(J, f)
diff --git a/tests/firedrake/adjoint/test_tlm.py b/tests/firedrake/adjoint/test_tlm.py
index ff6fbe7a12..14148a233d 100644
--- a/tests/firedrake/adjoint/test_tlm.py
+++ b/tests/firedrake/adjoint/test_tlm.py
@@ -54,9 +54,7 @@ def test_tlm_assemble():
     h = Function(V)
     h.vector()[:] = rand(h.dof_dset.size)
     g = f.copy(deepcopy=True)
-    f.block_variable.tlm_value = h
-    tape.evaluate_tlm()
-    assert (taylor_test(Jhat, g, h, dJdm=J.block_variable.tlm_value) > 1.9)
+    assert (taylor_test(Jhat, g, h, dJdm=Jhat.tlm(h)) > 1.9)
 
 
 @pytest.mark.skipcomplex
@@ -80,12 +78,7 @@ def test_tlm_bc():
     J = assemble(c ** 2 * u * dx)
     Jhat = ReducedFunctional(J, Control(c))
 
-    # 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 = Function(R, val=1)
-    tape.evaluate_tlm()
-
-    assert (taylor_test(Jhat, c, Function(R, val=1), dJdm=J.block_variable.tlm_value) > 1.9)
+    assert (taylor_test(Jhat, c, Function(R, val=1), dJdm=Jhat.tlm(Function(R, val=1))) > 1.9)
 
 
 @pytest.mark.skipcomplex
@@ -113,10 +106,8 @@ def test_tlm_func():
     h = Function(V)
     h.vector()[:] = rand(h.dof_dset.size)
     g = c.copy(deepcopy=True)
-    c.block_variable.tlm_value = h
-    tape.evaluate_tlm()
 
-    assert (taylor_test(Jhat, g, h, dJdm=J.block_variable.tlm_value) > 1.9)
+    assert (taylor_test(Jhat, g, h, dJdm=Jhat.tlm(h)) > 1.9)
 
 
 @pytest.mark.parametrize("solve_type",
@@ -170,9 +161,7 @@ def test_time_dependent(solve_type):
     Jhat = ReducedFunctional(J, control)
     h = Function(V)
     h.vector()[:] = rand(h.dof_dset.size)
-    u_1.tlm_value = h
-    tape.evaluate_tlm()
-    assert (taylor_test(Jhat, control.tape_value(), h, dJdm=J.block_variable.tlm_value) > 1.9)
+    assert (taylor_test(Jhat, control.tape_value(), h, dJdm=Jhat.tlm(h)) > 1.9)
 
 
 @pytest.mark.skipcomplex
@@ -224,9 +213,7 @@ def Dt(u, u_, timestep):
     h = Function(V)
     h.vector()[:] = rand(h.dof_dset.size)
     g = ic.copy(deepcopy=True)
-    ic.block_variable.tlm_value = h
-    tape.evaluate_tlm()
-    assert (taylor_test(Jhat, g, h, dJdm=J.block_variable.tlm_value) > 1.9)
+    assert (taylor_test(Jhat, g, h, dJdm=Jhat.tlm(h)) > 1.9)
 
 
 @pytest.mark.skipcomplex
@@ -255,9 +242,7 @@ def test_projection():
     J = assemble(u_**2*dx)
     Jhat = ReducedFunctional(J, Control(k))
 
-    k.block_variable.tlm_value = Constant(1)
-    tape.evaluate_tlm()
-    assert (taylor_test(Jhat, k, Function(R, val=1), dJdm=J.block_variable.tlm_value) > 1.9)
+    assert (taylor_test(Jhat, k, Function(R, val=1), dJdm=Jhat.tlm(Constant(1))) > 1.9)
 
 
 @pytest.mark.skipcomplex
@@ -268,7 +253,6 @@ def test_projection_function():
     V = FunctionSpace(mesh, "CG", 1)
 
     bc = DirichletBC(V, Constant(1), "on_boundary")
-    # g = Function(V)
     x, y = SpatialCoordinate(mesh)
     g = project(sin(x)*sin(y), V, annotate=False)
     expr = sin(g*x)
@@ -289,6 +273,4 @@ def test_projection_function():
     h = Function(V)
     h.vector()[:] = rand(h.dof_dset.size)
     m = g.copy(deepcopy=True)
-    g.block_variable.tlm_value = h
-    tape.evaluate_tlm()
-    assert (taylor_test(Jhat, m, h, dJdm=J.block_variable.tlm_value) > 1.9)
+    assert (taylor_test(Jhat, m, h, dJdm=Jhat.tlm(h)) > 1.9)
diff --git a/tests/firedrake/output/test_adjoint_disk_checkpointing.py b/tests/firedrake/output/test_adjoint_disk_checkpointing.py
new file mode 100644
index 0000000000..37747e36e1
--- /dev/null
+++ b/tests/firedrake/output/test_adjoint_disk_checkpointing.py
@@ -0,0 +1,105 @@
+from firedrake import *
+from firedrake.__future__ import *
+from pyadjoint import (ReducedFunctional, get_working_tape, stop_annotating,
+                       pause_annotation, Control)
+import numpy as np
+import os
+import pytest
+
+
+@pytest.fixture(autouse=True, scope="module")
+def handle_annotation():
+    from firedrake.adjoint import annotate_tape, continue_annotation
+    if not annotate_tape():
+        continue_annotation()
+    yield
+    # Ensure annotations are paused at the end of the module.
+    annotate = annotate_tape()
+    if annotate:
+        pause_annotation()
+
+
+def adjoint_example(mesh):
+    # This example is designed to exercise all the block types for which
+    # the disk checkpointer does something.
+    dg_space = FunctionSpace(mesh, "DG", 1)
+    cg_space = FunctionSpace(mesh, "CG", 2)
+    W = dg_space * cg_space
+
+    w = Function(W)
+
+    x, y = SpatialCoordinate(mesh)
+    # AssembleBlock
+    m = assemble(interpolate(sin(4*pi*x)*cos(4*pi*y), cg_space))
+
+    u, v = w.subfunctions
+    # FunctionAssignBlock, FunctionMergeBlock
+    v.assign(m)
+    # SubfunctionBlock, GenericSolveBlock
+    u.project(v)
+
+    # AssembleBlock
+    J = assemble((u - v)**2 * dx)
+    Jhat = ReducedFunctional(J, Control(m))
+
+    with stop_annotating():
+        m_new = assemble(interpolate(sin(4*pi*x)*cos(4*pi*y), cg_space))
+    checkpointer = get_working_tape()._package_data
+    init_file_timestamp = os.stat(checkpointer["firedrake"].init_checkpoint_file.name).st_mtime
+    current_file_timestamp = os.stat(checkpointer["firedrake"].current_checkpoint_file.name).st_mtime
+    Jnew = Jhat(m_new)
+    # Check that any new disk checkpoints are written to the correct file.
+    assert init_file_timestamp == os.stat(checkpointer["firedrake"].init_checkpoint_file.name).st_mtime
+    assert current_file_timestamp < os.stat(checkpointer["firedrake"].current_checkpoint_file.name).st_mtime
+
+    assert np.allclose(J, Jnew)
+
+    grad_Jnew = Jhat.derivative(apply_riesz=True)
+
+    return Jnew, grad_Jnew
+
+
+@pytest.mark.skipcomplex
+@pytest.mark.parallel(nprocs=3)
+def test_disk_checkpointing():
+    from firedrake.adjoint import enable_disk_checkpointing, \
+        checkpointable_mesh, pause_disk_checkpointing, continue_annotation
+    tape = get_working_tape()
+    tape.clear_tape()
+    enable_disk_checkpointing()
+    continue_annotation()
+    mesh = checkpointable_mesh(UnitSquareMesh(10, 10, name="mesh"))
+    J_disk, grad_J_disk = adjoint_example(mesh)
+    tape.clear_tape()
+    pause_disk_checkpointing()
+
+    J_mem, grad_J_mem = adjoint_example(mesh)
+
+    assert np.allclose(J_disk, J_mem)
+    assert np.allclose(assemble((grad_J_disk - grad_J_mem)**2*dx), 0.0)
+    tape.clear_tape()
+
+
+@pytest.mark.skipcomplex
+def test_disk_checkpointing_successive_writes():
+    from firedrake.adjoint import enable_disk_checkpointing, \
+        checkpointable_mesh, pause_disk_checkpointing
+    tape = get_working_tape()
+    tape.clear_tape()
+    enable_disk_checkpointing()
+
+    mesh = checkpointable_mesh(UnitSquareMesh(1, 1))
+
+    cg_space = FunctionSpace(mesh, "CG", 1)
+    u = Function(cg_space, name='u')
+    v = Function(cg_space, name='v')
+
+    u.assign(1.)
+    v.assign(v + 2.*u)
+    v.assign(v + 3.*u)
+
+    J = assemble(v*dx)
+    Jhat = ReducedFunctional(J, Control(u))
+    assert np.allclose(J, Jhat(Function(cg_space).interpolate(1.)))
+    pause_disk_checkpointing()
+    tape.clear_tape()
diff --git a/tests/firedrake/regression/test_adjoint_operators.py b/tests/firedrake/regression/test_adjoint_operators.py
index 6410d74957..f995295674 100644
--- a/tests/firedrake/regression/test_adjoint_operators.py
+++ b/tests/firedrake/regression/test_adjoint_operators.py
@@ -863,10 +863,10 @@ def test_assign_zero_cofunction():
     J = assemble(((sol + Constant(1.0)) ** 2) * dx)
     # The zero assignment should break the tape and hence cause a zero
     # gradient.
-    grad_l2 = compute_gradient(J, Control(k), options={"riesz_representation": "l2"})
-    grad_none = compute_gradient(J, Control(k), options={"riesz_representation": None})
-    grad_h1 = compute_gradient(J, Control(k), options={"riesz_representation": "H1"})
-    grad_L2 = compute_gradient(J, Control(k), options={"riesz_representation": "L2"})
+    grad_l2 = compute_derivative(J, Control(k, riesz_map="l2"), apply_riesz=True)
+    grad_none = compute_derivative(J, Control(k), apply_riesz=False)
+    grad_h1 = compute_derivative(J, Control(k, riesz_map="H1"), apply_riesz=True)
+    grad_L2 = compute_derivative(J, Control(k, riesz_map="L2"), apply_riesz=True)
     assert isinstance(grad_l2, Function) and isinstance(grad_L2, Function) \
         and isinstance(grad_h1, Function)
     assert isinstance(grad_none, Cofunction)
@@ -912,7 +912,6 @@ def test_riesz_representation_for_adjoints():
     space = FunctionSpace(mesh, "Lagrange", 1)
     f = Function(space).interpolate(SpatialCoordinate(mesh)[0])
     J = assemble((f ** 2) * dx)
-    rf = ReducedFunctional(J, Control(f))
     with stop_annotating():
         v = TestFunction(space)
         u = TrialFunction(space)
@@ -931,21 +930,27 @@ def test_riesz_representation_for_adjoints():
         dJdu_function_L2 = Function(space)
         solve(a == dJdu_cofunction, dJdu_function_L2)
 
-    dJdu_none = rf.derivative(options={"riesz_representation": None})
-    dJdu_l2 = rf.derivative(options={"riesz_representation": "l2"})
-    dJdu_H1 = rf.derivative(options={"riesz_representation": "H1"})
-    dJdu_L2 = rf.derivative(options={"riesz_representation": "L2"})
-    dJdu_default_L2 = rf.derivative()
-    assert (
-        isinstance(dJdu_none, Cofunction) and isinstance(dJdu_function_l2, Function)
-        and isinstance(dJdu_H1, Function) and isinstance(dJdu_default_L2, Function)
-        and isinstance(dJdu_L2, Function)
-        and np.allclose(dJdu_none.dat.data, dJdu_cofunction.dat.data)
-        and np.allclose(dJdu_l2.dat.data, dJdu_function_l2.dat.data)
-        and np.allclose(dJdu_H1.dat.data, dJdu_function_H1.dat.data)
-        and np.allclose(dJdu_default_L2.dat.data, dJdu_function_L2.dat.data)
-        and np.allclose(dJdu_L2.dat.data, dJdu_function_L2.dat.data)
+    dJdu_none = ReducedFunctional(J, Control(f)).derivative()
+    dJdu_l2 = ReducedFunctional(J, Control(f, riesz_map="l2")).derivative(
+        apply_riesz=True
+    )
+    dJdu_H1 = ReducedFunctional(J, Control(f, riesz_map="H1")).derivative(
+        apply_riesz=True
+    )
+    dJdu_L2 = ReducedFunctional(J, Control(f, riesz_map="L2")).derivative(
+        apply_riesz=True
+    )
+    dJdu_default_L2 = ReducedFunctional(J, Control(f)).derivative(
+        apply_riesz=True
     )
+    assert isinstance(dJdu_none, Cofunction) and isinstance(dJdu_function_l2, Function)
+    assert isinstance(dJdu_H1, Function) and isinstance(dJdu_default_L2, Function)
+    assert isinstance(dJdu_L2, Function)
+    assert np.allclose(dJdu_none.dat.data, dJdu_cofunction.dat.data)
+    assert np.allclose(dJdu_l2.dat.data, dJdu_function_l2.dat.data)
+    assert np.allclose(dJdu_H1.dat.data, dJdu_function_H1.dat.data)
+    assert np.allclose(dJdu_default_L2.dat.data, dJdu_function_L2.dat.data)
+    assert np.allclose(dJdu_L2.dat.data, dJdu_function_L2.dat.data)
 
 
 @pytest.mark.skipcomplex
@@ -970,13 +975,13 @@ def test_lvs_constant_jacobian(constant_jacobian):
 
     J_hat = ReducedFunctional(J, Control(u))
 
-    dJ = J_hat.derivative(options={"riesz_representation": None})
+    dJ = J_hat.derivative()
     assert np.allclose(dJ.dat.data_ro, 2 * assemble(inner(u_ref, test) * dx).dat.data_ro)
 
     u_ref = Function(space, name="u").interpolate(X[0] - 0.1)
     J_hat(u_ref)
 
-    dJ = J_hat.derivative(options={"riesz_representation": None})
+    dJ = J_hat.derivative()
     assert np.allclose(dJ.dat.data_ro, 2 * assemble(inner(u_ref, test) * dx).dat.data_ro)
 
 
@@ -997,7 +1002,7 @@ def test_cofunction_assign_functional():
     cof2.assign(cof)  # Test is checking that this is taped.
     J = assemble(action(cof2, f2))
     Jhat = ReducedFunctional(J, Control(f))
-    assert np.allclose(float(Jhat.derivative()), 1.0)
+    assert np.allclose(float(Jhat.derivative(apply_riesz=True)), 1.0)
     f.assign(2.0)
     assert np.allclose(Jhat(f), 2.0)
 
@@ -1035,7 +1040,7 @@ def u_analytical(x, a, b):
         (u_analytical(X[0], a, b)**2) * dx, b))
     J = assemble(sol * sol * dx)
     J_hat = ReducedFunctional(J, [Control(a), Control(b)])
-    adj_derivatives = J_hat.derivative(options={"riesz_representation": "l2"})
+    adj_derivatives = J_hat.derivative()
     assert np.allclose(adj_derivatives[0].dat.data_ro, der_analytical0.dat.data_ro)
     assert np.allclose(adj_derivatives[1].dat.data_ro, der_analytical1.dat.data_ro)
     a = Function(R, val=1.5)
diff --git a/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py b/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
index 9975aa6dd1..5ed29951e6 100644
--- a/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
+++ b/tests/firedrake/vertexonly/test_poisson_inverse_conductivity.py
@@ -59,8 +59,8 @@ def test_poisson_inverse_conductivity(num_points):
     # Compute the true solution of the PDE.
     u_true = Function(V)
     v = TestFunction(V)
-    f = Constant(1.0, domain=m)
-    k0 = Constant(0.5, domain=m)
+    f = Constant(1.0)
+    k0 = Constant(0.5)
     bc = DirichletBC(V, 0, 'on_boundary')
     F = (k0 * exp(q_true) * inner(grad(u_true), grad(v)) - f * v) * dx
     solve(F == 0, u_true, bc)
@@ -79,7 +79,7 @@ def test_poisson_inverse_conductivity(num_points):
     signal_to_noise = 20
     U = u_true.dat.data_ro[:]
     u_range = U.max() - U.min()
-    σ = Constant(u_range / signal_to_noise, domain=point_cloud)
+    σ = Constant(u_range / signal_to_noise)
     ζ = generator.standard_normal(len(xs))
     u_obs_vals = np.array(u_true.at(xs)) + float(σ) * ζ
 
@@ -102,7 +102,7 @@ def test_poisson_inverse_conductivity(num_points):
 
     # Two terms in the functional
     misfit_expr = 0.5 * ((u_obs - assemble(interpolate(u, P0DG))) / σ)**2
-    α = Constant(0.5, domain=m)
+    α = Constant(0.5)
     regularisation_expr = 0.5 * α**2 * inner(grad(q), grad(q))
 
     # Form functional and reduced functional
@@ -111,7 +111,7 @@ def test_poisson_inverse_conductivity(num_points):
     Ĵ = ReducedFunctional(J, q̂)
 
     # Estimate q using Newton-CG which evaluates the hessian action
-    minimize(Ĵ, method='Newton-CG', options={'maxiter': 3, 'disp': True})
+    minimize(Ĵ, method='Newton-CG', options={'maxiter': 10, 'xtol': 1e-1, 'disp': True})
 
 
 @pytest.mark.skipcomplex  # Taping for complex-valued 0-forms not yet done