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