<div dir="ltr">On Mon, Nov 11, 2013 at 1:08 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<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"><div class="gmail_extra"><div class="gmail_quote"><div><div class="h5">On Mon, Nov 4, 2013 at 5:25 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<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>On Mon, Nov 4, 2013 at 4:54 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><div>On Mon, Nov 4, 2013 at 4:49 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:</div>
</div><div class="gmail_extra"><div class="gmail_quote"><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">On Fri, Nov 1, 2013 at 1:30 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><div>On Fri, Nov 1, 2013 at 12:08 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:</div>
</div><div class="gmail_extra"><div class="gmail_quote"><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">On Thu, Oct 31, 2013 at 11:34 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><div>On Thu, Oct 31, 2013 at 5:02 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:</div>
</div><div class="gmail_extra"><div class="gmail_quote"><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">On Thu, Oct 31, 2013 at 7:43 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, Oct 31, 2013 at 1:25 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:</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 class="gmail_extra"><div class="gmail_quote"><div>I did not call DMCreateMatrix() before when using just dmda (and it is working). In any case, I added a call to DMCreateMatrix() now but still there are problems. Here are the code snippets and the error I get:<br>
</div></div></div></div></blockquote><div><br></div></div><div>Then you were allowing it to be called automatically by PETSc, and it would have been this time as well.</div><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 class="gmail_extra"><div class="gmail_quote"><div></div><div>//Setting up the section and linking it to DM:<br>ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s);CHKERRQ(ierr);<br> ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);<br>
ierr = PetscSectionSetChart(s, 0, nC);CHKERRQ(ierr);<br>
for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {<br> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {<br> PetscInt point;<br> if(isPosInDomain(&testPoisson,i,j,0)) {<br>
ierr = DMDAGetCellPoint(dm, i, j, 0, &point);CHKERRQ(ierr);<br> ierr = PetscSectionSetDof(s, point, testPoisson.mDof); // No. of dofs associated with the point.<br> }<br> }<br>
}<br> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);<br> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);<br> ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);<br> ierr = DMCreateMatrix(dm,&A);CHKERRQ(ierr);<br>
<br></div><div> //Set up KSP:<br> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);<br> ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr);<br> ierr = KSPSetComputeOperators(ksp,computeMatrix2dSection,(void*)&testPoisson);CHKERRQ(ierr);<br>
ierr = KSPSetComputeRHS(ksp,computeRhs2dSection,(void*)&testPoisson);CHKERRQ(ierr);<br> ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);<br> ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);<br><br></div><div>------------------------------------------------------<br>
The computeMatrix2dSection function has:<br><br>ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);<br> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);<br> ierr = DMGetDefaultGlobalSection(da,&gs);CHKERRQ(ierr);<br>
for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {<br> for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {<br> if (isPosInDomain(ctx,i,j)) {<br> ierr = DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);<br>
ierr = PetscSectionGetOffset(gs,point,&row);<br> ierr = PetscSectionGetDof(gs,point,&dof);<br> for(PetscInt cDof = 0; cDof < dof; ++cDof) {<br> num = 0;<br>
row+=cDof;<br> col[num] = row; //(i,j) position<br> v[num++] = -4;<br> if(isPosInDomain(ctx,i+1,j)) {<br> ierr = DMDAGetCellPoint(da,i+1,j,0,&point);CHKERRQ(ierr);<br>
ierr = PetscSectionGetOffset(gs,point,&col[num]);<br> col[num] += cDof;<br> v[num++] = 1;<br> }<br> if(isPosInDomain(ctx,i-1,j)) {<br>
ierr = DMDAGetCellPoint(da,i-1,j,0,&point);CHKERRQ(ierr);<br> ierr = PetscSectionGetOffset(gs,point,&col[num]);<br> col[num] += cDof;<br> v[num++] = 1;<br>
}<br> if(isPosInDomain(ctx,i,j+1)) {<br> ierr = DMDAGetCellPoint(da,i,j+1,0,&point);CHKERRQ(ierr);<br> ierr = PetscSectionGetOffset(gs,point,&col[num]);<br>
col[num] += cDof;<br> v[num++] = 1;<br> }<br> if(isPosInDomain(ctx,i,j-1)) {<br> ierr = DMDAGetCellPoint(da,i,j-1,0,&point);CHKERRQ(ierr);<br>
ierr = PetscSectionGetOffset(gs,point,&col[num]);<br> col[num] += cDof;<br> v[num++] = 1;<br> }<br> ierr = MatSetValues(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);<br>
}<br> }<br> }<br> }<br> ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br> ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br></div><div><br></div><div>
But this is not working. For e.g. for a 6X6 grid size I get the following error:<br></div></div></div></div></blockquote><div><br></div></div></div><div>Okay, here is how to debug this:</div><div><br></div><div> 0) Go to a single scalar equations to make things easier</div>
<div><br></div><div> 1) Use MatView(A, PETSC_VIEWER_STDOUT_WORLD), which should output a 36 row matrix with</div><div> the preallocated nonzero pattern, filled with zeros</div><div><br></div><div> 2) Make sure this is the pattern you want</div>
<div><br></div><div> 3) Run in the debugger with -start_in_debugger</div><div><br></div><div> 4) When you get the error, see</div><div><br></div><div> a) If the (i, j) is one that should be setting a value</div><div>
<br></div><div> b) Why this (i, j) was not preallocated</div></div></div></div></blockquote><div><br></div><div>Up to 4 (a), it is correct. There is a problem in the way DMDAGetCellPoint and PetscSectionGetOffset works in this case scenario. I will try to explain it below for the case of 4X4 grid.<br>
</div><div>First case: <br>If I set the computational domain to be all the points of the dmda grid, (i.e. isPosInDomain(..,i,j,..) returns true for all i and j in the dmda grid), the program runs fine and does not give any error.<br>
</div><div><br>Second case:<br>I want the computational domain to be some part of the whole grid. There is a problem in this case.<br>The following test is in a case where,<br>isPosInDomain(..,i,j,..) returns true ONLY for (i,j) pairs of (2,1) (1,2) (2,2) (3,2) and (3,2). The grid with its corresponding point number in the petscsection is shown below:<br>
<br>12 13 14 15<br></div><div>8 9 10 11<br>4 5 6 7<br>0 1 2 3 <br><br>of which 6, 9, 10, 11 and 14 correspond to the (i,j) points that returns true for isPosInDomain(..,i,j,..) <br>MatView gives me the 16-row matrix with the star stencil non-zero structure as expected.<br>
The error I get is: new non-zero at (0,2) caused a malloc.<br></div><div><br>This error is for the value of (i,j) = (2,1), i.e. point 6.<br></div><div>The loop to set values in the matrix starts as:<br>for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {<br>
for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {<br> if (isPosInDomain(ctx,i,j)) {<br> ierr = DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);<br> ierr = PetscSectionGetOffset(gs,point,&row);<br>
<br></div><div>Now, what I wanted from the beginning is to create the matrix containing the rows corresponding to only those points (i,j) which has true values for isPosInDomain(ctx,i,j), that means in this case 5 row matrix. In the current test, the first (i,j) that reaches DMDAGetCellPoint is (2,1), which gives sets point variable to 6. <br>
The line PetscSectionGetOffset(gs,point,&row) sets the value of row to 0.<br></div>So I think this where the inconsistency lies. <br>MatSetValues(A,1,&row,num,col,v,INSERT_VALUES) expects row variable and col[] array to correspond to (i,j) based on DMDA, but PetsSectionGetOffset gives result based on how we masked away some of the points from dmda grid. Or am I missing something very basic here ?</div>
</div></div></blockquote><div><br></div></div></div><div>You are missing something basic. The row is correctly identified as 0, and that means the 2 is</div><div>point 10. </div></div></div></div></blockquote><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 class="gmail_extra"><div class="gmail_quote"><div>The problem must be preallocation. First, you should send the result of MatView().</div>
</div></div></div></blockquote><div><br></div><div>MatView shows that the preallocation and the matrix created is for the full dmda grid which is probably not what we want here, right ? Here is the matview result:<br></div>
</div></div></div></blockquote><div><br></div></div></div><div>Thats right. Did you call</div><div><br></div><div> DMSetDefaultSection(da, s)</div><div><br></div><div>after you made the new section and before you called DMGetMatrix()?</div>
</div></div></div></blockquote><div>Do you mean before I called DMCreateMatrix() ? I called DMSetDefaultSection(da,s) as below:<br><br>...<br></div><div> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);<br> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);<br>
ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);<br>// ierr = DMGetMatrix(dm,&A);CHKERRQ(ierr); //It seems DMGetMatrix is changed to DMCreateMatrix().<br> ierr = DMCreateMatrix(dm,&A);CHKERRQ(ierr);<br>
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br></div></div></div></div></blockquote><div><br></div></div></div><div>Yes, this is correct. I do exactly this in the example code and it works fine. I suspect</div>
<div>
you are not creating the section you think you are. You can send the code if you want</div><div>and I will look at it.</div></div></div></div></blockquote><div><br></div></div></div><div>Thanks, here is the code: <a href="https://github.com/bishesh/petscPoissonIrregular/blob/master/poissonIrregular.cxx" target="_blank">https://github.com/bishesh/petscPoissonIrregular/blob/master/poissonIrregular.cxx</a><br>
</div><div>It works for 2D and 3D without using PetscSection. (The choice of using Petsc_section is on line 68). Using PetsSection currently it implements only the 2D case with dof = 1. It works when the domain consists of all the cells of the dmda, it works (to check you can uncomment line 129 in isPosInDomain(). But it does not when I mask out some of the cells while creating the petscsection.<br>
</div></div></div></div></blockquote><div> </div></div></div><div>Hi Matt, I'm afraid if my last email somehow escaped you. Could you please let me know what could be possibly going wrong with my code linked above ?</div>
</div></div></div></blockquote><div><br></div><div>There is a bug in DMCreateMatrix() for DMDA when using a PetscSection. I am fixing it. I will mail when its done.</div><div><br></div><div> Thanks,</div><div><br></div><div>
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><span style="color:rgb(80,0,80)"> </span></div>
<div><div class="h5"><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"><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><br></div><div>
Matt</div><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 class="gmail_extra"><div class="gmail_quote"><div></div><div>This MatView on A shows (result below) that it still has 16 rows for all DMDA points instead of only those petscsection points for which I set a non-zero dof. Is it what it should not be doing ?<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><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>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">Mat Object: 1 MPI processes</span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)"> type: seqaij</span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 0: (0, 0) (1, 0) (4, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 1: (0, 0) (1, 0) (2, 0) (5, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 2: (1, 0) (2, 0) (3, 0) (6, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 3: (2, 0) (3, 0) (7, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 4: (0, 0) (4, 0) (5, 0) (8, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 5: (1, 0) (4, 0) (5, 0) (6, 0) (9, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 6: (2, 0) (5, 0) (6, 0) (7, 0) (10, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 7: (3, 0) (6, 0) (7, 0) (11, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 8: (4, 0) (8, 0) (9, 0) (12, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 9: (5, 0) (8, 0) (9, 0) (10, 0) (13, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 10: (6, 0) (9, 0) (10, 0) (11, 0) (14, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 11: (7, 0) (10, 0) (11, 0) (15, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 12: (8, 0) (12, 0) (13, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 13: (9, 0) (12, 0) (13, 0) (14, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 14: (10, 0) (13, 0) (14, 0) (15, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(26,26,26)">row 15: (11, 0) (14, 0) (15, 0) </span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(178,8,8)">[0]PETSC ERROR: --------------------- Error Message ------------------------------------</span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(178,8,8)">[0]PETSC ERROR: Argument out of range!</span></p>
<p style="margin:0px;text-indent:0px"><span style="font-family:'Monospace';font-size:9pt;color:rgb(178,8,8)">[0]PETSC ERROR: New nonzero at (0,2) caused a malloc!</span></p> <br></div><div><br></div></div></div>
</div></blockquote></div></div></div></div></div></blockquote><div><br></div><div>For reference here is the code snippet again:<br>...<br>ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s);CHKERRQ(ierr);<br> ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);<br>
ierr = PetscSectionSetChart(s, 0, nC);CHKERRQ(ierr);<br> for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {<br> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {<br> PetscInt point;<br>
if(isPosInDomain(&testPoisson,i,j,0)) {<br> ierr = DMDAGetCellPoint(dm, i, j, 0, &point);CHKERRQ(ierr);<br> ierr = PetscSectionSetDof(s, point, testPoisson.mDof); // No. of dofs associated with the point.<br>
} <br> }<br> }<br> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);<br> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);<br> ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);<br> ierr = DMCreateMatrix(dm,&A);CHKERRQ(ierr);<br>
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br><br>ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);<br> ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr);<br><br> ierr = KSPSetComputeOperators(ksp,computeMatrix2dSection,(void*)&testPoisson);CHKERRQ(ierr);<br>
ierr = KSPSetComputeRHS(ksp,computeRhs2dSection,(void*)&testPoisson);CHKERRQ(ierr);<br><br>....<br></div><div>The function computeMatrix2dSection:<br><br>#undef __FUNCT__<br>#define __FUNCT__ "computeMatrix2dSection"<br>
PetscErrorCode computeMatrix2dSection(KSP ksp, Mat A, Mat B, MatStructure *str, void *context) {<br> PoissonModel *ctx = (PoissonModel*)context;<br> PetscErrorCode ierr;<br> DM da;<br>
DMDALocalInfo info;<br> PetscInt row, col[5];<br> PetscInt dof;<br> PetscScalar v[5]; //array to store 5 point stencil for each row<br> PetscInt num; //non-zero position in the current row<br>
// PetscScalar dScale = 1; //to scale Dirichlet identity rows if needed<br> PetscSection gs; //Global Section that keeps the grid info and indices<br> PetscInt point; //Current point of the petscSection<br>
<br> PetscFunctionBeginUser;<br> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);<br> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);<br> ierr = DMGetDefaultGlobalSection(da,&gs);CHKERRQ(ierr);<br> ierr = DMCreateMatrix(da,&A);CHKERRQ(ierr);<br>
// ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br> for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {<br> for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {<br> if (isPosInDomain(ctx,i,j)) {<br>
ierr = DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);<br> ierr = PetscSectionGetOffset(gs,point,&row);<br> ierr = PetscSectionGetDof(gs,point,&dof);<br> for(PetscInt cDof = 0; cDof < ctx->mDof; ++cDof) {<br>
num = 0;<br> row+=cDof;<br> col[num] = row; //(i,j) position<br> v[num++] = -4;<br> if(isPosInDomain(ctx,i+1,j)) {<br>
ierr = DMDAGetCellPoint(da,i+1,j,0,&point);CHKERRQ(ierr);<br> ierr = PetscSectionGetOffset(gs,point,&col[num]);<br> col[num] += cDof;<br> v[num++] = 1;<br>
}<br> if(isPosInDomain(ctx,i-1,j)) {<br> ierr = DMDAGetCellPoint(da,i-1,j,0,&point);CHKERRQ(ierr);<br> ierr = PetscSectionGetOffset(gs,point,&col[num]);<br>
col[num] += cDof;<br> v[num++] = 1;<br> }<br> if(isPosInDomain(ctx,i,j+1)) {<br> ierr = DMDAGetCellPoint(da,i,j+1,0,&point);CHKERRQ(ierr);<br>
ierr = PetscSectionGetOffset(gs,point,&col[num]);<br> col[num] += cDof;<br> v[num++] = 1;<br> }<br> if(isPosInDomain(ctx,i,j-1)) {<br>
ierr = DMDAGetCellPoint(da,i,j-1,0,&point);CHKERRQ(ierr);<br> ierr = PetscSectionGetOffset(gs,point,&col[num]);<br> col[num] += cDof;<br> v[num++] = 1;<br>
}<br> ierr = MatSetValues(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);<br> }<br> }<br> }<br> }<br> ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br><br></div><div><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><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><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><br></div><div> Matt</div><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 class="gmail_extra">
<div class="gmail_quote"><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><br></div><div> Thanks,</div><div><br></div><div> Matt</div><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 class="gmail_extra"><div class="gmail_quote"><div>[0]PETSC ERROR: --------------------- Error Message ------------------------------------<br>[0]PETSC ERROR: Argument out of range!<br>[0]PETSC ERROR: New nonzero at (0,6) caused a malloc!<br>
[0]PETSC ERROR: ------------------------------------------------------------------------<br>
[0]PETSC ERROR: Petsc Development GIT revision: 61e5e40bb2c5bf2423e94b71a15fef47e411b0da GIT Date: 2013-10-25 21:47:45 -0500<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: build/poissonIrregular on a arch-linux2-cxx-debug named edwards by bkhanal Thu Oct 31 19:23:58 2013<br>
[0]PETSC ERROR: Libraries linked from /home/bkhanal/Documents/softwares/petsc/arch-linux2-cxx-debug/lib<br>[0]PETSC ERROR: Configure run at Sat Oct 26 16:35:15 2013<br>[0]PETSC ERROR: Configure options --download-mpich -download-f-blas-lapack=1 --download-metis --download-parmetis --download-superlu_dist --download-scalapack --download-mumps --download-hypre --with-clanguage=cxx<br>
[0]PETSC ERROR: ------------------------------------------------------------------------<br>[0]PETSC ERROR: MatSetValues_SeqAIJ() line 413 in src/mat/impls/aij/seq/aij.c<br>[0]PETSC ERROR: MatSetValues() line 1130 in src/mat/interface/matrix.c<br>
[0]PETSC ERROR: computeMatrix2dSection() line 318 in /user/bkhanal/home/works/poissonIrregular/poissonIrregular.cxx<br>[0]PETSC ERROR: KSPSetUp() line 228 in src/ksp/ksp/interface/itfunc.c<br>[0]PETSC ERROR: KSPSolve() line 399 in src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: main() line 598 in /user/bkhanal/home/works/poissonIrregular/poissonIrregular.cxx<br>application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0<br>[unset]: aborting job:<br>application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0<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><br></div><div> Matt</div><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 class="gmail_extra"><div class="gmail_quote"><div>However, now since I do not want the rows in the matrix corresponding to the points in the grid which do not lie in the computational domain. This means the size of the matrix will not necessarily equal the total number of cells in DMDA. So in the following code:<br>
<br> ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);<br> ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s);CHKERRQ(ierr);<br> ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);<br> ierr = PetscSectionSetChart(s, 0, nC);CHKERRQ(ierr);<br>
for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {<br> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {<br> PetscInt point;<br> if(isPosInDomain(&testPoisson,i,j,0)) {<br>
ierr = DMDAGetCellPoint(dm, i, j, 0, &point);CHKERRQ(ierr);<br> ierr = PetscSectionSetDof(s, point, testPoisson.mDof); // No. of dofs associated with the point.<br> }<br><br>
}<br>
}<br> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);<br> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);<br> ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);<br><br></div><div>should I myself compute proper nC (i.e. total no. of points for which I want the matrix to have a corresponding row) before calling PetscSectionSetChart() or is the masking out of the points inside the loop with if(isPosInDomain(&testPoisson,i,j,0)) sufficient ?<br>
</div><div>And, when you write:<br> PetscSectionGetOffset(gs, p, &ind);<div> row = ind < 0 ? -(ind+1) : ind;</div>
<div><br></div><div>it seems you allow the possibility of getting a -ve ind when using PetscSectionGetOffset. When would an offset to a particular dof and point?<br><br></div><br></div><div><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> </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> for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {<br> for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {<br> num = 0;<br>
/*ierr = DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);*/ /*Get the PetscPoint*/<br>
</div><div> /*But I did now understand how I would emulate the row.c and col.c info with PetscSection*/<br></div><div> row.i = i; row.j = j;<br></div></div></div></div></blockquote><div><br>
</div></div><div>Here you would call</div><div><br></div><div> ierr = DMDAGetCellPoint(da, i, j, 0, &cell);CHKERRQ(ierr);</div><div> ierr = PetscSectionGetOffset(da, cell, &row);CHKERRQ(ierr);</div><div> ierr = PetscSectionGetDof(da, cell, &dof);CHKERRQ(ierr);</div>
<div><br></div><div>This means that this cell has global rows = [row, row+dof).</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 class="gmail_extra">
<div class="gmail_quote"><div> col[num].i = i; col[num].j = j;<br>
if (isPosInDomain(ctx,i,j)) { /*put the operator for only the values for only the points lying in the domain*/<br> v[num++] = -4;<br> if(isPosInDomain(ctx,i+1,j)) {<br>
col[num].i = i+1; col[num].j = j;<br> v[num++] = 1;<br> }<br> if(isPosInDomain(ctx,i-1,j)) {<br> col[num].i = i-1; col[num].j = j;<br>
v[num++] = 1;<br> }<br> if(isPosInDomain(ctx,i,j+1)) {<br> col[num].i = i; col[num].j = j+1;<br> v[num++] = 1;<br>
}<br> if(isPosInDomain(ctx,i,j-1)) {<br> col[num].i = i; col[num].j = j-1;<br> v[num++] = 1;<br> }<br> } else {<br>
v[num++] = dScale; /*set Dirichlet identity rows for points in the rectangle but outside the computational domain*/<br> }<br></div></div></div></div></blockquote><div><br></div></div><div>
You do the same thing for the columns, and then call</div>
<div><br></div><div> ierr = MatSetValues(A, dof, rows, num*dof, cols, v, INSERT_VALUES);CHKERRQ(ierr);</div><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 class="gmail_extra"><div class="gmail_quote"><div> ierr = MatSetValuesStencil(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);<br>
}<br> }<br> }<br></div><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 class="gmail_extra">
<div class="gmail_quote">
<div></div><div> Thanks,</div><div><br></div><div> Matt</div><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 class="gmail_extra">
<div class="gmail_quote"><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><br></div><div> Matt</div>
<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 class="gmail_extra"><div class="gmail_quote"><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><br></div><div> Matt</div><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 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 class="gmail_extra"><div class="gmail_quote">
<div><br></div><div> Thanks,</div><div><br></div><div> Matt </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><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><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"><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> If not, just put the identity into</div>
<div>the rows you do not use on the full cube. It will not hurt scalability or convergence.</div></div></div></div></blockquote><div> </div>In the case of Poisson with Dirichlet condition this might be the case. But is it always true that having identity rows in the system matrix will not hurt convergence ? I thought otherwise for the following reasons:<br>
</div><div class="gmail_quote">1) Having read Jed's answer here : <a href="http://scicomp.stackexchange.com/questions/3426/why-is-pinning-a-point-to-remove-a-null-space-bad/3427#3427" target="_blank">http://scicomp.stackexchange.com/questions/3426/why-is-pinning-a-point-to-remove-a-null-space-bad/3427#3427</a></div>
</div></div></blockquote><div><br></div></div><div>Jed is talking about a constraint on a the pressure at a point. This is just decoupling these unknowns from the rest</div><div>of the problem.</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 class="gmail_extra"><div class="gmail_quote">2) Some observation I am getting (but I am still doing more experiments to confirm) while solving my staggered-grid 3D stokes flow with schur complement and using -pc_type gamg for A00 matrix. Putting the identity rows for dirichlet boundaries and for ghost cells seemed to have effects on its convergence. I'm hoping once I know how to use PetscSection, I can get rid of using ghost cells method for the staggered grid and get rid of the identity rows too.<br>
</div></div></div></blockquote><div><br></div></div><div>It can change the exact iteration, but it does not make the matrix conditioning worse.</div><div><br></div><div> Matt</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 class="gmail_extra"><div class="gmail_quote">
</div><div class="gmail_quote">Anyway please provide me with some pointers so that I can start trying with petscsection on top of a dmda, in the beginning for non-staggered case.<br><br></div><div class="gmail_quote">Thanks,<br>
Bishesh<br></div><div class="gmail_quote"><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><br></div><div> Matt</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>Thanks,<br>Bishesh<span><font color="#888888"><br></font></span></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br><br clear="all"><span><font color="#888888"><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></font></span></div></div>
</blockquote></div><br></div></div>
</blockquote></div></div><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></div>
</blockquote></div></div><br></div></div>
</blockquote></div></div><div><br><br clear="all"><span><font color="#888888"><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></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><div><br><br clear="all"><span><font color="#888888"><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></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><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></font></span></div></div>
</blockquote></div><br></div></div>
</blockquote></div></div></div><div><div><br><br clear="all"><span><font color="#888888"><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></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><div><br><br clear="all"><span><font color="#888888"><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></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><div><br><br clear="all"><span><font color="#888888"><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></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><div><br><br clear="all"><span><font color="#888888"><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></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><div><br><br clear="all"><span><font color="#888888"><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></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br></font></span></div></div><span><font color="#888888">
</font></span></blockquote></div></div></div><span><font color="#888888"><div><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></font></span></div></div>
</blockquote></div><br></div></div>
</blockquote></div></div></div><div><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></div></div>
</blockquote></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>