Allow modifying coefficients of reduced functionals

[APaganini]

Describe the current issue
If a reduced functional depends on a coefficient that (i) belongs to a FunctionSpace and (ii) is not declared as a Control, modifying the value of the coefficient after the reduced functional has been created does not affect the reduced functional, that is, the dependence is lost. For example,

from firedrake import *
from firedrake.adjoint import *
continue_annotation()
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)
P = FunctionSpace(mesh, "DG", 0)
f = Function(P)
R = FunctionSpace(mesh, "Real", 0) # this does not work
c = Function(R)  # coefficient to be modified
F = (u-(f+c))*v*dx
solve(F == 0, u)
J = assemble(u**2*dx)
Jred = ReducedFunctional(J, Control(f))
stop_annotating()
print(Jred(f)) # output is zero 0
c.interpolate(1.)
print(Jred(f)) # output is still 0

Describe the solution you'd like
It would be nice to have an easy way to modify the coefficients. This would be useful in many applications.

Describe alternatives you've considered

1. The example above can be fixed by using c = Constant(0) and c.assign(1), but this does not resolve the problem for coefficients that are not constant.

2. Colin has a workaround based on declaring all coefficients as controls and then carefully wrapping pyadjoint.minimize.

3. One could create a different tape for each value of the coefficient.

Additional info
Daiane commented that enabling this would be unsafe. David probably has strong opinions about this.