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.