Yes, this works. Thanks a lot for your help. Regards, Edwin
________________________________ From: Barry Smith <bsm...@petsc.dev> Sent: Wednesday, November 6, 2024 7:03 PM To: Weide, Edwin van der (UT-ET) <e.t.a.vanderwe...@utwente.nl> Cc: petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov> Subject: Re: [petsc-users] Matrix free SNES with user provided matrix vector product and preconditioner operation You need to provide a callback function. Why? Otherwise your MatShell and PCshell have no way of knowing at what location the Jacobian is suppose to be evaluated at (in a matrix free way). That is the x for which J(x) is used. Normally one puts the x into the application context of mJac and accesses it every time the matmult is called. Similarly it needs to be accessed in application of your preconditioner. On Nov 6, 2024, at 12:01 PM, Weide, Edwin van der (UT-ET) <e.t.a.vanderwe...@utwente.nl> wrote: Barry, If I do that, I get the following error [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!dfJ_L4Ay2QlfGz0JU6uAfD7Zy_6a7x9D0bisHMG1T53fTy_gWL3Q8fzzLvWZ_TBK37FFmDYerLywkFOLtMnU1rISfF_VA3gxsDTm$ and https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dfJ_L4Ay2QlfGz0JU6uAfD7Zy_6a7x9D0bisHMG1T53fTy_gWL3Q8fzzLvWZ_TBK37FFmDYerLywkFOLtMnU1rISfF_VAwRW4W44$ [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 Jacobian [0]PETSC ERROR: #2 SNESComputeJacobian() at /home/vdweide/petsc/src/snes/interface/snes.c:2966 [0]PETSC ERROR: #3 SNESSolve_NEWTONLS() at /home/vdweide/petsc/src/snes/impls/ls/ls.c:218 [0]PETSC ERROR: #4 SNESSolve() at /home/vdweide/petsc/src/snes/interface/snes.c:4841 [0]PETSC ERROR: #5 SolveCurrentStage() at SolverClass.cpp:502 [0]PETSC ERROR: #6 main() at Condensation.cpp:20 -------------------------------------------------------------------------- So SNES tries to call the call back function for the Jacobian, but that is not provided. Hence the failure. Regards, Edwin ________________________________ From: Barry Smith <bsm...@petsc.dev<mailto:bsm...@petsc.dev>> Sent: Wednesday, November 6, 2024 5:52 PM To: Weide, Edwin van der (UT-ET) <e.t.a.vanderwe...@utwente.nl<mailto:e.t.a.vanderwe...@utwente.nl>> Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> <petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>> Subject: Re: [petsc-users] Matrix free SNES with user provided matrix vector product and preconditioner operation Just pass mJac, mJac instead of mJac, nullptr and it will be happy. In your case, the second mJac won't be used in your preconditioner it is just a place holder so other parts of SNES won't try to create a matrix. Barry On Nov 6, 2024, at 11:36 AM, Weide, Edwin van der (UT-ET) via petsc-users <petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>> wrote: Hi, I am trying to solve a nonlinear problem with matrix-free SNES where I would like to provide both the matrix vector product and the preconditioner myself. For that purpose, I use the following construction. // Set up the matrix free evaluation of the Jacobian times a vector // by setting the appropriate function in snes. PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, nEqns, nEqns, this, &mJac)); PetscCall(MatShellSetOperation(mJac, MATOP_MULT, (void (*)(void))JacobianTimesVector)); PetscCall(SNESSetJacobian(mSnes, mJac, nullptr, nullptr, nullptr)); // Set the function to be used as preconditioner for the krylov solver. KSP ksp; PC pc; PetscCall(SNESGetKSP(mSnes, &ksp)); PetscCall(KSPGetPC(ksp, &pc)); PetscCall(PCSetType(pc, PCSHELL)); PetscCall(PCSetApplicationContext(pc, this)); PetscCall(PCShellSetApply(pc, Preconditioner)); For small problems this construction works, and it does exactly what I expect it to do. However, when I increase the problem size, I get a memory allocation failure in SNESSolve, because it looks like SNES attempts to allocate memory for a full dense matrix for the preconditioner, which is not used. This is the call stack when the error occurs. [0]PETSC ERROR: #1 PetscMallocAlign() at /home/vdweide/petsc/src/sys/memory/mal.c:53 [0]PETSC ERROR: #2 PetscTrMallocDefault() at /home/vdweide/petsc/src/sys/memory/mtr.c:175 [0]PETSC ERROR: #3 PetscMallocA() at /home/vdweide/petsc/src/sys/memory/mal.c:421 [0]PETSC ERROR: #4 MatSeqDenseSetPreallocation_SeqDense() at /home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:3357 [0]PETSC ERROR: #5 MatSeqDenseSetPreallocation() at /home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:3338 [0]PETSC ERROR: #6 MatDuplicateNoCreate_SeqDense() at /home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:372 [0]PETSC ERROR: #7 MatDuplicate_SeqDense() at /home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:399 [0]PETSC ERROR: #8 MatDuplicate() at /home/vdweide/petsc/src/mat/interface/matrix.c:4964 [0]PETSC ERROR: #9 DMCreateMatrix_Shell() at /home/vdweide/petsc/src/dm/impls/shell/dmshell.c:195 [0]PETSC ERROR: #10 DMCreateMatrix() at /home/vdweide/petsc/src/dm/interface/dm.c:1501 [0]PETSC ERROR: #11 SNESSetUpMatrices() at /home/vdweide/petsc/src/snes/interface/snes.c:794 [0]PETSC ERROR: #12 SNESSetUp_NEWTONLS() at /home/vdweide/petsc/src/snes/impls/ls/ls.c:290 [0]PETSC ERROR: #13 SNESSetUp() at /home/vdweide/petsc/src/snes/interface/snes.c:3395 [0]PETSC ERROR: #14 SNESSolve() at /home/vdweide/petsc/src/snes/interface/snes.c:4831 [0]PETSC ERROR: #15 SolveCurrentStage() at SolverClass.cpp:502 In the function SNESSetUpMatrices the source looks as follows 784 } else if (!snes->jacobian_pre) { 785 PetscDS prob; 786 Mat J, B; 787 PetscBool hasPrec = PETSC_FALSE; 788 789 J = snes->jacobian; 790 PetscCall(DMGetDS(dm, &prob)); 791 if (prob) PetscCall(PetscDSHasJacobianPreconditioner(prob, &hasPrec)); 792 if (J) PetscCall(PetscObjectReference((PetscObject)J)); 793 else if (hasPrec) PetscCall(DMCreateMatrix(snes->dm, &J)); 794 PetscCall(DMCreateMatrix(snes->dm, &B)); 795 PetscCall(SNESSetJacobian(snes, J ? J : B, B, NULL, NULL)); 796 PetscCall(MatDestroy(&J)); 797 PetscCall(MatDestroy(&B)); 798 } It looks like in line 794 it is attempted to create the preconditioner, because it was (intentionally) not provided. Hence my question. Is it possible to use matrix-free SNES with a user provided matrix vector product (via MatShell) and a user provided preconditioner operation without SNES allocating the memory for a dense matrix? If so, what do I need to change in the construction above to make it work? If needed, I can provide the source code for which this problem occurs. Thanks, Edwin --------------------------------------------------- Edwin van der Weide Department of Mechanical Engineering University of Twente Enschede, the Netherlands
