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

[Weide, Edwin van der (UT-ET)]

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