[petsc-users] SNES seems not use my matrix-free operation

[Barry Smith]

   "formJacobian" should not be __creating__ the matrices. Here "form" means computing the numerical values in the matrix (or when using a shell matrix it means keeping a copy of X so that your custom matrix-free multiply knows the base location where the matrix free Jacobian-vector products are computed.)

   You create the shell matrices up in your main program and pass them in with SNESSetJacobian(). 

    Try first calling SNESSetJacobian() to provide the matrices (provide a dummy function argument) and then call DMSNESSetJacobianLocal() to provide your "formjacobian"  function (that does not create the matrices).

   Barry


   Yes, "form" is a bad word that should not have been used in our code.



> On Jan 8, 2024, at 12:24 PM, Yi Hu <y.hu at mpie.de> wrote:
> 
> Dear PETSc Experts,
>  
> I am implementing a matrix-free jacobian for my SNES solver in Fortran. (command line option -snes_type newtonls -ksp_type gmres)
>  
> In the main program, I define my residual and jacobian and matrix-free jacobian like the following,
>  
> …
> call DMDASNESSetFunctionLocal(DM_mech, INSERT_VALUES, formResidual, PETSC_NULL_SNES, err_PETSc)
> call DMSNESSetJacobianLocal(DM_mech, formJacobian, PETSC_NULL_SNES, err_PETSc)
> …
>  
> subroutine formJacobian(residual_subdomain,F,Jac_pre,Jac,dummy,err_PETSc)
>   
> #include <petsc/finclude/petscmat.h>
>   use petscmat
>   implicit None
>   DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
>     residual_subdomain                                                                              !< DMDA info (needs to be named "in" for macros like XRANGE to work)
>   real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: &
>     F                                                                                               !< deformation gradient field
>   Mat                                  :: Jac, Jac_pre
>   PetscObject                          :: dummy
>   PetscErrorCode                       :: err_PETSc
>   PetscInt                             :: N_dof ! global number of DoF, maybe only a placeholder
>  
>   N_dof = 9*product(cells(1:2))*cells3 
>  
>   print*, 'in my jac'
>   
>   call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac,err_PETSc)
>   CHKERRQ(err_PETSc)
>   call MatShellSetOperation(Jac,MATOP_MULT,GK_op,err_PETSc)
>   CHKERRQ(err_PETSc)
>   
>   print*, 'in my jac'
>  
>   ! for jac preconditioner
>   call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac_pre,err_PETSc)
>   CHKERRQ(err_PETSc)
>   call MatShellSetOperation(Jac_pre,MATOP_MULT,GK_op,err_PETSc)
>   CHKERRQ(err_PETSc)
>  
>   print*, 'in my jac'
>  
> end subroutine formJacobian
>  
> subroutine GK_op(Jac,dF,output,err_PETSc)
>   real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: &
>     dF                                                                                               
>   real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(out) :: &
>     output                                                                                               
>   real(pREAL),  dimension(3,3) :: &
>     deltaF_aim = 0.0_pREAL
>  
>   Mat                                  :: Jac
>   PetscErrorCode                       :: err_PETSc
>  
>   integer :: i, j, k, e
>  
>   … a lot of calculations …
>  
>   print*, 'in GK op'
>   
> end subroutine GK_op
>  
> The first question is that: it seems I still need to explicitly define the interface of MatCreateShell() and MatShellSetOperation() to properly use them, even though I include them via “use petscmat”. It is a little bit strange to me, since some examples do not perform this step. 
>  
> Then the main issue is that I can build my own Jacobian from my call back function formJacobian, and confirm my Jacobian is a shell matrix (by MatView). However, my customized operator GK_op is not called when solving the nonlinear system (not print my “in GK op”). When I try to monitor my SNES, it gives me some conventional output not mentioning my matrix-free operations. So I guess my customized MATOP_MULT may be not associated with Jacobian. Or my configuration is somehow wrong. Could you help me solve this issue? 
>  
> Thanks,
> Yi
>  
> 
> 
> -------------------------------------------------
> Stay up to date and follow us on LinkedIn, Twitter and YouTube.
> 
> Max-Planck-Institut für Eisenforschung GmbH
> Max-Planck-Straße 1
> D-40237 Düsseldorf
>  
> Handelsregister B 2533 
> Amtsgericht Düsseldorf
>  
> Geschäftsführung
> Prof. Dr. Gerhard Dehm
> Prof. Dr. Jörg Neugebauer
> Prof. Dr. Dierk Raabe
> Dr. Kai de Weldige
>  
> Ust.-Id.-Nr.: DE 11 93 58 514 
> Steuernummer: 105 5891 1000
> 
> 
> Please consider that invitations and e-mails of our institute are 
> only valid if they end with …@mpie.de. 
> If you are not sure of the validity please contact rco at mpie.de <mailto:rco at mpie.de>
> 
> Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails
> aus unserem Haus nur mit der Endung …@mpie.de gültig sind. 
> In Zweifelsfällen wenden Sie sich bitte an rco at mpie.de <mailto:rco at mpie.de>
> -------------------------------------------------

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20240108/1bcb1cde/attachment-0001.html>