Doubts about SNESMF and matrix-free preconditioners
Rafael Santos Coelho
rafaelsantoscoelho at gmail.com
Tue Nov 18 09:59:56 CST 2008
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:
1. Let J(x), F(x) and di be the jacobian matrix, the nonlinear function
and ith diagonal element of J(x), respectively
2. Suppose J(x) = L + D + U and consider J(x)s = -F(x), where 's' is the
Newton correction
3. Let M = (D + wL)D^(-1)(D + wU), where w >0
4. 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
5. Let y = M^(-1)s = (D + wU)^(-1)D(D + wL)^(-1)s and z = (D + wL)^(-1)s
<=> (D + wL)z = s
6. To find z, we follow steps 7 and 8
7. x' = x, z1 = s1 / d1, x'1 = x1 + hz1
8. for(i = 2; i <= n; i++) { zi = (si - w(Fi(x') - Fi(x)) / h) / di ; x'i
= xi + hzi ; }
9. Now, we have y = (D + wU)^(-1)Dz
10. Let v = Dz and we got y = (D + wU)^(-1)v
11. Let u = (D + wU)^(-1)v <=> (D + wU)u = v, and to find u we go through
steps 12 and 13
12. x' = x, un = vn / dn, x'n = xn + hvn
13. for(i = n - 1; i >= 1; i--) { ui = (vi - w(Fi(x') - Fi(x)) / h) / di
; x'i = xi + hui ; }
14. 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
1. ierr = PCSetType(pc, PCSHELL); CHKERRQ(ierr);
2. ierr = PCShellSetApply(pc, MatrixFreePreconditioner); CHKERRQ(ierr);
Is my reasoning correct?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20081118/9801bba7/attachment.htm>
More information about the petsc-users
mailing list