Add optional tensor arg to project

If supplied, write result to tensor rather than building a new
function.
Fix Boussinesq demo so it runs

To be able to take the derivative of F wrt out to compute the jacobian,
we need to use UFL's split function (rather than out.split) to build b,
u and p.  Additionally: project the buoyancy into P1xP1 for output;
crank the size of forcing term up by 10^6 (0.5**18 is not very big); and
fiddle a little with solver options.