<html>
<head>
<meta content="text/html; charset=utf-8" http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
Hi Timothee,<br>
<br>
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.<br>
<br>
However, now I am using only 1 procs.<br>
<br>
I was going to resend you my code but then I realised my mistake. I
used:<br>
<br>
<b><i>call
MatSetValuesStencil(A_mat,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)</i></b><b><i><br>
</i></b><br>
for all pts, including those at the boundary. Hence, those coupling
outside the boundary is also included.<br>
<br>
I changed to:<br>
<br>
<b><i>call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,7),value_insert(7),INSERT_VALUES,ierr)</i></b><br>
<br>
so I am now entering values individually.<br>
<br>
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?<br>
<pre class="moz-signature" cols="72">Thank you
Yours sincerely,
TAY wee-beng</pre>
<div class="moz-cite-prefix">On 26/8/2015 12:24 PM, Timothée Nicolas
wrote:<br>
</div>
<blockquote
cite="mid:CAGi1ndR78hB8M9_ybho-favYfdJnsissE_A45V1hFYZTQ0AC0w@mail.gmail.com"
type="cite">
<div dir="ltr">
<div>
<div>What is the definition of ksta, kend, jsta, jend ? Etc ?
You are parallelized only in j and k ?<br>
<br>
</div>
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)<br>
<br>
</div>
<div>Would you mind sending the complete routine with the
complete definitions of ksta,kend,jsta,jend,and size_x ?<br>
</div>
<div><br>
</div>
Timothee<br>
</div>
<div class="gmail_extra"><br>
<div class="gmail_quote">2015-08-26 13:12 GMT+09:00 TAY wee-beng
<span dir="ltr"><<a moz-do-not-send="true"
href="mailto:zonexo@gmail.com" target="_blank">zonexo@gmail.com</a>></span>:<br>
<blockquote class="gmail_quote" style="margin:0 0 0
.8ex;border-left:1px #ccc solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000"> Hi,<br>
<br>
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.<br>
<br>
size x/y/z = 4/8/10<br>
<br>
<b><i>MatStencil :: row(4,1),col(4,7)</i></b><b><i><br>
</i></b><b><i><br>
</i></b><b><i>PetscScalar :: value_insert(7)</i></b><b><i><br>
</i></b><b><i><br>
</i></b><b><i>PetscInt :: ione,iseven</i></b><b><i><br>
</i></b><b><i><br>
</i></b><b><i>ione = 1; iseven = 7</i></b><b><i><br>
</i></b><b><i><br>
</i></b><b><i>do k=ksta,kend</i></b><b><i><br>
</i></b><b><i><br>
</i></b><b><i> do j = jsta,jend</i></b><b><i><br>
</i></b><b><i><br>
</i></b><b><i> do i=1,size_x</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> row(MatStencil_i,1) = i -
1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> row(MatStencil_j,1) = j -
1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> row(MatStencil_k,1) = k -
1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> row(MatStencil_c,1) = 0 !
1 - 1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> value_insert = 0.d0</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> if (i /= size_x) then</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_i,3) =
i + 1 - 1 !east</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_j,3) =
j - 1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_k,3) =
k - 1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_c,3) =
0</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> value_insert(3) =
(cp_yz(j,k)%fc_E)/(cp_x(i)%pd_E+cp_x(i+1)%pd_W)</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> end if</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> if (i /= 1) then</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_i,5) =
i - 1 - 1 !west</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_j,5) =
j - 1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_k,5) =
k - 1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_c,5) =
0</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> value_insert(5) =
(cp_yz(j,k)%fc_E)/(cp_x(i)%pd_W+cp_x(i-1)%pd_E)</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> end if</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> if (j /= size_y) then</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_i,2) =
i - 1 !north</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_j,2) =
j + 1 - 1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_k,2) =
k - 1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_c,2) =
0</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> value_insert(2) =
(cp_zx(i,k)%fc_N)/(cp_y(j)%pd_N+cp_y(j+1)%pd_S)</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> end if</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> ...</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_i,1) = i -
1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_j,1) = j -
1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_k,1) = k -
1</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> col(MatStencil_c,1) = 0</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> value_insert(1) =
-value_insert(2) - value_insert(3) - value_insert(4) -
value_insert(5) - value_insert(6) - value_insert(7)</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> call
MatSetValuesStencil(A_mat,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> end do</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> end do</i></b><b><i><br>
</i></b><b><i> </i></b><b><i><br>
</i></b><b><i> end do</i></b><br>
<br>
but I got the error :<br>
<br>
[0]PETSC ERROR: Argument out of range<br>
[0]PETSC ERROR: Inserting a new nonzero at (3,0) in the
matrix.<br>
<br>
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?<br>
<pre cols="72">Thank you
Yours sincerely,
TAY wee-beng</pre>
<div>
<div class="h5">
<div>On 24/8/2015 5:54 PM, Timothée Nicolas wrote:<br>
</div>
<blockquote type="cite">
<div dir="ltr">
<div>
<div>
<div>
<div>
<div>Hi,<br>
<br>
</div>
ex5 of snes can give you an example of the
two routines.<br>
<br>
</div>
The C version ex5.c uses MatSetValuesStencil
whereas the Fortran90 version ex5f90.F uses
MatSetValuesLocal.<br>
<br>
</div>
However, I use MatSetValuesStencil also in
Fortran, there is no problem, and no need to
mess around with DMDAGetAO, I think.<br>
<br>
</div>
To input values in the matrix, you need to do
the following :<br>
<br>
</div>
! Declare the matstencils for matrix columns and
rows<br>
<div>MatStencil :: row(4,1),col(4,n)<br>
</div>
<div>! Declare the quantity which will store the
actual matrix elements<br>
</div>
<div>PetscScalar :: v(8)<br>
<br>
</div>
<div>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<br>
<br>
</div>
<div>Then you define the row number by naturally
doing the following, inside a local loop :<br>
<br>
row(MatStencil_i,1) = i -1<br>
row(MatStencil_j,1) = j -1<br>
row(MatStencil_c,1) = 1 -1<br>
<br>
</div>
<div>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.<br>
<br>
</div>
<div>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)<br>
<br>
col(MatStencil_i,1) = i -1<br>
col(MatStencil_j,1) = j -1<br>
col(MatStencil_c,1) = 1 -1<br>
v(1) = whatever_it_is<br>
<br>
col(MatStencil_i,2) = i-1 -1<br>
col(MatStencil_j,2) = j -1<br>
col(MatStencil_c,2) = 1 -1<br>
v(2) = whatever_it_is<br>
<br>
col(MatStencil_i,3) = i -1<br>
col(MatStencil_j,3) = j -1<br>
col(MatStencil_c,3) = 2 -1<br>
v(3) = whatever_it_is<br>
<br>
...<br>
...<br>
..<br>
<br>
...<br>
...<br>
...<br>
<br>
</div>
<div>Note that the index of the degree of freedom
(or what field you are coupling to), is
indicated by MatStencil_c<br>
<br>
</div>
<div><br>
</div>
<div>Finally use MatSetValuesStencil<br>
</div>
<div><br>
</div>
<div>ione = 1<br>
</div>
<div>isix = 6<br>
</div>
<div>call
MatSetValuesStencil(Matrix,ione,row,isix,col,v,INSERT_VALUES,ierr)<br>
<br>
</div>
<div>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.<br>
<br>
</div>
<div>Best<br>
<br>
</div>
<div>Timothee<br>
<br>
</div>
</div>
<div class="gmail_extra"><br>
<div class="gmail_quote">2015-08-24 18:09
GMT+09:00 Wee-Beng Tay <span dir="ltr"><<a
moz-do-not-send="true"
href="mailto:zonexo@gmail.com"
target="_blank"><a class="moz-txt-link-abbreviated" href="mailto:zonexo@gmail.com">zonexo@gmail.com</a></a>></span>:<br>
<blockquote class="gmail_quote" style="margin:0
0 0 .8ex;border-left:1px #ccc
solid;padding-left:1ex">
<div dir="ltr">Hi,
<div><br>
</div>
<div>I'm modifying my 3d fortran code from
MPI along 1 direction (z) to MPI along 2
directions (y,z)</div>
<div><br>
</div>
<div>Previously I was using MatSetValues
with global indices. However, now I'm
using DM and global indices is much more
difficult.</div>
<div><br>
</div>
<div>I come across MatSetValuesStencil or
MatSetValuesLocal.</div>
<div><br>
</div>
<div>So what's the difference bet the one
since they both seem to work locally?</div>
<div><br>
</div>
<div>Which is a simpler/better option?</div>
<div><br>
</div>
<div>Is there an example in Fortran for
MatSetValuesStencil?</div>
<div><br>
</div>
<div>Do I also need to use DMDAGetAO
together with MatSetValuesStencil or
MatSetValuesLocal?</div>
<div><br>
</div>
<div>Thanks!</div>
</div>
</blockquote>
</div>
<br>
</div>
</blockquote>
<br>
</div>
</div>
</div>
</blockquote>
</div>
<br>
</div>
</blockquote>
<br>
</body>
</html>