[petsc-users] Insert values into matrix using MatSetValuesStencil or MatSetValuesLocal
Timothée Nicolas
timothee.nicolas at gmail.com
Tue Aug 25 23:24:18 CDT 2015
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>:
> 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>:
>
>> 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/20150826/ebe41779/attachment-0001.html>
More information about the petsc-users
mailing list