Preconditioning Overview
 Krylov methods need preconditioning to be efficient (or even effective!).
 Reduces the total number of linear iterations
 Each linear iteration in MOOSE includes a residual evaluation (
computeQpResidual
)
 Krylov methods, in theory, converge in the number of linear iterations equal to the number of unknowns in the system
 Even though the Jacobian is never formed, JFNK methods still require preconditioning.
 MOOSE's automatic (without user intervention) preconditioning is fairly minimal.
 Many options exist for implementing improved preconditioning in MOOSE.
Preconditioned JFNK
 Using right preconditioning, solve
$$$$\mathbf{J} (\vec{u}_i) \mathbf{M}^{1} (\mathbf{M}\delta \vec{u}_{i+1}) = \vec{R}(\vec{u}_i)$$$$
 $$$\mathbf{M}$$$ symbolically represents the preconditioning matrix or process
 Inside GMRES, we only apply the action of $$$\mathbf{M}^{1}$$$ on a vector
 Recall the unpreconditioned JFNK approximation ($$$\mathbf{M}^{1} = \mathbf{I}$$$):
$$$$\mathbf{J} (\vec{u}_i) \vec{v} \approx \frac{\vec{R}(\vec{u}_i + \epsilon \vec{v})  \vec{R}(\vec{u}_i)}{\epsilon}$$$$
 Compare to the rightpreconditioned, matrixfree version:
$$$$\mathbf{J} (\vec{u}_i) \mathbf{M}^{1}\vec{v} \approx \frac{\vec{R}(\vec{u}_i + \epsilon \mathbf{M}^{1}\vec{v})  \vec{R}(\vec{u}_i)}{\epsilon}$$$$
Preconditioning Matrix vs Process
 On the previous slide $$$\mathbf{M}$$$ represented the "Preconditioning Matrix".
 The action of $$$\mathbf{M}^{1}$$$ on a vector represents the "Preconditioner" or "Preconditioning Process".
 In MOOSE the "matrix to build" and the "process to apply" with that matrix are separated.
 There are several different ways to build preconditioning matrices:
 Default: Block Diagonal Preconditioning
 Single Matrix Preconditioner (SMP)
 Finite Difference Preconditioner (FDP)
 Physics Based Preconditioner (PBP)
 Field Split Preconditioner
 After selecting how to build a preconditioning matrix you can then use solver options to select how to apply the Preconditioner.
Solve Type
 The default
solve_type
for MOOSE is "Preconditioned JFNK".
 An alternative
solve_type
can be set through either the [Executioner]
or [Preconditioner/*]
block.
 Valid options include:
PJFNK
(default)
JFNK
NEWTON
FD
(Finite Difference)
PETSc Preconditioning Options
 For specifying the preconditioning process we use solver options directly (i.e. PETSc options).
 Currently the options for preconditioning with PETSc are exposed to the applications.
 This will change in the future... there will be more generic ways of specifying preconditioning parameters.
 The best place to learn about all of the preconditioning options with PETSc is the user manual.
 We use the commandline syntax, but provide places to enter it into the input file.
http://www.mcs.anl.gov/petsc/petsccurrent/docs/manual.pdf
PETSc Specific Executioner Options
petsc_options 
Description 
snes_ksp_ew 
Variable linear solve tolerance, useful for transient solves 
help 
Show PETSc options during the solve 
petsc_options_iname 
petsc_options_value 
Description 
pctype 
ilu 
Default for serial 

bjacobi 
Default for parallel with sub_pc_type ilu 

asm 
Additive Schwartz with sub_pc_type ilu 

lu 
Full LU, serial only 

gamg 
Generalized GeometricAlgebric MultiGrid 

hypre 
Hypre, usually used with boomeramg 
sub_pc_type 
ilu, lu, hypre 
Can be used with bjacobi or asm 
pc_hypre_type 
boomeramg 
Algebraic Multigrid 
ksp_gmres_restart 
# (default = 30) 
Number of Krylov vectors to store 



Default Preconditioning Matrix

Consider the fully coupled system of equations:
$$$$\begin{align*} \nabla \cdot k(T,s) \nabla T &= 0 \\  \nabla \cdot D(T,s) \nabla s &= 0 \end{align*}$$$$

Fully coupled Jacobian approximation
$$$$\begin{equation*}
\mathbf{J}(T,s) =
\begin{bmatrix}
\frac{\partial (R_T)_i}{\partial T_j} & \frac{\partial (R_T)_i}{\partial s_j} \\
\frac{\partial (R_s)_i}{\partial T_j} & \frac{\partial (R_s)_i}{\partial s_j}
\end{bmatrix}
\approx
\begin{bmatrix}
\frac{\partial (R_T)_i}{\partial T_j} & \mathbf{0} \\
\mathbf{0} & \frac{\partial (R_s)_i}{\partial s_j} \\
\end{bmatrix}
\end{equation*}$$$$

For our example:
$$$$\begin{equation*}
\mathbf{M} \equiv
\begin{bmatrix}
(k(T,s) \nabla \phi_j, \nabla \psi_i) & \mathbf{0} \\
\mathbf{0} & (D(T,s) \nabla \phi_j, \nabla\psi_i)
\end{bmatrix} \approx \mathbf{J}(T,s)
\end{equation*}$$$$

This simple style of throwing away the offdiagonal blocks is the way MOOSE will precondition when using the default solve_type
.
The Preconditioning Block
[Preconditioning]
active = 'my_prec'
[./my_prec]
type = SMP
#SMP Options Go Here!
#Override PETSc Options Here!
[../]
[./other_prec]
type = PBP
#PBP Options Go Here!
#Override PETSc Options Here!
[../]
[]
 The
Preconditioning
block allows you to define which type of preconditioning matrix to build and what process to apply.
 You can define multiple blocks with different names, allowing you to quickly switch out preconditioning options.
 Each subblock takes a
type
parameter to specify the type of preconditioning matrix.
 Within the subblocks you can also provide other options specific to that type of preconditioning matrix.
 You can also override PETSc options here.
 Only one block can be active at a time.
Single Matrix Preconditioning (SMP)
 Single Matrix Preconditioner (SMP) builds one preconditioning matrix.
 You enable SMP with:
type = SMP
 You specify which blocks of the matrix to use with:
off_diag_row = 's'
off_diag_column = 'T'
 Which would produce an $$$\mathbf{M}$$$ like this:
$$$$\mathbf{M} \equiv
\begin{bmatrix}
\left(k(T,s) \nabla \phi_j, \nabla \psi_i\right) & \mathbf{0}
\\[3pt]
\left(\frac{\partial D(T,s)}{\partial T_j} \nabla s, \nabla \psi_i\right) & \left(D(T,s) \nabla \phi_j, \nabla\psi_i\right)
\end{bmatrix} \approx \mathbf{J}$$$$
 In order for this to work you must provide a
computeQpOffDiagJacobian()
function in your Kernels that computes the required partial derivatives.
 To use all off diagonal blocks set:
full = true
.
Finite Difference Preconditioning (FDP)
 The Finite Difference Preconditioner (FDP) allows you to form a "Numerical Jacobian" by doing direct finite differences of your residual statements.
 This is extremely slow and inefficient, but is a great debugging tool because it allows you to form a nearly perfect preconditioner.
 You specify it by using:
type = FDP
 You can use the same options for specifying offdiagonal blocks as SMP.
 Since FDP allows you to build the perfect approximate Jacobian it can be useful to use it directly to solve instead of using JFNK.
 The finite differencing is sensitive to the differencing parameter which can be specified using:
petsc_options_iname = 'mat_fd_coloring_err mat_fd_type'
petsc_options_value = '1e6 ds'
 NOTE: FDP currently works in serial only! This might change in the future, but FDP will always be meant for debugging purposes!
Examples
 Default Preconditioning Matrix, Preconditioned JFNK, monitor linear solver, variable linear solver tolerance.
 Use Hypre with algebraic multigrid and store 101 Krylov vectors.
[Executioner]
...
petsc_options_iname = 'pc_type pc_hypre_type ksp_gmres_restart'
petsc_options_value = 'hypre boomeramg 101'
...
[]
 Single Matrix Preconditioner, Fill in the (forced, diffused) block, Preconditioned JFNK, Full inverse with LU
[./SMP_jfnk]
type = SMP
off_diag_row = 'forced'
off_diag_column = 'diffused'
petsc_options_iname = 'pc_type'
petsc_options_value = 'lu'
[../]
Example 11
Physics Based Preconditioning
 Physics based preconditioning is an advanced concept used to more efficiently solve using JFNK.
 The idea is to create a preconditioning process that targets each physics individually.
 In this way you can create a more effective preconditioner...while also maintaining efficiency.
 In MOOSE there is a
PhysicsBasedPreconditioner
object.
 This object allows you to dial up a preconditioning matrix and the operations to be done on the different blocks of that matrix on the fly from the input file.
What the PBP Does
 The PBP works by partially inverting a preconditioning matrix (usually an approximation of the true Jacobian) by partially inverting each block row in a BlockGaussSeidel way.
$$$$\vec{R}(u,v) =
\begin{bmatrix}
\vec{r}_u
\\
\vec{r}_v
\end{bmatrix}$$$$
$$$$\mathbf{M} \equiv
\begin{bmatrix}
\frac{\partial (R_u)_i}{\partial u_j} & \mathbf{0} \\
\frac{\partial (R_v)_i}{\partial u_j} & \frac{\partial (R_v)_i}{\partial v_j}
\end{bmatrix} \approx \mathbf{J}$$$$
$$$$\mathbf{M} \vec{q} = \vec{p} \quad \Rightarrow \quad
\begin{cases}
\frac{\partial (R_u)_i}{\partial u_j} \vec{q}_u = \vec{p}_u \\[6pt]
\frac{\partial (R_v)_i}{\partial v_j} \vec{q}_v = \vec{p}_v  \frac{\partial (R_v)_i}{\partial u_j} \vec{q}_u
\end{cases}$$$$
Using the PBP
[Variables]
...
[]
[Preconditioning]
active = 'myPBP'
[./myPBP]
type = PBP
solve_order = 'u v'
preconditioner = 'ILU AMG'
off_diag_row = 'v'
off_diag_column = 'u'
[../]
[]
 Set up a PBP object for a two variable system (consisting of variables "u" and "v").
 Use ILU for the "u" block and AMG for "v".
 Use the lower diagonal (v, u) block.
 When using
type = PBP
, will set solve_type = JFNK
automatically.
Applying PBP
 Applying these ideas to a coupled thermomechanics problem
Example 12