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
>
