[petsc-users] Insert values into matrix using MatSetValuesStencil or MatSetValuesLocal

[Wee-Beng Tay]

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!

[Timothée Nicolas]

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

Post by Wee-Beng Tay
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!

[Timothée Nicolas]

Small erratum, in the declaration for v, it should be

PetscScalar :: v(n)

where n is the same as for col (6 in the example, not 8 which I copied from
my particular case)

Post by Timothée Nicolas
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.
! 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
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

Post by Wee-Beng Tay
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!

[Wee Beng Tay]

Sent using CloudMagic [https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2] On Mon, Aug 24, 2015 at 5:54 PM, Timothée Nicolas < ***@gmail.com [***@gmail.com] > 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


Hi Timothee,
Thanks for the help. So for boundary pts I will just leave blank for non
existent locations?
Also, can I use PETSc multigrid to solve this problem? This is a poisson eqn.
2015-08-24 18:09 GMT+09:00 Wee-Beng Tay < ***@gmail.com [***@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!

[TAY wee-beng]

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

Post by Timothée Nicolas
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.
! 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,
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
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!

[Timothée Nicolas]

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

Post by TAY wee-beng
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*
[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
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.
! 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
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

Post by Wee-Beng Tay
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!

[TAY wee-beng]

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

Post by Timothée Nicolas
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
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/*
[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

Post by Timothée Nicolas
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.
! 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,
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
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!

[Timothée Nicolas]

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

Post by Wee Beng Tay
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.
*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.
*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
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

Post by TAY wee-beng
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*
[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
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.
! 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
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

Post by Wee-Beng Tay
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!

[Wee Beng Tay]

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 < ***@gmail.com [***@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 < ***@gmail.com [***@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 < ***@gmail.com [***@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 < [***@gmail.com] ***@gmail.com [***@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!

[TAY wee-beng]

This post might be inappropriate. Click to display it.

[Timothée Nicolas]

The first thing that strikes me is your definition of the stencils

*MatStencil :: row(6,1),col(6,7)*

Why is it not defined with

*MatStencil :: row(4,1),col(4,7)*

instead ?

Where does the 6 come from ?

Timothée

Post by TAY wee-beng
Hi,
I have successfully used MatSetValuesStencil to insert values into a
Poisson eqn matrix which has 1 DOF (pressure). Now I'm trying to insert
values in a momentum eqn matrix which has 3 DOF (u,v,w)
*[0]PETSC ERROR: --------------------- Error Message
----------------------------*
*----------------------------------*
*[0]PETSC ERROR: Argument out of range*
*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in the matrix*
*[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
<http://www.mcs.anl.gov/petsc/documentation/faq.html> for trou*
*ble shooting.*
I wonder what's wrong. For the momentum eqn, for each DOF, at at node
(dof,i,j,k), I have coupling i +/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j = 2, k = 2,
which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*PetscScalar :: value_insert(7)*
*MatStencil :: row(6,1),col(6,7)*
*ione = 1; iseven = 7*
*if (cell_type == 'u') then*
* offset = 1*
*else if (cell_type == 'v') then*
* offset = 2*
*else if (cell_type == 'w') then*
* offset = 3*
*end if*
*do k=ksta2,kend2*
* do j = jsta2,jend2*
* do i=2,size_x-1*
* row(MatStencil_i,1) = i - 1*
* row(MatStencil_j,1) = j - 1*
* row(MatStencil_k,1) = k - 1*
* row(MatStencil_c,1) = offset - 1*
* value_insert = 0.d0*
* col(MatStencil_i,3) = i + 1 - 1 !east*
* col(MatStencil_j,3) = j - 1*
* col(MatStencil_k,3) = k - 1*
* col(MatStencil_c,3) = offset - 1*
* value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,5) = i - 1 - 1 !west*
* col(MatStencil_j,5) = j - 1*
* col(MatStencil_k,5) = k - 1*
* col(MatStencil_c,5) = offset - 1*
* value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,2) = i - 1 !north*
* col(MatStencil_j,2) = j + 1 - 1*
* col(MatStencil_k,2) = k - 1*
* col(MatStencil_c,2) = offset - 1*
* value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,4) = i - 1 !south*
* col(MatStencil_j,4) = j - 1 - 1*
* col(MatStencil_k,4) = k - 1*
* col(MatStencil_c,4) = offset - 1*
* value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,6) = i - 1 !front*
* col(MatStencil_j,6) = j - 1*
* col(MatStencil_k,6) = k + 1 - 1*
* col(MatStencil_c,6) = offset - 1*
* value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,7) = i - 1 !back*
* col(MatStencil_j,7) = j - 1*
* col(MatStencil_k,7) = k - 1 - 1*
* col(MatStencil_c,7) = offset - 1*
* value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,1) = i - 1*
* col(MatStencil_j,1) = j - 1*
* col(MatStencil_k,1) = k - 1*
* col(MatStencil_c,1) = offset - 1*
* value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))*
* call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)*
* end do*
* end do*
*end do *
Thank you
Yours sincerely,
TAY wee-beng

[Timothée Nicolas]

Can you also tell how you created the matrix ? Just in case you created it
with the 1 dof DMDA, it would not work if you try to input values at places
where it is not allocated (which could explain the error message)

Post by Timothée Nicolas
The first thing that strikes me is your definition of the stencils
*MatStencil :: row(6,1),col(6,7)*
Why is it not defined with
*MatStencil :: row(4,1),col(4,7)*
instead ?
Where does the 6 come from ?
Timothée

Post by TAY wee-beng
Hi,
I have successfully used MatSetValuesStencil to insert values into a
Poisson eqn matrix which has 1 DOF (pressure). Now I'm trying to insert
values in a momentum eqn matrix which has 3 DOF (u,v,w)
*[0]PETSC ERROR: --------------------- Error Message
----------------------------*
*----------------------------------*
*[0]PETSC ERROR: Argument out of range*
*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in the matrix*
*[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
<http://www.mcs.anl.gov/petsc/documentation/faq.html> for trou*
*ble shooting.*
I wonder what's wrong. For the momentum eqn, for each DOF, at at node
(dof,i,j,k), I have coupling i +/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j = 2, k = 2,
which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*PetscScalar :: value_insert(7)*
*MatStencil :: row(6,1),col(6,7)*
*ione = 1; iseven = 7*
*if (cell_type == 'u') then*
* offset = 1*
*else if (cell_type == 'v') then*
* offset = 2*
*else if (cell_type == 'w') then*
* offset = 3*
*end if*
*do k=ksta2,kend2*
* do j = jsta2,jend2*
* do i=2,size_x-1*
* row(MatStencil_i,1) = i - 1*
* row(MatStencil_j,1) = j - 1*
* row(MatStencil_k,1) = k - 1*
* row(MatStencil_c,1) = offset - 1*
* value_insert = 0.d0*
* col(MatStencil_i,3) = i + 1 - 1 !east*
* col(MatStencil_j,3) = j - 1*
* col(MatStencil_k,3) = k - 1*
* col(MatStencil_c,3) = offset - 1*
* value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,5) = i - 1 - 1 !west*
* col(MatStencil_j,5) = j - 1*
* col(MatStencil_k,5) = k - 1*
* col(MatStencil_c,5) = offset - 1*
* value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,2) = i - 1 !north*
* col(MatStencil_j,2) = j + 1 - 1*
* col(MatStencil_k,2) = k - 1*
* col(MatStencil_c,2) = offset - 1*
* value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,4) = i - 1 !south*
* col(MatStencil_j,4) = j - 1 - 1*
* col(MatStencil_k,4) = k - 1*
* col(MatStencil_c,4) = offset - 1*
* value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,6) = i - 1 !front*
* col(MatStencil_j,6) = j - 1*
* col(MatStencil_k,6) = k + 1 - 1*
* col(MatStencil_c,6) = offset - 1*
* value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,7) = i - 1 !back*
* col(MatStencil_j,7) = j - 1*
* col(MatStencil_k,7) = k - 1 - 1*
* col(MatStencil_c,7) = offset - 1*
* value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,1) = i - 1*
* col(MatStencil_j,1) = j - 1*
* col(MatStencil_k,1) = k - 1*
* col(MatStencil_c,1) = offset - 1*
* value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))*
* call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)*
* end do*
* end do*
*end do *
Thank you
Yours sincerely,
TAY wee-beng

[TAY wee-beng]

Hi Timothée,

The matrix is created with 3 dof - u,v,w. So for each i,j,k, there are 3
values. Actually I got 3 eqns from the u,v,w momentum eqns. They are not
coupled together so I can also solve them individually. But I was told
it's faster to group them together. The global indices is something like
this:

i, j, k, ijk global indices, dof
1 1 1 1 u
1 1 1 2 v
1 1 1 3 w
2 1 1 4 u

...

Hope it's clearer now.

Ok, I changed the stencil, it's now working.

Thank you

Yours sincerely,

TAY wee-beng

Post by Timothée Nicolas
Can you also tell how you created the matrix ? Just in case you
created it with the 1 dof DMDA, it would not work if you try to input
values at places where it is not allocated (which could explain the
error message)
2015-09-23 17:18 GMT+09:00 Timothée Nicolas
The first thing that strikes me is your definition of the stencils
*/MatStencil :: row(6,1),col(6,7)/**/
/*
Why is it not defined with
*/MatStencil :: row(4,1),col(4,7)/**/
/*
instead ?
Where does the 6 come from ?
Timothée
Hi,
I have successfully used MatSetValuesStencil to insert values
into a Poisson eqn matrix which has 1 DOF (pressure). Now I'm
trying to insert values in a momentum eqn matrix which has 3
DOF (u,v,w)
/*[0]PETSC ERROR: --------------------- Error Message
----------------------------*//*
*//*----------------------------------*//*
*//*[0]PETSC ERROR: Argument out of range*//*
*//*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in the matrix*//*
*//*[0]PETSC ERROR: See
http://www.mcs.anl.gov/petsc/documentation/faq.html for trou*//*
*//*ble shooting.*/
I wonder what's wrong. For the momentum eqn, for each DOF, at
at node (dof,i,j,k), I have coupling i +/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j = 2,
k = 2, which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*/
/**/PetscScalar :: value_insert(7)/**/
/**/
/**/MatStencil :: row(6,1),col(6,7)/**/
/**/
/**/ione = 1; iseven = 7/**/
/**/
/**/if (cell_type == 'u') then/**/
/**/
/**/ offset = 1/**/
/**//**/
/**/else if (cell_type == 'v') then/**/
/**/
/**/ offset = 2/**/
/**//**/
/**/else if (cell_type == 'w') then/**/
/**/
/**/ offset = 3/**/
/**//**/
/**/end if/**/
/**/
/**/do k=ksta2,kend2/**/
/**/
/**/ do j = jsta2,jend2/**/
/**/
/**/ do i=2,size_x-1/**/
/**//**/
/**/ row(MatStencil_i,1) = i - 1/**/
/**//**/
/**/ row(MatStencil_j,1) = j - 1/**/
/**//**/
/**/ row(MatStencil_k,1) = k - 1/**/
/**//**/
/**/ row(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/ value_insert = 0.d0/**/
/**//**/
/**/ col(MatStencil_i,3) = i + 1 - 1 !east/**/
/**//**/
/**/ col(MatStencil_j,3) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,3) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,3) = offset - 1/**/
/**//**/
/**/ value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,5) = i - 1 - 1 !west/**/
/**//**/
/**/ col(MatStencil_j,5) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,5) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,5) = offset - 1/**/
/**//**/
/**/ value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,2) = i - 1 !north/**/
/**//**/
/**/ col(MatStencil_j,2) = j + 1 - 1/**/
/**//**/
/**/ col(MatStencil_k,2) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,2) = offset - 1/**/
/**//**/
/**/ value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,4) = i - 1 !south/**/
/**//**/
/**/ col(MatStencil_j,4) = j - 1 - 1/**/
/**//**/
/**/ col(MatStencil_k,4) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,4) = offset - 1/**/
/**//**/
/**/ value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,6) = i - 1 !front/**/
/**//**/
/**/ col(MatStencil_j,6) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,6) = k + 1 - 1/**/
/**//**/
/**/ col(MatStencil_c,6) = offset - 1/**/
/**//**/
/**/ value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,7) = i - 1 !back/**/
/**//**/
/**/ col(MatStencil_j,7) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,7) = k - 1 - 1/**/
/**//**/
/**/ col(MatStencil_c,7) = offset - 1/**/
/**//**/
/**/ value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,1) = i - 1/**/
/**//**/
/**/ col(MatStencil_j,1) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,1) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/ value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))/**/
/**//**/
/**/ call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)/**/
/**//**/
/**/ end do/**/
/**//**/
/**/ end do/**/
/**/
/**/end do /*
Thank you
Yours sincerely,
TAY wee-beng

[Timothée Nicolas]

Yes, I had understood that, I am doing the same, but with 8 dof. This does
not change the declaration for row and col.

Can you send (i) the commands you use to create the DMDA, (ii) the commands
to create the matrix and (iii) those for the definitions of
ksta2,kend2,jsta2,jend2,size_x ?

Timothée

Post by TAY wee-beng
Hi Timothée,
The matrix is created with 3 dof - u,v,w. So for each i,j,k, there are 3
values. Actually I got 3 eqns from the u,v,w momentum eqns. They are not
coupled together so I can also solve them individually. But I was told it's
i, j, k, ijk global indices, dof
1 1 1 1 u
1 1 1 2 v
1 1 1 3 w
2 1 1 4 u
...
Hope it's clearer now.
Ok, I changed the stencil, it's now working.
Thank you
Yours sincerely,
TAY wee-beng
Can you also tell how you created the matrix ? Just in case you created it
with the 1 dof DMDA, it would not work if you try to input values at places
where it is not allocated (which could explain the error message)

Post by Timothée Nicolas
The first thing that strikes me is your definition of the stencils
*MatStencil :: row(6,1),col(6,7)*
Why is it not defined with
*MatStencil :: row(4,1),col(4,7)*
instead ?
Where does the 6 come from ?
Timothée

Post by TAY wee-beng
Hi,
I have successfully used MatSetValuesStencil to insert values into a
Poisson eqn matrix which has 1 DOF (pressure). Now I'm trying to insert
values in a momentum eqn matrix which has 3 DOF (u,v,w)
*[0]PETSC ERROR: --------------------- Error Message
----------------------------*
*----------------------------------*
*[0]PETSC ERROR: Argument out of range*
*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in the matrix*
*[0]PETSC ERROR: See
<http://www.mcs.anl.gov/petsc/documentation/faq.html>http://www.mcs.anl.gov/petsc/documentation/faq.html
<http://www.mcs.anl.gov/petsc/documentation/faq.html> for trou*
*ble shooting.*
I wonder what's wrong. For the momentum eqn, for each DOF, at at node
(dof,i,j,k), I have coupling i +/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j = 2, k = 2,
which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*PetscScalar :: value_insert(7)*
*MatStencil :: row(6,1),col(6,7)*
*ione = 1; iseven = 7*
*if (cell_type == 'u') then*
* offset = 1*
*else if (cell_type == 'v') then*
* offset = 2*
*else if (cell_type == 'w') then*
* offset = 3*
*end if*
*do k=ksta2,kend2*
* do j = jsta2,jend2*
* do i=2,size_x-1*
* row(MatStencil_i,1) = i - 1*
* row(MatStencil_j,1) = j - 1*
* row(MatStencil_k,1) = k - 1*
* row(MatStencil_c,1) = offset - 1*
* value_insert = 0.d0*
* col(MatStencil_i,3) = i + 1 - 1 !east*
* col(MatStencil_j,3) = j - 1*
* col(MatStencil_k,3) = k - 1*
* col(MatStencil_c,3) = offset - 1*
* value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,5) = i - 1 - 1 !west*
* col(MatStencil_j,5) = j - 1*
* col(MatStencil_k,5) = k - 1*
* col(MatStencil_c,5) = offset - 1*
* value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,2) = i - 1 !north*
* col(MatStencil_j,2) = j + 1 - 1*
* col(MatStencil_k,2) = k - 1*
* col(MatStencil_c,2) = offset - 1*
* value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,4) = i - 1 !south*
* col(MatStencil_j,4) = j - 1 - 1*
* col(MatStencil_k,4) = k - 1*
* col(MatStencil_c,4) = offset - 1*
* value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,6) = i - 1 !front*
* col(MatStencil_j,6) = j - 1*
* col(MatStencil_k,6) = k + 1 - 1*
* col(MatStencil_c,6) = offset - 1*
* value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,7) = i - 1 !back*
* col(MatStencil_j,7) = j - 1*
* col(MatStencil_k,7) = k - 1 - 1*
* col(MatStencil_c,7) = offset - 1*
* value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,1) = i - 1*
* col(MatStencil_j,1) = j - 1*
* col(MatStencil_k,1) = k - 1*
* col(MatStencil_c,1) = offset - 1*
* value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))*
* call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)*
* end do*
* end do*
*end do *
Thank you
Yours sincerely,
TAY wee-beng

[TAY wee-beng]

Hi Timothée,

Maybe I can send you part of it 1st. I'm trying to pinpoint why my
matrix using MatView shows zero for a lot of the values

For i=1,j=1,k=1,

It should be :

Mat Object: 1 MPI processes
type: seqaij
row 0: (0, 2) (12, -2)

but now it's:

row 0: (0, 0) (1, 0) (2, 0) (3, 0) (4, 0) (5, 0) (6, 0) (7, 0)
(8, 0) (12, 0) (13, 0) (14, 0) (24, 0) (25, 0) (26, 0) (96, 0)
(97, 0) (98, 0) (192, 0) (193, 0) (194, 0)

I used:

/call
DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&//
//
//size_z,1,PETSC_DECIDE,PETSC_DECIDE,3,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_uvw_ast,ierr)//
//
//call DMSetMatType(da_uvw_ast,MATAIJ,ierr)//
////
// call DMCreateMatrix(da_uvw_ast,A_semi_xyz,ierr)//
//
// call MatSetFromOptions(A_semi_xyz,ierr)//
//
//i = 1, j = 1, k = 1//
//
//row(MatStencil_i,1) = i - 1//
////
// row(MatStencil_j,1) = j - 1//
////
// row(MatStencil_k,1) = k - 1//
////
// row(MatStencil_c,1) = 0//
////
// 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) = 2.//
//
// call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,1),value_insert(1),INSERT_VALUES,ierr)//
//
// 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) = -2.//
//
// call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,2),value_insert(2),INSERT_VALUES,ierr)/

Thank you

Yours sincerely,

TAY wee-beng

Post by Timothée Nicolas
Yes, I had understood that, I am doing the same, but with 8 dof. This
does not change the declaration for row and col.
Can you send (i) the commands you use to create the DMDA, (ii) the
commands to create the matrix and (iii) those for the definitions of
ksta2,kend2,jsta2,jend2,size_x ?
Timothée
Hi Timothée,
The matrix is created with 3 dof - u,v,w. So for each i,j,k, there
are 3 values. Actually I got 3 eqns from the u,v,w momentum eqns.
They are not coupled together so I can also solve them
individually. But I was told it's faster to group them together.
i, j, k, ijk global indices, dof
1 1 1 1 u
1 1 1 2 v
1 1 1 3 w
2 1 1 4 u
...
Hope it's clearer now.
Ok, I changed the stencil, it's now working.
Thank you
Yours sincerely,
TAY wee-beng

Post by Timothée Nicolas
Can you also tell how you created the matrix ? Just in case you
created it with the 1 dof DMDA, it would not work if you try to
input values at places where it is not allocated (which could
explain the error message)
2015-09-23 17:18 GMT+09:00 Timothée Nicolas
The first thing that strikes me is your definition of the stencils
*/MatStencil :: row(6,1),col(6,7)/**/
/*
Why is it not defined with
*/MatStencil :: row(4,1),col(4,7)/**/
/*
instead ?
Where does the 6 come from ?
Timothée
Hi,
I have successfully used MatSetValuesStencil to insert
values into a Poisson eqn matrix which has 1 DOF
(pressure). Now I'm trying to insert values in a momentum
eqn matrix which has 3 DOF (u,v,w)
/*[0]PETSC ERROR: --------------------- Error Message
----------------------------*//*
*//*----------------------------------*//*
*//*[0]PETSC ERROR: Argument out of range*//*
*//*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in
the matrix*//*
*//*[0]PETSC ERROR: See
http://www.mcs.anl.gov/petsc/documentation/faq.html for trou*//*
*//*ble shooting.*/
I wonder what's wrong. For the momentum eqn, for each
DOF, at at node (dof,i,j,k), I have coupling i +/- 1, j
+/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j
= 2, k = 2, which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*/
/**/PetscScalar :: value_insert(7)/**/
/**/
/**/MatStencil :: row(6,1),col(6,7)/**/
/**/
/**/ione = 1; iseven = 7/**/
/**/
/**/if (cell_type == 'u') then/**/
/**/
/**/ offset = 1/**/
/**//**/
/**/else if (cell_type == 'v') then/**/
/**/
/**/ offset = 2/**/
/**//**/
/**/else if (cell_type == 'w') then/**/
/**/
/**/ offset = 3/**/
/**//**/
/**/end if/**/
/**/
/**/do k=ksta2,kend2/**/
/**/
/**/ do j = jsta2,jend2/**/
/**/
/**/ do i=2,size_x-1/**/
/**//**/
/**/row(MatStencil_i,1) = i - 1/**/
/**//**/
/**/row(MatStencil_j,1) = j - 1/**/
/**//**/
/**/row(MatStencil_k,1) = k - 1/**/
/**//**/
/**/row(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/value_insert = 0.d0/**/
/**//**/
/**/col(MatStencil_i,3) = i + 1 - 1 !east/**/
/**//**/
/**/col(MatStencil_j,3) = j - 1/**/
/**//**/
/**/col(MatStencil_k,3) = k - 1/**/
/**//**/
/**/col(MatStencil_c,3) = offset - 1/**/
/**//**/
/**/value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,5) = i - 1 - 1 !west/**/
/**//**/
/**/col(MatStencil_j,5) = j - 1/**/
/**//**/
/**/col(MatStencil_k,5) = k - 1/**/
/**//**/
/**/col(MatStencil_c,5) = offset - 1/**/
/**//**/
/**/value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,2) = i - 1 !north/**/
/**//**/
/**/col(MatStencil_j,2) = j + 1 - 1/**/
/**//**/
/**/col(MatStencil_k,2) = k - 1/**/
/**//**/
/**/col(MatStencil_c,2) = offset - 1/**/
/**//**/
/**/value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,4) = i - 1 !south/**/
/**//**/
/**/col(MatStencil_j,4) = j - 1 - 1/**/
/**//**/
/**/col(MatStencil_k,4) = k - 1/**/
/**//**/
/**/col(MatStencil_c,4) = offset - 1/**/
/**//**/
/**/value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,6) = i - 1 !front/**/
/**//**/
/**/col(MatStencil_j,6) = j - 1/**/
/**//**/
/**/col(MatStencil_k,6) = k + 1 - 1/**/
/**//**/
/**/col(MatStencil_c,6) = offset - 1/**/
/**//**/
/**/value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,7) = i - 1 !back/**/
/**//**/
/**/col(MatStencil_j,7) = j - 1/**/
/**//**/
/**/col(MatStencil_k,7) = k - 1 - 1/**/
/**//**/
/**/col(MatStencil_c,7) = offset - 1/**/
/**//**/
/**/value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,1) = i - 1/**/
/**//**/
/**/col(MatStencil_j,1) = j - 1/**/
/**//**/
/**/col(MatStencil_k,1) = k - 1/**/
/**//**/
/**/col(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))/**/
/**//**/
/**/ call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)/**/
/**//**/
/**/ end do/**/
/**//**/
/**/ end do/**/
/**/
/**/end do /*
Thank you
Yours sincerely,
TAY wee-beng

[Timothée Nicolas]

Hum, sorry, I don't know. I asked you to provide the definitions of start
and end values of i,j,k, because I was concerned whether you may mess up
the boundaries. Especially because you seem to treat x differently from y
and z. You have the problem also on only 1 process ?

Post by TAY wee-beng
Hi Timothée,
Maybe I can send you part of it 1st. I'm trying to pinpoint why my matrix
using MatView shows zero for a lot of the values
For i=1,j=1,k=1,
Mat Object: 1 MPI processes
type: seqaij
row 0: (0, 2) (12, -2)
row 0: (0, 0) (1, 0) (2, 0) (3, 0) (4, 0) (5, 0) (6, 0) (7, 0) (8,
0) (12, 0) (13, 0) (14, 0) (24, 0) (25, 0) (26, 0) (96, 0) (97, 0)
(98, 0) (192, 0) (193, 0) (194, 0)
*call
DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&*
*size_z,1,PETSC_DECIDE,PETSC_DECIDE,3,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_uvw_ast,ierr)*
*call DMSetMatType(da_uvw_ast,MATAIJ,ierr)*
* call DMCreateMatrix(da_uvw_ast,A_semi_xyz,ierr)*
* call MatSetFromOptions(A_semi_xyz,ierr)*
*i = 1, j = 1, k = 1*
*row(MatStencil_i,1) = i - 1*
* row(MatStencil_j,1) = j - 1*
* row(MatStencil_k,1) = k - 1*
* row(MatStencil_c,1) = 0*
* 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) = 2.*
* call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,1),value_insert(1),INSERT_VALUES,ierr)*
* 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) = -2.*
* call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,2),value_insert(2),INSERT_VALUES,ierr)*
Thank you
Yours sincerely,
TAY wee-beng
Yes, I had understood that, I am doing the same, but with 8 dof. This does
not change the declaration for row and col.
Can you send (i) the commands you use to create the DMDA, (ii) the
commands to create the matrix and (iii) those for the definitions of
ksta2,kend2,jsta2,jend2,size_x ?
Timothée

Post by TAY wee-beng
Hi Timothée,
The matrix is created with 3 dof - u,v,w. So for each i,j,k, there are 3
values. Actually I got 3 eqns from the u,v,w momentum eqns. They are not
coupled together so I can also solve them individually. But I was told it's
i, j, k, ijk global indices, dof
1 1 1 1 u
1 1 1 2 v
1 1 1 3 w
2 1 1 4 u
...
Hope it's clearer now.
Ok, I changed the stencil, it's now working.
Thank you
Yours sincerely,
TAY wee-beng
Can you also tell how you created the matrix ? Just in case you created
it with the 1 dof DMDA, it would not work if you try to input values at
places where it is not allocated (which could explain the error message)
2015-09-23 17:18 GMT+09:00 Timothée Nicolas <

Post by Timothée Nicolas
The first thing that strikes me is your definition of the stencils
*MatStencil :: row(6,1),col(6,7)*
Why is it not defined with
*MatStencil :: row(4,1),col(4,7)*
instead ?
Where does the 6 come from ?
Timothée

Post by TAY wee-beng
Hi,
I have successfully used MatSetValuesStencil to insert values into a
Poisson eqn matrix which has 1 DOF (pressure). Now I'm trying to insert
values in a momentum eqn matrix which has 3 DOF (u,v,w)
*[0]PETSC ERROR: --------------------- Error Message
----------------------------*
*----------------------------------*
*[0]PETSC ERROR: Argument out of range*
*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in the matrix*
*[0]PETSC ERROR: See
http://www.mcs.anl.gov/petsc/documentation/faq.html
<http://www.mcs.anl.gov/petsc/documentation/faq.html> for trou*
*ble shooting.*
I wonder what's wrong. For the momentum eqn, for each DOF, at at node
(dof,i,j,k), I have coupling i +/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j = 2, k = 2,
which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*PetscScalar :: value_insert(7)*
*MatStencil :: row(6,1),col(6,7)*
*ione = 1; iseven = 7*
*if (cell_type == 'u') then*
* offset = 1*
*else if (cell_type == 'v') then*
* offset = 2*
*else if (cell_type == 'w') then*
* offset = 3*
*end if*
*do k=ksta2,kend2*
* do j = jsta2,jend2*
* do i=2,size_x-1*
* row(MatStencil_i,1) = i - 1*
* row(MatStencil_j,1) = j - 1*
* row(MatStencil_k,1) = k - 1*
* row(MatStencil_c,1) = offset - 1*
* value_insert = 0.d0*
* col(MatStencil_i,3) = i + 1 - 1 !east*
* col(MatStencil_j,3) = j - 1*
* col(MatStencil_k,3) = k - 1*
* col(MatStencil_c,3) = offset - 1*
* value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,5) = i - 1 - 1 !west*
* col(MatStencil_j,5) = j - 1*
* col(MatStencil_k,5) = k - 1*
* col(MatStencil_c,5) = offset - 1*
* value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,2) = i - 1 !north*
* col(MatStencil_j,2) = j + 1 - 1*
* col(MatStencil_k,2) = k - 1*
* col(MatStencil_c,2) = offset - 1*
* value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,4) = i - 1 !south*
* col(MatStencil_j,4) = j - 1 - 1*
* col(MatStencil_k,4) = k - 1*
* col(MatStencil_c,4) = offset - 1*
* value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,6) = i - 1 !front*
* col(MatStencil_j,6) = j - 1*
* col(MatStencil_k,6) = k + 1 - 1*
* col(MatStencil_c,6) = offset - 1*
* value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,7) = i - 1 !back*
* col(MatStencil_j,7) = j - 1*
* col(MatStencil_k,7) = k - 1 - 1*
* col(MatStencil_c,7) = offset - 1*
* value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,1) = i - 1*
* col(MatStencil_j,1) = j - 1*
* col(MatStencil_k,1) = k - 1*
* col(MatStencil_c,1) = offset - 1*
* value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))*
* call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)*
* end do*
* end do*
*end do *
Thank you
Yours sincerely,
TAY wee-beng

[TAY wee-beng]

Hi Timothée,

I got the error sorted out. I was inserting into the wrong matrix. But
the key error was the defn of row and col.

Thank you for pointing it out.

Yours sincerely,

TAY wee-beng

Post by Timothée Nicolas
Hum, sorry, I don't know. I asked you to provide the definitions of
start and end values of i,j,k, because I was concerned whether you may
mess up the boundaries. Especially because you seem to treat x
differently from y and z. You have the problem also on only 1 process ?
Hi Timothée,
Maybe I can send you part of it 1st. I'm trying to pinpoint why my
matrix using MatView shows zero for a lot of the values
For i=1,j=1,k=1,
Mat Object: 1 MPI processes
type: seqaij
row 0: (0, 2) (12, -2)
row 0: (0, 0) (1, 0) (2, 0) (3, 0) (4, 0) (5, 0) (6, 0) (7,
0) (8, 0) (12, 0) (13, 0) (14, 0) (24, 0) (25, 0) (26, 0)
(96, 0) (97, 0) (98, 0) (192, 0) (193, 0) (194, 0)
/call
DMDACreate3d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,size_x,size_y,&//
//
//size_z,1,PETSC_DECIDE,PETSC_DECIDE,3,stencil_width,lx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da_uvw_ast,ierr)//
//
//call DMSetMatType(da_uvw_ast,MATAIJ,ierr)//
////
// call DMCreateMatrix(da_uvw_ast,A_semi_xyz,ierr)//
//
// call MatSetFromOptions(A_semi_xyz,ierr)//
//
//i = 1, j = 1, k = 1//
//
//row(MatStencil_i,1) = i - 1//
////
// row(MatStencil_j,1) = j - 1//
////
// row(MatStencil_k,1) = k - 1//
////
// row(MatStencil_c,1) = 0//
////
// 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) = 2.//
//
// call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,1),value_insert(1),INSERT_VALUES,ierr)//
//
// 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) = -2.//
//
// call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,2),value_insert(2),INSERT_VALUES,ierr)/
Thank you
Yours sincerely,
TAY wee-beng

Post by Timothée Nicolas
Yes, I had understood that, I am doing the same, but with 8 dof.
This does not change the declaration for row and col.
Can you send (i) the commands you use to create the DMDA, (ii)
the commands to create the matrix and (iii) those for the
definitions of ksta2,kend2,jsta2,jend2,size_x ?
Timothée
Hi Timothée,
The matrix is created with 3 dof - u,v,w. So for each i,j,k,
there are 3 values. Actually I got 3 eqns from the u,v,w
momentum eqns. They are not coupled together so I can also
solve them individually. But I was told it's faster to group
i, j, k, ijk global indices, dof
1 1 1 1 u
1 1 1 2 v
1 1 1 3 w
2 1 1 4 u
...
Hope it's clearer now.
Ok, I changed the stencil, it's now working.
Thank you
Yours sincerely,
TAY wee-beng

Post by Timothée Nicolas
Can you also tell how you created the matrix ? Just in case
you created it with the 1 dof DMDA, it would not work if you
try to input values at places where it is not allocated
(which could explain the error message)
2015-09-23 17:18 GMT+09:00 Timothée Nicolas
The first thing that strikes me is your definition of
the stencils
*/MatStencil :: row(6,1),col(6,7)/**/
/*
Why is it not defined with
*/MatStencil :: row(4,1),col(4,7)/**/
/*
instead ?
Where does the 6 come from ?
Timothée
2015-09-23 17:14 GMT+09:00 TAY wee-beng
Hi,
I have successfully used MatSetValuesStencil to
insert values into a Poisson eqn matrix which has 1
DOF (pressure). Now I'm trying to insert values in a
momentum eqn matrix which has 3 DOF (u,v,w)
/*[0]PETSC ERROR: --------------------- Error
Message ----------------------------*//*
*//*----------------------------------*//*
*//*[0]PETSC ERROR: Argument out of range*//*
*//*[0]PETSC ERROR: Inserting a new nonzero at
(111,5) in the matrix*//*
*//*[0]PETSC ERROR: See
http://www.mcs.anl.gov/petsc/documentation/faq.html
for trou*//*
*//*ble shooting.*/
I wonder what's wrong. For the momentum eqn, for
each DOF, at at node (dof,i,j,k), I have coupling i
+/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i =
2, j = 2, k = 2, which is an internal node.
Here's part of my code below. Hope someone can help.
Thanks!
*/
/**/PetscScalar :: value_insert(7)/**/
/**/
/**/MatStencil :: row(6,1),col(6,7)/**/
/**/
/**/ione = 1; iseven = 7/**/
/**/
/**/if (cell_type == 'u') then/**/
/**/
/**/offset = 1/**/
/**//**/
/**/else if (cell_type == 'v') then/**/
/**/
/**/offset = 2/**/
/**//**/
/**/else if (cell_type == 'w') then/**/
/**/
/**/offset = 3/**/
/**//**/
/**/end if/**/
/**/
/**/do k=ksta2,kend2/**/
/**/
/**/do j = jsta2,jend2/**/
/**/
/**/ do i=2,size_x-1/**/
/**//**/
/**/row(MatStencil_i,1) = i - 1/**/
/**//**/
/**/row(MatStencil_j,1) = j - 1/**/
/**//**/
/**/row(MatStencil_k,1) = k - 1/**/
/**//**/
/**/row(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/value_insert = 0.d0/**/
/**//**/
/**/col(MatStencil_i,3) = i + 1 - 1 !east/**/
/**//**/
/**/col(MatStencil_j,3) = j - 1/**/
/**//**/
/**/col(MatStencil_k,3) = k - 1/**/
/**//**/
/**/col(MatStencil_c,3) = offset - 1/**/
/**//**/
/**/value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,5) = i - 1 - 1 !west/**/
/**//**/
/**/col(MatStencil_j,5) = j - 1/**/
/**//**/
/**/col(MatStencil_k,5) = k - 1/**/
/**//**/
/**/col(MatStencil_c,5) = offset - 1/**/
/**//**/
/**/value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,2) = i - 1 !north/**/
/**//**/
/**/col(MatStencil_j,2) = j + 1 - 1/**/
/**//**/
/**/col(MatStencil_k,2) = k - 1/**/
/**//**/
/**/col(MatStencil_c,2) = offset - 1/**/
/**//**/
/**/value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,4) = i - 1 !south/**/
/**//**/
/**/col(MatStencil_j,4) = j - 1 - 1/**/
/**//**/
/**/col(MatStencil_k,4) = k - 1/**/
/**//**/
/**/col(MatStencil_c,4) = offset - 1/**/
/**//**/
/**/value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,6) = i - 1 !front/**/
/**//**/
/**/col(MatStencil_j,6) = j - 1/**/
/**//**/
/**/col(MatStencil_k,6) = k + 1 - 1/**/
/**//**/
/**/col(MatStencil_c,6) = offset - 1/**/
/**//**/
/**/value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,7) = i - 1 !back/**/
/**//**/
/**/col(MatStencil_j,7) = j - 1/**/
/**//**/
/**/col(MatStencil_k,7) = k - 1 - 1/**/
/**//**/
/**/col(MatStencil_c,7) = offset - 1/**/
/**//**/
/**/value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/col(MatStencil_i,1) = i - 1/**/
/**//**/
/**/col(MatStencil_j,1) = j - 1/**/
/**//**/
/**/col(MatStencil_k,1) = k - 1/**/
/**//**/
/**/col(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))/**/
/**//**/
/**/call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)/**/
/**//**/
/**/end do/**/
/**//**/
/**/end do/**/
/**/
/**/end do /*
Thank you
Yours sincerely,
TAY wee-beng

[TAY wee-beng]

Thank you

Yours sincerely,

TAY wee-beng

Post by Timothée Nicolas
The first thing that strikes me is your definition of the stencils
*/MatStencil :: row(6,1),col(6,7)/**/
/*
Why is it not defined with
*/MatStencil :: row(4,1),col(4,7)/**/
/*
instead ?
Where does the 6 come from ?
Timothée

Hi,

Maybe I misunderstood. But in your email tutorial, last time, with 1
dof, it's


*/MatStencil :: row(4,1),col(4,7)/**/
/*
So now I have 3 dof, shouldn't it be

*/MatStencil :: row(6,1),col(6,7)/**/
?

/*

Post by Timothée Nicolas
Hi,
I have successfully used MatSetValuesStencil to insert values into
a Poisson eqn matrix which has 1 DOF (pressure). Now I'm trying to
insert values in a momentum eqn matrix which has 3 DOF (u,v,w)
/*[0]PETSC ERROR: --------------------- Error Message
----------------------------*//*
*//*----------------------------------*//*
*//*[0]PETSC ERROR: Argument out of range*//*
*//*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in the matrix*//*
*//*[0]PETSC ERROR: See
http://www.mcs.anl.gov/petsc/documentation/faq.html for trou*//*
*//*ble shooting.*/
I wonder what's wrong. For the momentum eqn, for each DOF, at at
node (dof,i,j,k), I have coupling i +/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j = 2, k =
2, which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*/
/**/PetscScalar :: value_insert(7)/**/
/**/
/**/MatStencil :: row(6,1),col(6,7)/**/
/**/
/**/ione = 1; iseven = 7/**/
/**/
/**/if (cell_type == 'u') then/**/
/**/
/**/ offset = 1/**/
/**//**/
/**/else if (cell_type == 'v') then/**/
/**/
/**/ offset = 2/**/
/**//**/
/**/else if (cell_type == 'w') then/**/
/**/
/**/ offset = 3/**/
/**//**/
/**/end if/**/
/**/
/**/do k=ksta2,kend2/**/
/**/
/**/ do j = jsta2,jend2/**/
/**/
/**/ do i=2,size_x-1/**/
/**//**/
/**/ row(MatStencil_i,1) = i - 1/**/
/**//**/
/**/ row(MatStencil_j,1) = j - 1/**/
/**//**/
/**/ row(MatStencil_k,1) = k - 1/**/
/**//**/
/**/ row(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/ value_insert = 0.d0/**/
/**//**/
/**/ col(MatStencil_i,3) = i + 1 - 1 !east/**/
/**//**/
/**/ col(MatStencil_j,3) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,3) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,3) = offset - 1/**/
/**//**/
/**/ value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,5) = i - 1 - 1 !west/**/
/**//**/
/**/ col(MatStencil_j,5) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,5) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,5) = offset - 1/**/
/**//**/
/**/ value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,2) = i - 1 !north/**/
/**//**/
/**/ col(MatStencil_j,2) = j + 1 - 1/**/
/**//**/
/**/ col(MatStencil_k,2) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,2) = offset - 1/**/
/**//**/
/**/ value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,4) = i - 1 !south/**/
/**//**/
/**/ col(MatStencil_j,4) = j - 1 - 1/**/
/**//**/
/**/ col(MatStencil_k,4) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,4) = offset - 1/**/
/**//**/
/**/ value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,6) = i - 1 !front/**/
/**//**/
/**/ col(MatStencil_j,6) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,6) = k + 1 - 1/**/
/**//**/
/**/ col(MatStencil_c,6) = offset - 1/**/
/**//**/
/**/ value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,7) = i - 1 !back/**/
/**//**/
/**/ col(MatStencil_j,7) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,7) = k - 1 - 1/**/
/**//**/
/**/ col(MatStencil_c,7) = offset - 1/**/
/**//**/
/**/ value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re/**/
/**//**/
/**/ col(MatStencil_i,1) = i - 1/**/
/**//**/
/**/ col(MatStencil_j,1) = j - 1/**/
/**//**/
/**/ col(MatStencil_k,1) = k - 1/**/
/**//**/
/**/ col(MatStencil_c,1) = offset - 1/**/
/**//**/
/**/ value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))/**/
/**//**/
/**/ call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)/**/
/**//**/
/**/ end do/**/
/**//**/
/**/ end do/**/
/**/
/**/end do /*
Thank you
Yours sincerely,
TAY wee-beng

[Timothée Nicolas]

No, you misunderstood. The 4 refers to

1 for i direction
+
1 for j direction
+
1 for k direction
+
1 for dof direction
=
4

Does it solve your problem ?

Best

Timothee

Post by TAY wee-beng
Thank you
Yours sincerely,
TAY wee-beng
The first thing that strikes me is your definition of the stencils
*MatStencil :: row(6,1),col(6,7)*
Why is it not defined with
*MatStencil :: row(4,1),col(4,7)*
instead ?
Where does the 6 come from ?
Timothée
Hi,
Maybe I misunderstood. But in your email tutorial, last time, with 1 dof,
it's
*MatStencil :: row(4,1),col(4,7)*
So now I have 3 dof, shouldn't it be
*MatStencil :: row(6,1),col(6,7)*
* ? *

Post by TAY wee-beng
Hi,
I have successfully used MatSetValuesStencil to insert values into a
Poisson eqn matrix which has 1 DOF (pressure). Now I'm trying to insert
values in a momentum eqn matrix which has 3 DOF (u,v,w)
*[0]PETSC ERROR: --------------------- Error Message
----------------------------*
*----------------------------------*
*[0]PETSC ERROR: Argument out of range*
*[0]PETSC ERROR: Inserting a new nonzero at (111,5) in the matrix*
*[0]PETSC ERROR: See
<http://www.mcs.anl.gov/petsc/documentation/faq.html>http://www.mcs.anl.gov/petsc/documentation/faq.html
<http://www.mcs.anl.gov/petsc/documentation/faq.html> for trou*
*ble shooting.*
I wonder what's wrong. For the momentum eqn, for each DOF, at at node
(dof,i,j,k), I have coupling i +/- 1, j +/- 1 and k +/- 1.
The error happens at 111,5, which corresponds to i = 2, j = 2, k = 2,
which is an internal node.
Here's part of my code below. Hope someone can help. Thanks!
*PetscScalar :: value_insert(7)*
*MatStencil :: row(6,1),col(6,7)*
*ione = 1; iseven = 7*
*if (cell_type == 'u') then*
* offset = 1*
*else if (cell_type == 'v') then*
* offset = 2*
*else if (cell_type == 'w') then*
* offset = 3*
*end if*
*do k=ksta2,kend2*
* do j = jsta2,jend2*
* do i=2,size_x-1*
* row(MatStencil_i,1) = i - 1*
* row(MatStencil_j,1) = j - 1*
* row(MatStencil_k,1) = k - 1*
* row(MatStencil_c,1) = offset - 1*
* value_insert = 0.d0*
* col(MatStencil_i,3) = i + 1 - 1 !east*
* col(MatStencil_j,3) = j - 1*
* col(MatStencil_k,3) = k - 1*
* col(MatStencil_c,3) = offset - 1*
* value_insert(3) = -(
1./(cell_x(i)%pd_E+cell_x(i+1)%pd_W))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,5) = i - 1 - 1 !west*
* col(MatStencil_j,5) = j - 1*
* col(MatStencil_k,5) = k - 1*
* col(MatStencil_c,5) = offset - 1*
* value_insert(5) = -(
1./(cell_x(i)%pd_W+cell_x(i-1)%pd_E))*(c_yz(j,k)%fc_E)*inv_Re*
* col(MatStencil_i,2) = i - 1 !north*
* col(MatStencil_j,2) = j + 1 - 1*
* col(MatStencil_k,2) = k - 1*
* col(MatStencil_c,2) = offset - 1*
* value_insert(2) = -(
1./(cell_y(j)%pd_N+cell_y(j+1)%pd_S))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,4) = i - 1 !south*
* col(MatStencil_j,4) = j - 1 - 1*
* col(MatStencil_k,4) = k - 1*
* col(MatStencil_c,4) = offset - 1*
* value_insert(4) = -(
1./(cell_y(j)%pd_S+cell_y(j-1)%pd_N))*(c_zx(i,k)%fc_N)*inv_Re*
* col(MatStencil_i,6) = i - 1 !front*
* col(MatStencil_j,6) = j - 1*
* col(MatStencil_k,6) = k + 1 - 1*
* col(MatStencil_c,6) = offset - 1*
* value_insert(6) = -(
1./(cell_z(k)%pd_F+cell_z(k+1)%pd_B))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,7) = i - 1 !back*
* col(MatStencil_j,7) = j - 1*
* col(MatStencil_k,7) = k - 1 - 1*
* col(MatStencil_c,7) = offset - 1*
* value_insert(7) = -(
1./(cell_z(k)%pd_B+cell_z(k-1)%pd_F))*(c_xy(i,j)%fc_F)*inv_Re*
* col(MatStencil_i,1) = i - 1*
* col(MatStencil_j,1) = j - 1*
* col(MatStencil_k,1) = k - 1*
* col(MatStencil_c,1) = offset - 1*
* value_insert(1) = 2.*c(i,j,k)%vol/del_t -
(value_insert(2)+value_insert(3)+value_insert(4)+value_insert(5)+value_insert(6)+value_insert(7))*
* call
MatSetValuesStencil(A_semi_xyz,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)*
* end do*
* end do*
*end do *
Thank you
Yours sincerely,
TAY wee-beng

[Matthew Knepley]

Post by Wee-Beng Tay
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?

No. MatSetValuesLocal() takes local indices. MatSetValuesStencil() takes
global vertex numbers.

Post by Wee-Beng Tay
Which is a simpler/better option?

MatSetValuesStencil()

Post by Wee-Beng Tay
Is there an example in Fortran for MatSetValuesStencil?

Timothée Nicolas shows one in his reply.

Do I also need to use DMDAGetAO together with MatSetValuesStencil or

Post by Wee-Beng Tay
MatSetValuesLocal?

No.

Thanks,

Matt

Post by Wee-Beng Tay
Thanks!

[Wee Beng Tay]

Sent using CloudMagic [https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2] On Mon, Aug 24, 2015 at 6:21 PM, Matthew Knepley < ***@gmail.com [***@gmail.com] > wrote:
On Mon, Aug 24, 2015 at 4:09 AM, Wee-Beng Tay < ***@gmail.com [***@gmail.com] > wrote:
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?
No. MatSetValuesLocal() takes local indices. MatSetValuesStencil() takes global
vertex numbers.
So MatSetValuesStencil() takes global vertex numbers. Do you mean the natural or
petsc ordering?
Which is a simpler/better option?
MatSetValuesStencil() Is there an example in Fortran for MatSetValuesStencil?
Timothée Nicolas shows one in his reply.

Do I also need to use DMDAGetAO together with MatSetValuesStencil or
MatSetValuesLocal?
No.
Thanks,
Matt Thanks!


--
What most experimenters take for granted before they begin their experiments is
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

[Matthew Knepley]

Post by Wee Beng Tay
Sent using CloudMagic
<https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2>

Post by Wee-Beng Tay
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?

No. MatSetValuesLocal() takes local indices. MatSetValuesStencil() takes
global vertex numbers.
So MatSetValuesStencil() takes global vertex numbers. Do you mean the
natural or petsc ordering?

There is no PETSc ordering for vertices, only the natural ordering.

Thanks,

Matt

Post by Wee Beng Tay
Which is a simpler/better option?
MatSetValuesStencil()

Post by Wee-Beng Tay
Is there an example in Fortran for MatSetValuesStencil?

Timothée Nicolas shows one in his reply.
Do I also need to use DMDAGetAO together with MatSetValuesStencil or

Post by Wee-Beng Tay
MatSetValuesLocal?

No.
Thanks,
Matt

Post by Wee-Beng Tay
Thanks!

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

[Wee Beng Tay]

Hi,

So can I use multigrid directly if using matsetvalues stencil?

Thanks



Sent using CloudMagic [https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2] On Tue, Aug 25, 2015 at 10:11 AM, Matthew Knepley < ***@gmail.com [***@gmail.com] > wrote:
On Mon, Aug 24, 2015 at 9:01 PM, Wee Beng Tay < ***@gmail.com [***@gmail.com] > wrote:


Sent using CloudMagic [https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2] On Mon, Aug 24, 2015 at 6:21 PM, Matthew Knepley < ***@gmail.com [***@gmail.com] > wrote:
On Mon, Aug 24, 2015 at 4:09 AM, Wee-Beng Tay < ***@gmail.com [***@gmail.com] > wrote:
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?
No. MatSetValuesLocal() takes local indices. MatSetValuesStencil() takes global
vertex numbers.
So MatSetValuesStencil() takes global vertex numbers. Do you mean the natural or
petsc ordering?
There is no PETSc ordering for vertices, only the natural ordering.
Thanks,
Matt Which is a simpler/better option?
MatSetValuesStencil() Is there an example in Fortran for MatSetValuesStencil?
Timothée Nicolas shows one in his reply.

Do I also need to use DMDAGetAO together with MatSetValuesStencil or
MatSetValuesLocal?
No.
Thanks,
Matt Thanks!


--
What most experimenters take for granted before they begin their experiments is
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener


--
What most experimenters take for granted before they begin their experiments is
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

[Matthew Knepley]

Post by Wee Beng Tay
Hi,
So can I use multigrid directly if using matsetvalues stencil?

Do you mean, if you use MatSetStencil() then you statements will work no
matter what grid comes in to your residual function?
That is true.

Thanks,

Matt

Post by Wee Beng Tay
Thanks
Sent using CloudMagic
<https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2>

Post by Wee Beng Tay
Sent using CloudMagic
<https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2>

Post by Wee-Beng Tay
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?

No. MatSetValuesLocal() takes local indices. MatSetValuesStencil() takes
global vertex numbers.
So MatSetValuesStencil() takes global vertex numbers. Do you mean the
natural or petsc ordering?

There is no PETSc ordering for vertices, only the natural ordering.
Thanks,
Matt

Post by Wee Beng Tay
Which is a simpler/better option?
MatSetValuesStencil()

Post by Wee-Beng Tay
Is there an example in Fortran for MatSetValuesStencil?

Timothée Nicolas shows one in his reply.
Do I also need to use DMDAGetAO together with MatSetValuesStencil or

Post by Wee-Beng Tay
MatSetValuesLocal?

No.
Thanks,
Matt

Post by Wee-Beng Tay
Thanks!

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener