[petsc-users] Preconditioning in matrix free SNES

[Travis C. Fisher]

I am trying to precondition a matrix free SNES solution. The application 
is a high order compressible Navier Stokes solver. My implementation is:

I use the -snes_mf_operator option and create a matrix free matrix via 
MatCreateSNESMF. The function SNESSetJacobian points to (FormJacobian) 
uses MatMFFDComputeJacobian to calculate the matrix free jacobian. I 
then use SNESSetFromOptions and extract the linear solver and 
preconditioner. The linear solver is GMRES. The matrix free jacobian 
solve with no preconditioner works pretty well for the smooth test 
problem I am currently working on, but in general I don't expect it to 
perform that well without preconditioning. I have attempted to perform 
the preconditioning in two ways:

1) Simply setting the preconditioner extracted from the SNES to PCLU. 
The preconditioner matrix I calculate in my FormJacobian routine is a 
frozen Jacobian matrix, so for the first linear solve, the 
preconditioned operator is essentially the identity operator. The code 
performs an LU decomposition, but the effect of the preconditioner is a 
larger linear residual, so clearly something is wrong.

2) I set the preconditioner from the SNES context to PCSHELL. I then 
make calls to PCShellSetSetup and PCShellSetApply to assign the routines 
for setting up and applying the PC. The routine for PCSetUp creates a 
new LU preconditioner and sets the operator the frozen jacobian matrix 
as described above via PCSetOperators. The PCApply just calls PCApply. 
Again the code performs the LU decomposition, but I get the same result 
as above.

I realize this may point to my preconditioner matrix being completely 
wrong, but it "looks" right when I check values. It would take a lot of 
effort for me to set up the coloring to have petsc calculate the 
jacobian via finite differences since I am using a high order stencil 
with full boundary closures.

My question is am I obviously doing something incorrectly? Am I somehow 
failing to direct the SNES context to apply the LU decomposition and not 
assume that I have given it a preconditioner matrix to simply perform a 
matrix multiply? I appreciate any direction you may be able to give me.

Thanks,

Travis Fisher