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

Mon Nov 4 04:23:35 CST 2024

Dear all,

The problem was the context. The context was not properly set, and yet 
for some reason it was running correctly in Debug mode.

Now this problem has been solved.

Thank you Barry for your help.

Best,

   Daniel

On 29/10/2024 23:04, Barry Smith wrote:
>    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!dkbZ6vbH8gLYOkcINufnJg_JyJcRHtNge8M_c8f-iuQ3J8EP-KVV58zSZ1cBNyyStNlGaXrJUcwnob5EdBnPDpfAiOKcpEClnwS21YZuJw$  and https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dkbZ6vbH8gLYOkcINufnJg_JyJcRHtNge8M_c8f-iuQ3J8EP-KVV58zSZ1cBNyyStNlGaXrJUcwnob5EdBnPDpfAiOKcpEClnwTJP97sJQ$ 
>>>> [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
>>>>>>>>