<div dir="ltr"><div><div>I don't really understand what you say, but it does not sound right.<br><br></div>You can enter the boundary points separately and then the points outside the boundary on separate calls, like this :<br><br>  do j=user%ys,user%ye<br>     do i=user%xs,user%xe<br>        if (i.eq.1 .or. i.eq.user%mx .or. j .eq. 1 .or. j .eq. user%my) then<br></div>        ! boundary point <br><div>           row(MatStencil_i,1) = i          -1<br>           row(MatStencil_j,1) = j          -1<br>           row(MatStencil_c,1) = 1          -1<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)   = one<br>           call MatSetValuesStencil(jac_prec,ione,row,ione,col,v,          &<br>                &                           INSERT_VALUES,ierr)<br><br></div><div>         else<br>           row(MatStencil_i,1) = i          -1<br>           row(MatStencil_j,1) = j          -1<br>           row(MatStencil_c,1) = 1          -1<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) = undemi*dxm1*(vx_ip1j-vx_im1j) + two*user%nu*(dxm1**2+dym1**2)<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) = undemi*dxm1*(vx_ij-vx_ip1j) - user%nu*dxm1**2<br>           col(MatStencil_i,3) = i-1       -1<br>           col(MatStencil_j,3) = j         -1<br>           col(MatStencil_c,3) = 1         -1<br>           v(3) = -undemi*dxm1*(vx_ij-vx_im1j)  - user%nu*dxm1**2<br>           col(MatStencil_i,4) = i         -1<br>           col(MatStencil_j,4) = j+1       -1<br>           col(MatStencil_c,4) = 1         -1<br>           v(4) = undemi*dym1*vy_ij           - user%nu*dym1**2<br>           col(MatStencil_i,5) = i         -1<br>           col(MatStencil_j,5) = j-1       -1<br>           col(MatStencil_c,5) = 1         -1<br></div><div>           v(5) =  -undemi*dym1*vy_ij           - user%nu*dym1**2<br>           col(MatStencil_i,6) = i         -1<br>           col(MatStencil_j,6) = j         -1<br>           col(MatStencil_c,6) = 2         -1<br>           v(6) = undemi*dym1*(vx_ijp1-vx_ijm1)<br>           col(MatStencil_i,7) = i+1       -1<br>           col(MatStencil_j,7) = j         -1<br>           col(MatStencil_c,7) = 2         -1<br>           v(7) = -undemi*dxm1*vy_ip1j<br>           col(MatStencil_i,8) = i-1       -1<br>           col(MatStencil_j,8) = j         -1<br>           col(MatStencil_c,8) = 2         -1<br>           v(8) = undemi*dxm1*vy_im1j<br><br>          call MatSetValuesStencil(jac_prec,ione,row,ieight,col,v,         &<br>                &                                INSERT_VALUES,ierr)<br><br>       endif<br>     enddo<br>  enddo<br><br></div><div>Timothee<br></div><div><br><br><br><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-08-26 14:02 GMT+09:00 TAY wee-beng <span dir="ltr"><<a 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 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?<span class=""><br>
    <pre cols="72">Thank you

Yours sincerely,

TAY wee-beng</pre>
    </span><div><div class="h5"><div>On 26/8/2015 12:24 PM, Timothée Nicolas
      wrote:<br>
    </div>
    <blockquote 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 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>
                  <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 href="mailto:zonexo@gmail.com" target="_blank"></a><a 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 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>
  </div></div></div>

</blockquote></div><br></div>