[firedrake] Pass approximate Jacobian to derivative() call

[Buesing, Henrik]

> For example here:
> 
> delta = Constant(1e-8)
> tmp = pw
> eps = delta + delta*tmp
> eps = assemble(eps)
> pw = tmp + eps
> val1 = F
> pw = tmp - eps
> val2 = F
> pw = tmp
> tmp00 = Constant(0.5)*(val1-val2)/eps
> 
> it looks like you want val1 to be F evaluated at tmp + eps and val2 to
> be F evaluated a tmp - eps, but in fact, this is just symbolic
> assignment so both val1 and val2 are just different names for F.
> 
> To write your approximate symbolic jacobian, you need to write J in
> terms of functions of the coefficients and test and trial functions
> (just like you did F).  Then it will have the appropriate block
> structure, I think.
> 
[Buesing, Henrik] Ok. Let's assume u=(pw,enth) are my unknowns. F(u)=0 and G(u)=0 are my two pdes formulated as weak forms. Now my Jacobian would have the form

J=[dF(u)/dpw, dF(u)/denth

 dG(u)/dpw, dG(u)/denth]

I guess, I will manage to write down the individual blocks in functional form. But how do I write down the Jacobian in terms of these blocks?

Or maybe I define H = F + G. Then J = [dH(u)/dpw, dH(u)/d_enth] and I get the appropriate blocks automatically. 

Or maybe even J = [dH(u)/du], but here I do not know anymore how to write this down in functional form...

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. 

Thank you!


> 
> I don't quite understand what all the bits are doing, but it looks
> like, for example the dF/dpw block is just zero (since F doesn't
> depend on pw?)

 [Buesing, Henrik] This is just a truncated example. In the end there will be non-zeros in every block.