<div dir="ltr">On Thu, Jul 4, 2013 at 10:39 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br><div class="gmail_extra"><div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote"><div><div class="h5">On Thu, Jul 4, 2013 at 3:11 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@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="gmail_extra"><br><div class="gmail_quote"><div>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>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>
<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><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 */<div>
<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></div><div><div> for(i = info.ys; i < info.ys+info.ym; ++i){<br> for(j = info.xs; j < info.xs+info.xm; ++j){<br></div></div>
</div></div></div></div></blockquote><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 class="gmail_extra"><div class="gmail_quote">
<div><div><div> 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></div></div></div></div></div></div></blockquote><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 class="gmail_extra"><div class="gmail_quote"><div><div><div> 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></div></div></div></div></div></div></blockquote></div></div><div><br><div>It seems to me that the problem is caused because of this kind of
assignment col[1].j = j+1 when j takes the value "info.xs+info.xm-1" and
for similar places below with "i+1" when "i = info.ys+info.ym".<br></div>My
guess is that the each processor contains only the local portion of the
matrix "mat1" corresponding to the grid area it is taking care of. It
probably does not store the ghost values. But when for e.g.<br></div><div>j = info.xs+info.xm-1, we have col[0].j = j+1 and MatSetValuesStencil needs accessing the region of another processor. Am I thinking correctly ? Iis this the case where I need each processor allocating a space for ghost values too?<br>
</div><div>If so, how do I do it?</div></div></div></div></blockquote><div><br></div><div style>If you want to step outside the mesh, you either need DMDA_BOUNDARY_PERIODIC or DMDA_BOUNDARY_GHOSTED.</div><div style><br></div>
<div style> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><div class="h5">
<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><div> 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></div></div> 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><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 class="gmail_extra">
<div class="gmail_quote">
<div><div><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><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></div></div><br></div></div>
</blockquote></div></div></div><br></div></div>
</blockquote></div><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
</div></div>