Multigrid preconditionner for GMRES in a Jacobian-Free Newton-Krylov method

[bouloumag]

I am developping a Jacobian-Free Newton-Krylov code to solve the
Navier Stokes equations. In this code, the product of the Jacobian
matrix with a given vector is represented by the matvec method of a
"LinearOperator" object.

The pyamg.krylov.gmres method can takes an optional preconditioner as
parameter which is either an object "LinearOperator" or a matrix.
However, the pyamg documentation and examples does not explain how one
can form this preconditioner in a matrix-free manner given only the
matvec method of the jacobian. It is not possible for me to store the
whole Jacobian as in the example page, because I am working on a large
grid.

Any suggestions / example code for this would be very appreciated.

Thanks,

Christine

[Luke]

Christine,

You're correct. AMG does require the matrix to be built. The reason
is that coarse grid operators are not explicitly known and need to be
formed through the triple matrix product A_c = R*A*P.

One thing you can do is force the grid hierarchy (so that it's no
longer algebraic) and then define your LinearOperator for each level
of the hierarchy (a matrix-free geometric multigrid-type approach).
PyAMG is not really built for this, so you may get more mileage by
building your own cycling or just ripping the cycling and multilevel
containers out of PyAMG for use in your own specialized matrix-free
code.

hth,
Luke

[bouloumag]

Thank you for your answer. According to some papers, it seems that
gmres could be used to precondition fgmres in a way that is completely
matrix-free. Is it also possible with the fgmres implementation
provided by pyamg ?

Christine

[Luke]

On Jan 10, 10:15 pm, bouloumag <boulou...@gmail.com> wrote:
> Thank you for your answer. According to some papers, it seems that
> gmres could be used to precondition fgmres in a way that is completely
> matrix-free. Is it also possible with the fgmres implementation
> provided by pyamg ?
>

Yep. The fGMRES in PyAMG works just fine as a stand alone
preconditioned Krylov method outside of PyAMG. In fact, many use the
Krylov solvers in PyAMG just for the Krylov solver (and not PyAMG). I
would look at building your own hierarchy with multilevel.py (defining
your own mat-vec) and then using the fGMRES.

Good luck. And keep us posted...

[Jacob Schroder]

Hi Christine,

Luke is right that your best hope for a robust preconditioner is to
manually construct a geometric multigrid type hierarchy where the
matrix is implicitly defined at each level by a LinearOperator.
However, if you simply want to precondition fGMRES with GMRES, then
you'll have to construct a LinearOperator whose mat-vec routine is
defined by the action of GMRES. You'll then pass this new
LinearOperator in as the parameter M to fGMRES. Since it's somewhat
tricky, here's a code snippet showing exactly how to define such an M.

--------------------------------------------------------------------
--------------------------------------------------------------------
##
# Define your own matrix A and right-hand-side b
##

##
# GMRES parameters
zz = numpy.zeros_like(b)
gmres_tol = scipy.finfo(A.dtype).tiny

##
# Define GMRES precondtioner as a Linear Operator
# --> maxiter for GMRES will likely need to be tuned. I use
# maxiter=10 below as only a guess
def precon_matvec(x):
# Basically, GMRES approximates inverse(A) with a polynomial
# P(A) so that x = x0 + P(A)*r, where r is the residual. Thus in the
# call to GMRES, x0 is a zero vector and the right-hand-side is x.
return pyamg.krylov.gmres(A, x, x0=zz, tol=gmres_tol, maxiter=10)[0]
#
M=scipy.sparse.linalg.interface.LinearOperator(A.shape,precon_matvec)

##
# Now call fGMRES with the above M
##
--------------------------------------------------------------------
--------------------------------------------------------------------

Hope this helps,

Jacob

> --
> You received this message because you are subscribed to the Google Groups "pyamg-user" group.
> To post to this group, send email to pyamg...@googlegroups.com.
> To unsubscribe from this group, send email to pyamg-user+...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/pyamg-user?hl=en.
>
>

[bouloumag]

Thank you very much for the code example. I will try this. Really
appreciated.