I encourage you to take a look at src/snes/examples/tutorials/ex28.c.  It solves the following coupled system [*]<div><div><br></div><div><div> PDE (U):</div><div>     -(k u_x)_x = 1 on (0,1), subject to u(0) = 0, u(1) = 1</div>
<div> Algebraic (K):</div><div>     exp(k) + k = u + 1/(1 + u_x^2)</div></div><div><br></div><div>and each part separately, with identical assembly code.  An analytic Jacobian is available in all cases.  Change between solving individual physics and coupled with -problem_type 0,1,2.  The assembly interface is particularly notable because the same code can also assemble into a MatNest, which is based, on Dave May's "Block" matrix from PETScExt.  The primary advantage of MatNest is that getting the submatrices used in a fieldsplit preconditioner is trivially cheap.  Currently DMGetMatrix_Composite_Nest does not allocate off-diagonal blocks so we are limited to additive fieldsplits, but when we sort out the preallocation API, this code will work with the other fieldsplit variants.</div>
<div><br></div><div>  cd src/snes/examples/tutorials && make ex28</div><div>  mpiexec -n 2 ./ex28 -da_grid_x 100 -{snes,ksp}_converged_reason -snes_monitor_short -problem_type 2 -pack_dm_mat_type nest -snes_mf_operator -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_u_pc_type ml -fieldsplit_k_pc_type jacobi</div>
<div><br></div><div>This assembles directly into the "best" structure for FieldSplit to work with, there are no intermediate copies.  Note that if the submatrices have extra structure, like (S)BAIJ, you can use them inside the blocks and use MatSetValuesBlockedLocal for assembly.  This blocked assembly will also work when you want to assemble the whole thing and use a direct solver, as in</div>
</div><div><br></div><div>  mpiexec -n 2 ./ex28 -da_grid_x 100 -{snes,ksp}_converged_reason -snes_monitor_short -problem_type 2 -pc_type lu -pc_factor_mat_solver_package mumps</div><div><br></div><div>Jed</div><div><div><br>
</div><div>[*] I have no idea if this "physics" is stable.  I think it is not based on floating-point errors at very high resolution, but if you keep problem size under 1000 or so, the solution does not blow up.  Suggestions for something of similar form, with interesting nonlinearities and a not completely boring solution profile, are certainly welcome.</div>
</div>