Preconditioning multiphysics problems is a complicated issue... but I'll try to shed some light. For further reading there is a fairly lengthy explanation of all of this in the MOOSE Workshop Manual which you can snag here: http://mooseframework.org/documentation/
Note that in MOOSE we separate out the ideas of the "Preconditioning Matrix" (which is what we call a
Preconditioner in MOOSE) and the "Preconditioning Process" which is the approximate solve done by something like PETSc to invert the Preconditioning Matrix... for specifying that you utilize PETSc options (look at the Workshop Manual)
The first rule of Preconditioning in MOOSE is to do the simplest thing possible that allows you to solve your problem while you're developing your application. Only after your application is mostly "complete" (meaning you have all of the relevant physics in there and you are basically getting the answers you want) should you start delving into the more complicated preconditioning options.
Because of this we have a hierarchy of Preconditioners spanning from simple to more complicated in this order:
FDP) - only use for debugging
As noted above, the default Preconditioner is a "Block Diagonal Preconditioner". This will basically call your
computeQpJacobian() statements and fill in the diagonal block for each variable in your problem. This basically creates a "decoupled" preconditioner that works well on a variety of problems. You can then use PETSc options (as described in the Workshop Manual) to apply PETSc preconditioning processes (like ILU or HYPRE-boomeramg) to that matrix.
Single Matrix Preconditioner (SMP) will fill the block diagonals by default and then allow you to specify some off-diagonal blocks to also fill in. This is accomplished by specifying
off_diag_column. Let's say you have two variables in your system: u and v. And v shows up strongly in the u equation (like a v^2 is in the u PDE somewhere). It might be important to tell the preconditioner about that relationship so you could do it like so:
[Preconditioning] [./some_smp] type = SMP off_diag_row = 'u' off_diag_column = 'v' [../] 
off_diag_column are matched "up and down" as in the first thing in
off_diag_row goes with the first thing in
off_diag_column... the second entries in both go together... the third entries go together... etc.
This will cause the SMP to call your
computeQpOffDiagJacobian() statements for your "u" Kernels and ask them to compute the derivative of your PDE terms with respect to "v".... filling in the u-v block of your matrix.
You can then use PETSc options to apply a PETSc preconditioning process to that matrix (see the Workshop Manual for more information).
Finite Difference Preconditioner (FDP) takes the same options as SMP but does not call your
computeQpOffDiagJacobian() statements.... instead it does a numerical finite differencing to fill in your preconditioning matrix. This only works in serial and is just there to help you track down problems.
Split are both "advanced" preconditioners. Each of them allows you to dial up a different PETSc preconditioning process to apply to each equation separately and then try to combine them back together. This allows you to, for instance, use ILU to target a convection dominated PDE in your system while utilizing an Algebraic Multigrid (AMG) preconditioner on your more elliptic looking PDEs.
PBP takes the same options as SMP and FDP for specifying off-diagonal blocks but also takes a
solve_order option and a
preconditioner option that are matched together.
solve_order is the order you want to solve each equation while
preconditioner is the type of preconditioning process (like ILU or AMG) that you want to use for that solve. solve is highlighted here because it's really a approximate solve. PBP basically gives you the ability to dial up a "block Gauss-Seidel" solve for your preconditioner.
For examples of all of these you can look in the
examples directory at example 11 and 12. You can also look in
moose/test/tests/preconditioners for additional examples.