Consider a problem where space is filled with a liquid and a solid phase, with a large, complicated geometry. On each of the phases, there's an electrical potential/poisson equation. The equations are coupled by a nonlinear surface current through the phase boundary.
Currently, an FEM is constructed for both phases and assembled into one large Jacobian which is then solved as a black-box in a Newton method. This leads to bad matrix condition and bad performance.
To improve performance, I plan to use PETSc. For multi-physics problems, they suggest the usage of nested matrices and then to use a separate preconditioner for each of the equations. BUT, for all the examples I found so far, there's one big difference: The examples always consider different physical equations on one phase which are coupled in all points in space. In my case, I have just one equation per phase, but a coupling only on the boundaries of the phase, and nonlinear.
How would you deal with that kind of problem in PETSc? (Emphasis on the combination of matrices, solvers and preconditioners.)