The following assumes you are not using the shell matrix context for some other purpose
> subroutine formJacobian(snes,F_global,Jac,Jac_pre,dummy,err_PETSc) > > SNES :: snes > Vec :: F_global > > ! real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: & > ! F > Mat :: Jac, Jac_pre > PetscObject :: dummy > PetscErrorCode :: err_PETSc > > print*, '@@ start build my jac' > > PetscCall(MatShellSetContext(Jac,F_global,ierr)) ! record the current > base vector where the Jacobian is to be applied > print*, '@@ end build my jac' > > end subroutine formJacobian subroutine Gk_op ... Vec base PetscCall(MatShellGetContext(Jac,base,ierr)) ! use base in the computation of your matrix-free Jacobian vector product .... > On Jan 11, 2024, at 5:55 AM, Yi Hu <y...@mpie.de> wrote: > > Now I understand a bit more about the workflow of set jacobian. It seems that > the SNES can be really fine-grained. As you point out, J is built via > formJacobian() callback, and can be based on previous solution (or the base > vector u, as you mentioned). And then KSP can use a customized MATOP_MULT to > solve the linear equations J(u)*x=rhs. > > So I followed your idea about removing DMSNESSetJacobianLocal() and did the > following. > > …… > call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,& > > 9*product(cells(1:2))*cells3,9*product(cells(1:2))*cells3,& > 0,Jac_PETSc,err_PETSc) > CHKERRQ(err_PETSc) > call MatShellSetOperation(Jac_PETSc,MATOP_MULT,GK_op,err_PETSc) > CHKERRQ(err_PETSc) > call SNESSetJacobian(SNES_mech,Jac_PETSc,Jac_PETSc,formJacobian,0,err_PETSc) > CHKERRQ(err_PETSc) > …… > > And my formJacobian() is > > subroutine formJacobian(snes,F_global,Jac,Jac_pre,dummy,err_PETSc) > > SNES :: snes > Vec :: F_global > > ! real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: & > ! F > Mat :: Jac, Jac_pre > PetscObject :: dummy > PetscErrorCode :: err_PETSc > > print*, '@@ start build my jac' > > call MatCopy(Jac_PETSc,Jac,SAME_NONZERO_PATTERN,err_PETSc) > CHKERRQ(err_PETSc) > call MatCopy(Jac_PETSc,Jac_pre,SAME_NONZERO_PATTERN,err_PETSc) > CHKERRQ(err_PETSc) > ! Jac = Jac_PETSc > ! Jac_pre = Jac_PETSc > > print*, '@@ end build my jac' > > end subroutine formJacobian > > it turns out that no matter by a simple assignment or MatCopy(), the compiled > program gives me the same error as before. So I guess the real jacobian is > still not set. I wonder how to get around this and let this built jac in > formJacobian() to be the same as my shell matrix. > > Yi > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Wednesday, January 10, 2024 4:27 PM > To: Yi Hu <y...@mpie.de <mailto:y...@mpie.de>> > Cc: petsc-users <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> > Subject: Re: [petsc-users] SNES seems not use my matrix-free operation > > > By default if SNESSetJacobian() is not called with a function pointer PETSc > attempts to compute the Jacobian matrix explicitly with finite differences > and coloring. This doesn't makes sense with a shell matrix. Hence the error > message below regarding MatFDColoringCreate(). > > DMSNESSetJacobianLocal() calls SNESSetJacobian() with a function pointer of > SNESComputeJacobian_DMLocal() so preventing the error from triggering in your > code. > > You can provide your own function to SNESSetJacobian() and thus not need to > call DMSNESSetJacobianLocal(). What you do depends on how you want to record > the "base" vector that tells your matrix-free multiply routine where the > Jacobian matrix vector product is being applied, that is J(u)*x. u is the > "base" vector which is passed to the function provided with SNESSetJacobian(). > > Barry > > > > On Jan 10, 2024, at 6:20 AM, Yi Hu <y...@mpie.de <mailto:y...@mpie.de>> wrote: > > Thanks for the clarification. It is more clear to me now about the global to > local processes after checking the examples, e.g. > ksp/ksp/tutorials/ex14f.F90. > > And for using Vec locally, I followed your advice of VecGet.. and VecRestore… > In fact I used DMDAVecGetArrayReadF90() and some other relevant subroutines. > > For your comment on DMSNESSetJacobianLocal(). It seems that I need to use > both SNESSetJacobian() and then DMSNESSetJacobianLocal() to get things > working. When I do only SNESSetJacobian(), it does not work, meaning the > following does not work > > …… > call > DMDASNESsetFunctionLocal(DM_mech,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) > > CHKERRQ(err_PETSc) > call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,& > > 9*product(cells(1:2))*cells3,9*product(cells(1:2))*cells3,& > 0,Jac_PETSc,err_PETSc) > CHKERRQ(err_PETSc) > call MatShellSetOperation(Jac_PETSc,MATOP_MULT,GK_op,err_PETSc) > CHKERRQ(err_PETSc) > call > SNESSetJacobian(SNES_mech,Jac_PETSc,Jac_PETSc,PETSC_NULL_FUNCTION,0,err_PETSc) > CHKERRQ(err_PETSc) > !call DMSNESsetJacobianLocal(DM_mech,formJacobian,PETSC_NULL_SNES,err_PETSc) > !CHKERRQ(err_PETSc) > call > SNESsetConvergenceTest(SNES_mech,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) > > CHKERRQ(err_PETSc) > call SNESSetDM(SNES_mech,DM_mech,err_PETSc) > CHKERRQ(err_PETSc) > …… > > It gives me the message > [0]PETSC ERROR: No support for this operation for this object type > > [0]PETSC ERROR: Code not yet written for matrix type shell > > [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting. > > [0]PETSC ERROR: Petsc Release Version 3.16.4, Feb 02, 2022 > [0]PETSC ERROR: Configure options PETSC_ARCH=linux-gnu > --with-fortran-bindings --with-mpi-f90module-visibility=0 --download-fftw > --download-hdf5 --download-hdf5-fortran-bindings --download-fblaslapack > --download-ml --download-zlib > > [0]PETSC ERROR: #1 MatFDColoringCreate() at > /home/yi/app/petsc-3.16.4/src/mat/matfd/fdmatrix.c:471 > [0]PETSC ERROR: #2 SNESComputeJacobian_DMDA() at > /home/yi/app/petsc-3.16.4/src/snes/utils/dmdasnes.c:173[0]PETSC ERROR: #3 > SNESComputeJacobian() at > /home/yi/app/petsc-3.16.4/src/snes/interface/snes.c:2864 > [0]PETSC ERROR: #4 SNESSolve_NEWTONLS() at > /home/yi/app/petsc-3.16.4/src/snes/impls/ls/ls.c:222 > [0]PETSC ERROR: #5 SNESSolve() at > /home/yi/app/petsc-3.16.4/src/snes/interface/snes.c:4809 > [0]PETSC ERROR: #6 User provided function() at User file:0 > > [0]PETSC ERROR: #7 VecSetErrorIfLocked() at > /home/yi/app/petsc-3.16.4/include/petscvec.h:623 > [0]PETSC ERROR: #8 VecGetArray() at > /home/yi/app/petsc-3.16.4/src/vec/vec/interface/rvector.c:1769 > [0]PETSC ERROR: #9 User provided function() at User file:0 > > [0]PETSC ERROR: #10 MatFDColoringCreate() at > /home/yi/app/petsc-3.16.4/src/mat/matfd/fdmatrix.c:471 > [0]PETSC ERROR: #11 SNESComputeJacobian_DMDA() at > /home/yi/app/petsc-3.16.4/src/snes/utils/dmdasnes.c:173 > > > [0]PETSC ERROR: #12 SNESComputeJacobian() at > /home/yi/app/petsc-3.16.4/src/snes/interface/snes.c:2864 > [0]PETSC ERROR: #13 SNESSolve_NEWTONLS() at > /home/yi/app/petsc-3.16.4/src/snes/impls/ls/ls.c:222 > [0]PETSC ERROR: #14 SNESSolve() at > /home/yi/app/petsc-3.16.4/src/snes/interface/snes.c:4809 > > It seems that I have to use a DMSNESSetJacobianLocal() to “activate” the use > of my shell matrix, although the formJacobian() in the > DMSNESSetJacobianLocal() is doing nothing. > > Best wishes, > Yi > > > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Tuesday, January 9, 2024 4:49 PM > To: Yi Hu <y...@mpie.de <mailto:y...@mpie.de>> > Cc: petsc-users <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> > Subject: Re: [petsc-users] SNES seems not use my matrix-free operation > > However, my GK_op (which is the reloaded MATOP_MULT) gives me some problem. > It is entered but then crashed with Segmentation Violation error. So I guess > my indices may be wrong. I wonder do I need to use the local Vec (of dF), and > should my output Vec also be in the correct shape (i.e. after calculation I > need to transform back into a Vec)? As you can see here, my dF is a tensor > defined on every grid point. > > > The input for the matrix-vector product is a global vector, as is the > result. (Not like the arguments to DMSNESSetJacobianLocal). > > This means that your MATOP_MULT function needs to do the > DMGlobalToLocal() vector operation first then the "unwrapping" from the > vector to the array format at the beginning of the routine. Similarly it > needs to "unwrap" the result vector as an array. See > src/snes/tutorials/ex14f.F90 and in particular the code block > > PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)) > PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)) > > ! Get pointers to vector data > > PetscCall(VecGetArrayReadF90(localX,xx,ierr)) > PetscCall(VecGetArrayF90(F,ff,ierr)) > > Barry > > You really shouldn't be using DMSNESSetJacobianLocal() for your code. > Basically all the DMSNESSetJacobianLocal() gives you is that it automatically > handles the global to local mapping and unwrapping of the vector to an array, > but it doesn't work for shell matrices. > > > > > On Jan 9, 2024, at 6:30 AM, Yi Hu <y...@mpie.de <mailto:y...@mpie.de>> wrote: > > Dear Barry, > > Thanks for your help. > > It works when doing first SNESSetJacobian() with my created shell matrix Jac > in the main (or module) and then DMSNESSetJacobianLocal() to associate with > my DM and an dummy formJacobian callback (which is doing nothing). My SNES > can now recognize my shell matrix and do my customized operation. > > However, my GK_op (which is the reloaded MATOP_MULT) gives me some problem. > It is entered but then crashed with Segmentation Violation error. So I guess > my indices may be wrong. I wonder do I need to use the local Vec (of dF), and > should my output Vec also be in the correct shape (i.e. after calculation I > need to transform back into a Vec)? As you can see here, my dF is a tensor > defined on every grid point. > > Best wishes, > Yi > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Monday, January 8, 2024 6:41 PM > To: Yi Hu <y...@mpie.de <mailto:y...@mpie.de>> > Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> > Subject: Re: [petsc-users] SNES seems not use my matrix-free operation > > > "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 <mailto: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> > ------------------------------------------------- > > > > > ------------------------------------------------- > 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> > ------------------------------------------------- > > > > ------------------------------------------------- > 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> > ------------------------------------------------- > > > > ------------------------------------------------- > 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> > -------------------------------------------------
