[petsc-users] Insert values into matrix using MatSetValuesStencil or MatSetValuesLocal
Wee Beng Tay
zonexo at gmail.com
Thu Aug 27 10:45:12 CDT 2015
Hi Timothee,
That's a better way. Thanks
Sent using CloudMagic [https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2] On Wed, Aug 26, 2015 at 1:15 PM, Timothée Nicolas < timothee.nicolas at gmail.com [timothee.nicolas at gmail.com] > wrote:
I don't really understand what you say, but it does not sound right.
You can enter the boundary points separately and then the points outside the
boundary on separate calls, like this :
do j=user%ys,user%ye
do i=user%xs,user%xe
if (i.eq.1 .or. i.eq.user%mx .or. j .eq. 1 .or. j .eq. user%my) then
! boundary point
row(MatStencil_i,1) = i -1
row(MatStencil_j,1) = j -1
row(MatStencil_c,1) = 1 -1
col(MatStencil_i,1) = i -1
col(MatStencil_j,1) = j -1
col(MatStencil_c,1) = 1 -1
v(1) = one
call MatSetValuesStencil(jac_prec,ione,row,ione,col,v, &
& INSERT_VALUES,ierr)
else
row(MatStencil_i,1) = i -1
row(MatStencil_j,1) = j -1
row(MatStencil_c,1) = 1 -1
col(MatStencil_i,1) = i -1
col(MatStencil_j,1) = j -1
col(MatStencil_c,1) = 1 -1
v(1) = undemi*dxm1*(vx_ip1j-vx_im1j) + two*user%nu*(dxm1**2+dym1**2)
col(MatStencil_i,2) = i+1 -1
col(MatStencil_j,2) = j -1
col(MatStencil_c,2) = 1 -1
v(2) = undemi*dxm1*(vx_ij-vx_ip1j) - user%nu*dxm1**2
col(MatStencil_i,3) = i-1 -1
col(MatStencil_j,3) = j -1
col(MatStencil_c,3) = 1 -1
v(3) = -undemi*dxm1*(vx_ij-vx_im1j) - user%nu*dxm1**2
col(MatStencil_i,4) = i -1
col(MatStencil_j,4) = j+1 -1
col(MatStencil_c,4) = 1 -1
v(4) = undemi*dym1*vy_ij - user%nu*dym1**2
col(MatStencil_i,5) = i -1
col(MatStencil_j,5) = j-1 -1
col(MatStencil_c,5) = 1 -1
v(5) = -undemi*dym1*vy_ij - user%nu*dym1**2
col(MatStencil_i,6) = i -1
col(MatStencil_j,6) = j -1
col(MatStencil_c,6) = 2 -1
v(6) = undemi*dym1*(vx_ijp1-vx_ijm1)
col(MatStencil_i,7) = i+1 -1
col(MatStencil_j,7) = j -1
col(MatStencil_c,7) = 2 -1
v(7) = -undemi*dxm1*vy_ip1j
col(MatStencil_i,8) = i-1 -1
col(MatStencil_j,8) = j -1
col(MatStencil_c,8) = 2 -1
v(8) = undemi*dxm1*vy_im1j
call MatSetValuesStencil(jac_prec,ione,row,ieight,col,v, &
& INSERT_VALUES,ierr)
endif
enddo
enddo
Timothee
2015-08-26 14:02 GMT+09:00 TAY wee-beng < zonexo at gmail.com [zonexo at gmail.com] > :
Hi Timothee,
Yes, I only parallelized in j and k. ksta,jsta are the starting k and j values.
kend,jend are the ending k and j values.
However, now I am using only 1 procs.
I was going to resend you my code but then I realised my mistake. I used:
call
MatSetValuesStencil(A_mat,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)
for all pts, including those at the boundary. Hence, those coupling outside the
boundary is also included.
I changed to:
call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,7),value_insert(7),INSERT_VALUES,ierr)
so I am now entering values individually.
Is there anyway I can use the 1st option to enter all the values together even
those some pts are invalid. I think it should be faster. Can I somehow tell
PETSc to ignore them?
Thank you
Yours sincerely,
TAY wee-beng
On 26/8/2015 12:24 PM, Timothée Nicolas wrote:
What is the definition of ksta, kend, jsta, jend ? Etc ? You are parallelized
only in j and k ?
What I said about the "-1" holds only if you have translated the start and end
points to FORTRAN numbering after getting the corners and ghost corners from the
DMDA (see ex ex5f90.F from snes)
Would you mind sending the complete routine with the complete definitions of
ksta,kend,jsta,jend,and size_x ?
Timothee
2015-08-26 13:12 GMT+09:00 TAY wee-beng < zonexo at gmail.com [zonexo at gmail.com] > :
Hi,
I have wrote the routine for my Poisson eqn. I have only 1 DOF, which is for
pressure. The center cell is coupled with 6 other cells (north, south, east,
west, front, back), so together 7 couplings.
size x/y/z = 4/8/10
MatStencil :: row(4,1),col(4,7)
PetscScalar :: value_insert(7)
PetscInt :: ione,iseven
ione = 1; iseven = 7
do k=ksta,kend
do j = jsta,jend
do i=1,size_x
row(MatStencil_i,1) = i - 1
row(MatStencil_j,1) = j - 1
row(MatStencil_k,1) = k - 1
row(MatStencil_c,1) = 0 ! 1 - 1
value_insert = 0.d0
if (i /= size_x) then
col(MatStencil_i,3) = i + 1 - 1 !east
col(MatStencil_j,3) = j - 1
col(MatStencil_k,3) = k - 1
col(MatStencil_c,3) = 0
value_insert(3) = (cp_yz(j,k)%fc_E)/(cp_x(i)%pd_E+cp_x(i+1)%pd_W)
end if
if (i /= 1) then
col(MatStencil_i,5) = i - 1 - 1 !west
col(MatStencil_j,5) = j - 1
col(MatStencil_k,5) = k - 1
col(MatStencil_c,5) = 0
value_insert(5) = (cp_yz(j,k)%fc_E)/(cp_x(i)%pd_W+cp_x(i-1)%pd_E)
end if
if (j /= size_y) then
col(MatStencil_i,2) = i - 1 !north
col(MatStencil_j,2) = j + 1 - 1
col(MatStencil_k,2) = k - 1
col(MatStencil_c,2) = 0
value_insert(2) = (cp_zx(i,k)%fc_N)/(cp_y(j)%pd_N+cp_y(j+1)%pd_S)
end if
...
col(MatStencil_i,1) = i - 1
col(MatStencil_j,1) = j - 1
col(MatStencil_k,1) = k - 1
col(MatStencil_c,1) = 0
value_insert(1) = -value_insert(2) - value_insert(3) - value_insert(4) -
value_insert(5) - value_insert(6) - value_insert(7)
call
MatSetValuesStencil(A_mat,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)
end do
end do
end do
but I got the error :
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at (3,0) in the matrix.
The error happens at i = 4, j = 1, k = 1. So I guess it has something to do with
the boundary condition. However, I can't figure out what's wrong. Can someone
help?
Thank you
Yours sincerely,
TAY wee-beng
On 24/8/2015 5:54 PM, Timothée Nicolas wrote:
Hi,
ex5 of snes can give you an example of the two routines.
The C version ex5.c uses MatSetValuesStencil whereas the Fortran90 version
ex5f90.F uses MatSetValuesLocal.
However, I use MatSetValuesStencil also in Fortran, there is no problem, and no
need to mess around with DMDAGetAO, I think.
To input values in the matrix, you need to do the following :
! Declare the matstencils for matrix columns and rows
MatStencil :: row(4,1),col(4,n)
! Declare the quantity which will store the actual matrix elements
PetscScalar :: v(8)
The first dimension in row and col is 4 to allow for 3 spatial dimensions (even
if you use only 2) plus one degree of freedom if you have several fields in your
DMDA. The second dimension is 1 for row (you input one row at a time) and n for
col, where n is the number of columns that you input. For instance, if at node
(1,i,j) (1 is the index of the degree of freedom), you have, say, 6 couplings,
with nodes (1,i,j), (1,i+1,j), (1,i-1,j), (1,i,j-1), (1,i,j+1), (2,i,j) for
example, then you need to set n=6
Then you define the row number by naturally doing the following, inside a local
loop :
row(MatStencil_i,1) = i -1
row(MatStencil_j,1) = j -1
row(MatStencil_c,1) = 1 -1
the -1 are here because FORTRAN indexing is different from the native C
indexing. I put them on the right to make this more apparent.
Then the column information. For instance to declare the coupling with node
(1,i,j), (1,i-1,j) and (2,i,j) (you can make up for the rest) you will have to
write (still within the same local loop on i and j)
col(MatStencil_i,1) = i -1
col(MatStencil_j,1) = j -1
col(MatStencil_c,1) = 1 -1
v(1) = whatever_it_is
col(MatStencil_i,2) = i-1 -1
col(MatStencil_j,2) = j -1
col(MatStencil_c,2) = 1 -1
v(2) = whatever_it_is
col(MatStencil_i,3) = i -1
col(MatStencil_j,3) = j -1
col(MatStencil_c,3) = 2 -1
v(3) = whatever_it_is
...
...
..
...
...
...
Note that the index of the degree of freedom (or what field you are coupling
to), is indicated by MatStencil_c
Finally use MatSetValuesStencil
ione = 1
isix = 6
call MatSetValuesStencil(Matrix,ione,row,isix,col,v,INSERT_VALUES,ierr)
If it is not clear don't hesitate to ask more details. For me it worked that
way, I succesfully computed a Jacobian that way. It is very sensitive. If you
slightly depart from the right jacobian, you will see a huge difference compared
to using matrix free with -snes_mf, so you can hardly make a mistake because you
would see it. That's how I finally got it to work.
Best
Timothee
2015-08-24 18:09 GMT+09:00 Wee-Beng Tay < [zonexo at gmail.com] zonexo at gmail.com [zonexo at gmail.com] > :
Hi,
I'm modifying my 3d fortran code from MPI along 1 direction (z) to MPI along 2
directions (y,z)
Previously I was using MatSetValues with global indices. However, now I'm using
DM and global indices is much more difficult.
I come across MatSetValuesStencil or MatSetValuesLocal.
So what's the difference bet the one since they both seem to work locally?
Which is a simpler/better option?
Is there an example in Fortran for MatSetValuesStencil?
Do I also need to use DMDAGetAO together with MatSetValuesStencil or
MatSetValuesLocal?
Thanks!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150827/ccf39e86/attachment-0001.html>
More information about the petsc-users
mailing list