<div dir="ltr"><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Jul 4, 2013 at 2:59 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="im">On Thu, Jul 4, 2013 at 7:32 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</div><div class="gmail_extra"><div class="gmail_quote"><div class="im">
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi all,<br>I'm trying to use DMCreateMatrix and MatStencil to fill up a matrix that results from a finite difference discretization of a PDE (objective is to solve the resulting linear system in the form of Ax=b).<br>


<br>However there are some errors and I'm not sure about the issue!<br>A short demo:<br>A 2D mXn grid, with two variables at each node (dof=2),<br>so the resulting A would be 2mn X 2mn.<br>Let's say the two variables are vx and vy, and the two associated <br>


equations discretized are x-eq and y-eq.<br><br>Here is the relevant part of the code: (mat1 variable for the A matrix):<br></div></div></blockquote><div><br></div></div><div>Just sending the whole code is better. I suspect the logic is wrong for selecting the boundary.</div>

<div><br></div><div>   Matt</div></div></div></div></blockquote><div><br></div><div>Thanks, here is the complete code:<br><br>#include <petscsys.h><br>#include <petscdmda.h><br><br>#undef __FUNCT__<br>#define __FUNCT__ "main"<br>
int main(int argc,char **argv)<br>{<br>    PetscErrorCode ierr;<br>    DM             da;<br>    PetscInt       m=10,n=10;<br><br>    ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);<br><br>    /* Create a DMDA and an associated vector */<br>
    ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,<br>                        PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);<br><br>    Mat mat1;<br>    ierr = DMCreateMatrix(da,MATMPIAIJ,&mat1); CHKERRQ(ierr);<br>
    MatStencil row, col[4]; //let's test 4 non-zeros in one row.<br>    PetscScalar val[4];<br>    PetscScalar coeff = 2.; //just a constant coefficient for testing.<br>    PetscInt i,j;<br>    DMDALocalInfo info;<br>
<br>    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);<br><br>    for(i = info.ys; i < info.ys+info.ym; ++i){<br>        for(j = info.xs; j < info.xs+info.xm; ++j){<br>            row.i = i;  row.j = j;   //one node at a time<br>
            if (j == 0 || j == info.mx-1){ //left and right borders<br>                //vx(i,j) = 0;<br>                row.c = 0;<br>                col[0].c = 0;   val[0] = 1;<br>                col[0].i = i;   col[0].j = j;<br>
                ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br><br>                //vy: //vy(i,j) - c*vy(i,j+-1) = 0;<br>                row.c = 1;<br>                col[0].c = 1;   val[1] = coeff;<br>
                col[1].c = 1;<br>                col[1].i = i;<br>                if(j == 0) //vy(i,j) - c*vy(i,j+1) = 0;<br>                    col[1].j = j+1;<br>                else //vy(i,j) - c*vy(i,j-1) = 0;<br>                    col[1].j = j-1;<br>
<br>                ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>            }<br>            else if (i == 0 || i == info.my-1){  //top and bottom borders<br>                //vx: vx(i,j) - c* vx(i+-1,j) = 0;<br>
                row.c = 0;<br>                col[0].c = 0;   val[0] = 1;<br>                col[0].i = i;   col[0].j = j;<br>                col[1].c = 0;   val[1] = coeff;<br>                if (i == 0) //vx(i,j) - c*vx(i+1,j) = 0;<br>
                    col[1].i = i+1;<br>                else //vx(i,j) - c*vx(i-1,j) = 0;<br>                    col[1].i = i-1;<br>                col[1].j = j;<br>                ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
<br>                //vy(i,j) = 0;<br>                row.c = 1;<br>                col[0].c = 1;<br>                ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>            }<br>            else { //Interior points:<br>
                row.c = 0;//x-eq<br>                col[0].c = 0;   val[0] = 2*coeff;<br>                col[0].i = i;   col[0].j = j+1;<br>                col[1].c = 0;   val[1] = -val[0] - coeff;<br>                col[1].i = i;   col[1].j = j;<br>
                col[2].c = 1;   val[2] = 4*coeff;<br>                col[2].i = i+1; col[2].j = j;<br>                col[3].c = 1;   val[3] = 4*coeff;<br>                col[3].i = i;   col[3].j = j;<br>                ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
<br>                row.c = 1; //y-eq<br>                col[0].c = 1;   val[0] = 2*coeff;<br>                col[0].i = i;   col[0].j = j+1;<br>                col[1].c = 1;   val[1] = -val[0] - coeff;<br>                col[1].i = i;   col[1].j = j;<br>
                col[2].c = 0;   val[2] = 4*coeff;<br>                col[2].i = i+1; col[2].j = j;<br>                col[3].c = 0;   val[3] = 4*coeff;<br>                col[3].i = i;   col[3].j = j;<br>                ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
            }<br>        }<br>    }<br><br>    MatAssemblyBegin(mat1,MAT_FINAL_ASSEMBLY);<br>    MatAssemblyEnd(mat1,MAT_FINAL_ASSEMBLY);<br>    ierr = MatSetOption(mat1,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);<br>
    ierr = PetscObjectSetName((PetscObject)mat1,"mat1");CHKERRQ(ierr);<br><br>    /* clean up and exit */<br>    ierr = DMDestroy(&da);CHKERRQ(ierr);<br>    ierr = MatDestroy(&mat1);<br>    ierr = PetscFinalize();<br>
    return 0;<br>}<br><br><br> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">
<div><div class="h5"><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div></div><div>    PetscInt m = 10, n=10;<br></div>

<div>     ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE,         DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,<br>
                        PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);<br><br>    Mat mat1;<br>    ierr = DMCreateMatrix(da,MATMPIAIJ,&mat1); CHKERRQ(ierr);<br>    MatStencil row, col[4]; //let's test 4 non-zeros in one row.<br>


    PetscScalar val[4];<br>    PetscScalar coeff = 2.; //just a constant coefficient for testing.<br>    PetscInt i,j;<br>    DMDALocalInfo info;<br>    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);<br>//Now fill up the matrix:<br>


for(i = info.ys; i < info.ys+info.ym; ++i){<br>        for(j = info.xs; j < info.xs+info.xm; ++j){<br>            row.i = i;  row.j = j;   //one node at a time<br>            if (j == 0 || j == info.mx-1){ //left and right borders<br>


                //vx(i,j) = 0;<br>                row.c = 0;<br>                col[0].c = 0;   val[0] = 1;<br>                col[0].i = i;   col[0].j = j;<br>                ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>


<br>                //vy: //vy(i,j) - c*vy(i,j+-1) = 0;<br>                row.c = 1;<br>                col[0].c = 1;   val[1] = coeff;<br>                col[1].c = 1;<br>                col[1].i = i;<br>                if(j == 0) //vy(i,j) - c*vy(i,j+1) = 0;<br>


                    col[1].j = j+1;<br>                else //vy(i,j) - c*vy(i,j-1) = 0;<br>                    col[1].j = j-1;<br><br>                ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>


            }<br>            else if (i == 0 || i == info.my-1){  //top and bottom borders<br>                //vx: vx(i,j) - c* vx(i+-1,j) = 0;<br>                row.c = 0;<br>                col[0].c = 0;   val[0] = 1;<br>


                col[0].i = i;   col[0].j = j;<br>                col[1].c = 0;   val[1] = coeff;<br>                if (i == 0) //vx(i,j) - c*vx(i+1,j) = 0;<br>                    col[1].i = i+1;<br>                else //vx(i,j) - c*vx(i-1,j) = 0;<br>


                    col[1].i = i-1;<br>                col[1].j = j;<br>                ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br><br>                //vy(i,j) = 0;<br>                row.c = 1;<br>


                col[0].c = 1;<br>                ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>            }<br>            else { //Interior points:<br>                row.c = 0;//x-eq<br>


                col[0].c = 0;   val[0] = 2*coeff;<br>                col[0].i = i;   col[0].j = j+1;<br>                col[1].c = 0;   val[1] = -val[0] - coeff;<br>                col[1].i = i;   col[1].j = j;<br>                col[2].c = 1;   val[2] = 4*coeff;<br>


                col[2].i = i+1; col[2].j = j;<br>                col[3].c = 1;   val[3] = 4*coeff;<br>                col[3].i = i;   col[3].j = j;<br>                ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>


<br>                row.c = 1; //y-eq<br>                col[0].c = 1;   val[0] = 2*coeff;<br>                col[0].i = i;   col[0].j = j+1;<br>                col[1].c = 1;   val[1] = -val[0] - coeff;<br>                col[1].i = i;   col[1].j = j;<br>


                col[2].c = 0;   val[2] = 4*coeff;<br>                col[2].i = i+1; col[2].j = j;<br>                col[3].c = 0;   val[3] = 4*coeff;<br>                col[3].i = i;   col[3].j = j;<br>                ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>


            }<br>        }<br>    }<br><br>    MatAssemblyBegin(mat1,MAT_FINAL_ASSEMBLY);<br>    MatAssemblyEnd(mat1,MAT_FINAL_ASSEMBLY);<br><br>However I get the following error when run with 2 processors.<br>[0]PETSC ERROR: --------------------- Error Message ------------------------------------<br>


[0]PETSC ERROR: Argument out of range!<br>[0]PETSC ERROR: Local index 120 too large 120 (max) at 0!<br>[0]PETSC ERROR: ------------------------------------------------------------------------<br>[0]PETSC ERROR: Petsc Release Version 3.4.1, Jun, 10, 2013 <br>


[0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>[0]PETSC ERROR: See docs/index.html for manual pages.<br>[0]PETSC ERROR: ------------------------------------------------------------------------<br>


[0]PETSC ERROR: examples/FDdmda_test on a arch-linux2-cxx-debug named edwards by bkhanal Thu Jul  4 14:20:11 2013<br>[0]PETSC ERROR: Libraries linked from /home/bkhanal/Documents/softwares/petsc-3.4.1/arch-linux2-cxx-debug/lib<br>


[0]PETSC ERROR: Configure run at Wed Jun 19 11:04:51 2013<br>[0]PETSC ERROR: Configure options --with-cc=gcc --with-fc=g77 --with-cxx=g++ --download-f-blas-lapack=1 --download-mpich=1 -with-clanguage=cxx --download-hypre=1<br>


[0]PETSC ERROR: ------------------------------------------------------------------------<br>[0]PETSC ERROR: ISLocalToGlobalMappingApply() line 444 in /home/bkhanal/Documents/softwares/petsc-3.4.1/src/vec/is/utils/isltog.c<br>


[0]PETSC ERROR: MatSetValuesLocal() line 1967 in /home/bkhanal/Documents/softwares/petsc-3.4.1/src/mat/interface/matrix.c<br>[0]PETSC ERROR: MatSetValuesStencil() line 1339 in /home/bkhanal/Documents/softwares/petsc-3.4.1/src/mat/interface/matrix.c<br>


[0]PETSC ERROR: main() line 75 in "unknowndirectory/"/user/bkhanal/home/works/cmake_tuts/petsc_test/examples/FDdmda_test.cxx<br>application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0<br>[cli_0]: aborting job:<br>


application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0<br><br>===================================================================================<br>=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES<br>=   EXIT CODE: 63<br>


=   CLEANING UP REMAINING PROCESSES<br>=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES<br>===================================================================================<br><br>I do not understand how the index (i,j) are out of range when set to corresponding fields in row and col variables. What could be the possible problem ? And any suggestions on the way to debug this sort of issues ?<br>


<br>Thanks,<br>Bishesh <br></div><br><br><div><br><br></div></div>
</blockquote></div></div></div><span class=""><font color="#888888"><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>

-- Norbert Wiener
</font></span></div></div>
</blockquote></div><br></div></div>