How do we deal with bounds constraints with the adjoint?

[JHopeCollins]

At the moment we just raise an error:

firedrake/firedrake/adjoint_utils/variational_solver.py

Lines 68 to 71 in 1faf6e6

	bounds = kwargs.pop("bounds", None)
	if bounds is not None:
	raise ValueError(
	"MissingMathsError: we do not know how to differentiate through a variational inequality")

It would be nice if we could deal with bounds, but someone needs to sit down and work out the maths so it DTRT.