[firedrake] Pass approximate Jacobian to derivative() call

[Buesing, Henrik]

> > Is there a multi-physics example where the user supplies a Jacobian to
> > solve()? I think if I knew how this Jacobian should look like in my
> > 2x2 block case, I would manage with the rest.
> 
> firedrake/tests/regression/test_vector_laplace_on_quadrilaterals.py
> does something similar.  It solves with one specified Jacobian, and uses
> a different one as the preconditioner.
> 
> The idea is that the blocks are "implicit" in the fact that your test
> and trial functions are from a mixed space.
> 
> For example:
> 
> W = DG1*DG2
> 
> u, p = TrialFunctions(W)
> v, q = TestFunctions(W)
> 
> J = u*v*dx + p*q*dx
> 
> Builds a jacobian that is, block-wise:
> 
> J = [DG1-mass-matrix, 0;
>      0, DG2-mass-matrix]
> 
> So you just need to write down the symbolic form in terms of things that
> have come from the mixed space.

[Buesing, Henrik] So in my nonlinear case I would use

u=Function(W)
pw, enth=split(u)

and write down my J in terms of pw,enth. Nevertheless I will have four contributions (for the four blocks) in my Jacobian, which I just add together as in your example. 

Every contribution has a (pw + eps) and (pw - eps), or with enth. Can I write down my form as a function dependent on the primary variables. Then I do not need to define all these different +/- eps...

Henrik