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

Thu Nov 7 11:21:19 CST 2024

Yes, this works. Thanks a lot for your help.
Regards,

Edwin

________________________________
From: Barry Smith <bsmith at petsc.dev>
Sent: Wednesday, November 6, 2024 7:03 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

  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!dfJ_L4Ay2QlfGz0JU6uAfD7Zy_6a7x9D0bisHMG1T53fTy_gWL3Q8fzzLvWZ_TBK37FFmDYerLywkFOLtMnU1rISfF_VA3gxsDTm$  and https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!dfJ_L4Ay2QlfGz0JU6uAfD7Zy_6a7x9D0bisHMG1T53fTy_gWL3Q8fzzLvWZ_TBK37FFmDYerLywkFOLtMnU1rISfF_VAwRW4W44$ 
[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/20241107/8149a50c/attachment-0001.html>