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

[Barry Smith]

  Run in the debugger, even though there will not be full debugging support with the optimizations it should still provide you some information


> On Oct 29, 2024, at 4:35 PM, Daniel Pino Munoz <daniel.pino_munoz at mines-paristech.fr> wrote:
> 
> That's what I thought, so I replaced its content by:
> 
> VecZeroEntries(f);
> 
> and the result is the same...
> 
> 
> On 29/10/2024 21:31, Barry Smith wrote:
>>   This
>>   [0]PETSC ERROR: #1 SNES callback function
>>    indicates the crash is in your computeResidual function and thus you need to debug your function
>> 
>>    Barry
>> 
>> 
>>> On Oct 29, 2024, at 4:28 PM, Daniel Pino Munoz <daniel.pino_munoz at mines-paristech.fr> wrote:
>>> 
>>> I ran it with -malloc_debug and it does not change anything.
>>> 
>>> The output is the following:
>>> 
>>> he absolute tolerance is 0.001
>>> The relative tolerance is 0.001
>>> The divergence tolerance is 10000
>>> The maximum iterations is 10000
>>> Initial load !
>>> [0]PETSC ERROR: ------------------------------------------------------------------------
>>> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
>>> [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
>>> [0]PETSC ERROR: or see https://urldefense.us/v3/__https://petsc.org/release/faq/*valgrind__;Iw!!G_uCfscf7eWS!dibDYt_ZFa-6FBchlzMndWYg1NqvsEIZOJlfOEQkpzMOqyVUFiiumr5JERY5e-jjXiKd8sI5HQxEi-65JQwqWxc$  and https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dibDYt_ZFa-6FBchlzMndWYg1NqvsEIZOJlfOEQkpzMOqyVUFiiumr5JERY5e-jjXiKd8sI5HQxEi-65GdSV720$ 
>>> [0]PETSC ERROR: ---------------------  Stack Frames ------------------------------------
>>> [0]PETSC ERROR: The line numbers in the error traceback are not always exact.
>>> [0]PETSC ERROR: #1 SNES callback function
>>> [0]PETSC ERROR: #2 SNESComputeFunction() at /home/daniel-pino/Software/Dependencies/petsc/src/snes/interface/snes.c:2489
>>> [0]PETSC ERROR: #3 SNESSolve_KSPONLY() at /home/daniel-pino/Software/Dependencies/petsc/src/snes/impls/ksponly/ksponly.c:27
>>> [0]PETSC ERROR: #4 SNESSolve() at /home/daniel-pino/Software/Dependencies/petsc/src/snes/interface/snes.c:4841
>>> --------------------------------------------------------------------------
>>> MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
>>>   Proc: [[27669,1],0]
>>>   Errorcode: 59
>>> 
>>> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
>>> You may or may not see output from other processes, depending on
>>> exactly when Open MPI kills them.
>>> --------------------------------------------------------------------------
>>> 
>>> 
>>> 
>>> On 29/10/2024 20:17, Barry Smith wrote:
>>>>    Hmm, cut and paste the output when it crashes.
>>>> 
>>>>     Also run with -malloc_debug and see what happens
>>>> 
>>>> 
>>>>     Barry
>>>> 
>>>> 
>>>>> On Oct 29, 2024, at 3:13 PM, Daniel Pino Munoz <daniel.pino_munoz at mines-paristech.fr> wrote:
>>>>> 
>>>>> 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_munoz at 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
>>>>>>>