Re: [petsc-users] Example of SNES using matrix free jacobian

[Daniel Pino Munoz]

Hi Barry,

Thanks for getting back to me!
I tried replacing KSPSetOperators(ksp, J, J); by 
SNESSetJacobian(snes,J,J, MatMFFDComputeJacobian)

and I get the same result = It works in Debug mode but not in Release. I 
also ran valgrind and it did not catch any memory problem.

Any ideas?

PS : You are right regarding the number of iterations of the non 
preconditioned problem. In the previous version of the code that only 
used a KSP, I already had to set -ksp_gmres_restart 100. But thanks for 
the heads up.

Best,

  Daniel

On 29/10/2024 20:01, Barry Smith wrote:
    Don't call

    KSPSetOperators(ksp, J, J);


     instead call

     SNESSetJacobian(snes,J,J, MatMFFDComputeJacobian)

     but I am not sure that would explain the crash.

     BTW: since you are applying no preconditioner if the matrix is 
ill-conditioned it may take many iterations or not converge. You can try 
something like -ksp_gmres_restart 100 or similar value to try to improve 
convergence (default is 30).

     Barry




On Oct 29, 2024, at 12:37 PM, Daniel Pino Munoz 
<daniel.pino_mu...@mines-paristech.fr> wrote:

Dear all,

I have a linear problem that I am currently solving with a KSP matrix-free.

I would like to move on to a non linear problem, so figure I could start by 
solving the same linear problem using SNES. So I am setting the problem as 
follows:

SNESCreate(PETSC_COMM_WORLD, &snes);
MatCreateShell(PETSC_COMM_WORLD, n_dofs, n_dofs, PETSC_DETERMINE, PETSC_DETERMINE, 
&ctx, &J);
MatCreateShell(PETSC_COMM_WORLD, n_dofs, n_dofs, PETSC_DETERMINE, PETSC_DETERMINE, 
&ctx, &B);
MatShellSetOperation(J, MATOP_MULT, (void (*)(void))(Multiplication));
MatCreateVecs(J, &x_sol, &b);
VecDuplicate(x_sol, &r);
SNESSetFromOptions(snes);
SNESSetFunction(snes, r, &(computeResidual), &ctx);
SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE);
SNESGetLineSearch(snes, &linesearch);
SNESGetKSP(snes, &ksp);
KSPSetOperators(ksp, J, J);
KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);

I tested it with a small problem (compiled in debug) and it works.

When I compiled it in Release, it crashes with a segfault. I tried running the 
Debug version through valgrind, but even for a small problem, it is too slow. 
So I was wondering if you guys could see any rocky mistake on the lines I used 
above?

Otherwise, is there any example that uses a SNES combined with a matrix free 
KSP operator?

Thank you,

   Daniel