Doubts about SNESMF and matrix-free preconditioners

[Matthew Knepley]

On Tue, Nov 18, 2008 at 9:59 AM, Rafael Santos Coelho
<rafaelsantoscoelho at gmail.com> wrote:
> Hi, everyone (:
>
> For a considerable amount of time, I've used PETSc a lot in my numerical
> experiments, but then I stopped for a while and got rusty at it. So now, I'm
> trying to get back on track by implementing my own serial matrix-free
> preconditioner for the JFNK (Jacobian-free Newton-Krylov) method. Basicly,
> it's a uniparametric LU-SGS (LU Symmetric Gauss-Seidel) preconditioner as
> bellow:
>
> Let J(x), F(x) and di be the jacobian matrix, the nonlinear function and ith
> diagonal element of J(x), respectively
> Suppose J(x) = L + D + U and consider J(x)s = -F(x), where 's' is the Newton
> correction
> Let M = (D + wL)D^(-1)(D + wU), where w >0
> Now, we have M^(-1)J(x)s = -M^(-1)F(x), knowing that M^(-1)J(x)s = (F(x +
> hM^(-1)s) - F(x)) / h
> Let y = M^(-1)s = (D + wU)^(-1)D(D + wL)^(-1)s and z = (D + wL)^(-1)s <=> (D
> + wL)z = s
> To find z, we follow steps 7 and 8
> x' = x, z1 = s1 / d1,  x'1 = x1 + hz1
> for(i = 2; i <= n; i++) { zi = (si - w(Fi(x') - Fi(x)) / h) / di ; x'i = xi
> + hzi ; }
> Now, we have y = (D + wU)^(-1)Dz
> Let v = Dz and we got y = (D + wU)^(-1)v
> Let u = (D + wU)^(-1)v <=> (D + wU)u = v, and to find u we go through steps
> 12 and 13
> x' = x, un = vn / dn,  x'n = xn + hvn
> for(i = n - 1; i >= 1; i--) { ui = (vi - w(Fi(x') - Fi(x)) / h) / di ; x'i =
> xi + hui ; }
> Finally, we compute M^(-1)J(x)s = (F(x + hy) - F(x)) / h
>
> As for the right-hand side, a similar approach can be done. I only need to
> store diag[J(x)]. Well, I'm not sure how to code the steps above. However,
> looking at "src/snes/examples/tutorials/ex6.c", things got cleared out a
> little bit for me. I think it suffices to type "-snes_mf" on the
> command-line, to write my own "MatrixFreePreconditioner" with all the steps
> above and then
>
> ierr = PCSetType(pc, PCSHELL); CHKERRQ(ierr);
> ierr = PCShellSetApply(pc, MatrixFreePreconditioner); CHKERRQ(ierr);
>
> Is my reasoning correct?

I did not read the PC discussion closely, however, you can do two things:

1) Use -snes_mf_operator described here:

  http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/docs/manualpages/SNES/SNESCreate.html

What this means is that matrix-vector products with the Jacobian are performed
using finite differences with F, and then a preconditioner is computed
using a simpler
matrix you specify using SNESSetJacobian(). This is easy, but not what
you describe above.

2) Use -snes_mf, which again does matvecs with FD, but now no
preconditioner is set by
default. Then you would have to extract the KSP and set your
matrix-free PC there. The
matrix it will receive will be the MFFD that is created by -mat_mf.

  Matt

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener