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

[Yi Hu]

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

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
-------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20240108/b3dd41e6/attachment-0001.html>