Dear Barry,
Thanks for your reply. In fact, this still does not work, giving the error
“code not yet written for matrix type shell”. I also tried recording F_global
as base in my global shell matrix Jac_PETSc. And fetch it in my GK_op, but the
same error occurs. So the question is still, how can I assign the Jac in my
formJacobian() to my global defined Jac_PETSc when set up with
SNESSetJacobian(snes,Jac_PETSc,Jac_PETSc,formJacobian,0,ierr). As far as I
understand, when using SNESSetJacobian() the jac argument is only for reserving
memory. And the matrix is not written yet. Only a valid formJacobian() can lead
to a write to the matrix. My formJacobian maybe not valid (because the output
not gives me print), so the Jac_PETSc is never written for the SNES. Maybe I am
wrong about it.
Best wishes,
Yi
From: Barry Smith <bsm...@petsc.dev>
Sent: Thursday, January 11, 2024 4:47 PM
To: Yi Hu <y...@mpie.de>
Cc: petsc-users <petsc-users@mcs.anl.gov>
Subject: Re: [petsc-users] SNES seems not use my matrix-free operation
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>
Sent: Wednesday, January 10, 2024 4:27 PM
To: Yi Hu <y...@mpie.de>
Cc: petsc-users <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> 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>
Sent: Tuesday, January 9, 2024 4:49 PM
To: Yi Hu <y...@mpie.de>
Cc: petsc-users <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> 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>
Sent: Monday, January 8, 2024 6:41 PM
To: Yi Hu <y...@mpie.de>
Cc: 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> 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
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
-------------------------------------------------
-------------------------------------------------
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
-------------------------------------------------
-------------------------------------------------
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
-------------------------------------------------
-------------------------------------------------
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
-------------------------------------------------
-------------------------------------------------
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
-------------------------------------------------
