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 r...@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 r...@mpie.de -------------------------------------------------
