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!elK6sfqwXIIaLqDzUlKw8Wy--UjvmZXYI0gRqu_nvUNMbIx2rMEX4aHIYXWikS4p4zHPTXyicwT9SEY8-JQt10s$ > and > https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!elK6sfqwXIIaLqDzUlKw8Wy--UjvmZXYI0gRqu_nvUNMbIx2rMEX4aHIYXWikS4p4zHPTXyicwT9SEY8mXmEtXk$ > > [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