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