[petsc-users] Matrix free SNES with user provided matrix vector product and preconditioner operation

[Barry Smith]

   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/4f719e7d/attachment-0001.html>