[petsc-users] Matrix free SNES with user provided matrix vector product and preconditioner operation
Weide, Edwin van der (UT-ET)
e.t.a.vanderweide at utwente.nl
Wed Nov 6 11:01:48 CST 2024
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!ZZB4be6WT2eisTFLrrHknlukvw3s7eibKfF0XMaUNTO9axbKKdE3l3C_RbgSVoOHvWH97FEJKHr6mY18J8kMv23zN-9CL9kUhgje$ and https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!ZZB4be6WT2eisTFLrrHknlukvw3s7eibKfF0XMaUNTO9axbKKdE3l3C_RbgSVoOHvWH97FEJKHr6mY18J8kMv23zN-9CLwFG0b6V$
[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>
Sent: Wednesday, November 6, 2024 5:52 PM
To: Weide, Edwin van der (UT-ET) <e.t.a.vanderweide at utwente.nl>
Cc: petsc-users at mcs.anl.gov <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> 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/2de4c023/attachment-0001.html>
More information about the petsc-users
mailing list