[petsc-users] Matrix free SNES with user provided matrix vector product and preconditioner operation
Barry Smith
bsmith at petsc.dev
Wed Nov 6 12:03:47 CST 2024
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.vanderweide at 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 <bsmith at petsc.dev <mailto:bsmith at petsc.dev>>
> Sent: Wednesday, November 6, 2024 5:52 PM
> To: Weide, Edwin van der (UT-ET) <e.t.a.vanderweide at utwente.nl <mailto:e.t.a.vanderweide at utwente.nl>>
> Cc: petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov <mailto:petsc-users at 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 at mcs.anl.gov <mailto:petsc-users at 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20241106/505f3c66/attachment-0001.html>
More information about the petsc-users
mailing list