Matrix-free solution using PetscNonlinearSolver

[brightzhang91]

Hi :)
I'm building a new NonlinearFormintegrator for nonlinear material modeling. Due to the difficulty in producing the analytical Jacobian, I want to try the PetscNonlinearSolver to perform the matrix-free solution. I had a look at the petsc/ex10p, which can be activated to be matrix-free with -snes_mf or -snes_mf_operator. However, ex10p already has the analytical Jacobian built in the AssembleElementGrad().
By using the matrix-free method, do I still need to explicitly build the GetGradient() of the operator ? How should I construct the AssembleElementGrad() of the new NonlinearFormintegrator, as I try to use the matrix-free Newton-Krylov to avoid deriving the analytical Jacobian ?

Thanks !
Bright

[jandrej]

By using the matrix-free method, do I still need to explicitly build the GetGradient() of the operator ?

Yes. If you need the matrix for preconditioning, you need to assemble it.

You can setup petsc to compute jacobian-vector products on the fly using finite differences and coloring, but this option is not available in combination with mfem right now.

[brightzhang91]

By using the matrix-free method, do I still need to explicitly build the GetGradient() of the operator ?

Yes. If you need the matrix for preconditioning, you need to assemble it.

You can setup petsc to compute jacobian-vector products on the fly using finite differences and coloring, but this option is not available in combination with mfem right now.

Many thanks for your quick reply !
By using -snes_mf_operator, PETSc can evaluate the correct the Jacobian even with incorrect or incomplete analytic Jacobian (https://www.mcs.anl.gov/petsc/petsc-current/docs/manual.pdf 5.12). Therefore I wonder if I can perform matrix-free solution to my nonlinear modeling in MFEM without knowing the exact analytic Jacobian ?
I heard the finite difference Jacobian -snes_fd is very slow, so the matrix-free method is closer to what I want.

Thanks
Bright

[tangqi]

I think snes_mf will use JFNK and it will not use GetGradient() as preconditioner. So for simple problems which has weak nonlinearity, it may be enough. I believe in example 10, the case of rc_ex10p_mf should demonstrate that.

In general, you need to provide the preconditioner for JFNK through GetGradient(). You can start with providing a simple approximation to the Jacobian matrix through GetGradient(), like a block diagonal matrix that only keeps the diagonal small matrices. The case of rc_ex10p_mf demonstrates how to use it.

Yes, snes_fd can be very slow. The normal way to get it working is through snes_fd_color. But it looks like this option is not available in mfem yet.

[jandrej]

You get the matrix-vector-product of the Jacobian with -snes_mf_operator, but this will not allow the preconditioning with matrix based solvers. Jacobi, Krylov methods and the like will work fine regardless.

This works out of the box with MFEM by just setting the option on the command line (e.g. ex10). The gradient that is implemented only uses (if desired) the assembled gradient for the PC.

[jandrej]

Yes, snes_fd can be very slow. The normal way to get it working is through snes_fd_color. But it looks like this option is not available in mfem yet.

The limiting factor in MFEM is that we don't provide a DM object for PETSc from which coloring of matrix could be determined.

[brightzhang91]

Thanks for the comments.
It seems even for the matrix-free method, I'll need to produce an approximate analytic Jacobian matrix in the GetGradient(), with or without the preconditioner. Am I understanding it correctly ?
As for the approximated analytic Jacobian, how does the degree of approximation matter ? such as some terms missing in the Jacobian equations.

thanks

[jandrej]

We can't answer the the question without knowing the details of the equations you are solving. In coupled equations it might be OK to just provide diagonal blocks for the approximation, although it might not be enough in strongly nonlinear cases.

[stefanozampini]

@brightzhang91 If you use -snes_mf, then you do not need to provide a Jacobian via the GetGradient method.
In this case, the Krylov method will be without preconditioner. Instead, if you use -snes_mf_operator you need to provide a Jacobian matrix to be used to setup the preconditioner. The mat-vec within the Krylov method will be matrix free, but you will be able to construct a preconditioner using whatever operator you provide via GetGradient

@jandrej Regarding the coloring: well, the interface from MFEM to PETSc is completely algebraic and specifying mesh information is somewhat out-of-synch with the MFEM approach, since GetGradient is already designed to give access back to the user to provide the Jacobian. On the other hand, users can typecast the PetscNonlinearSolver to SNES, get the DM out of it, and, using private petsc headers, set the dm->ops->getcoloring callback. Cumbersome, I know, but PETSc is designed to be extensible. I can add a specific interface in DMSHELL (like DMShellSetColoring) if you need it

[tzanio]

This issue has not been active for more than two weeks, so I will close it now. Feel free to reopen if you still have questions.