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
