<div dir="ltr">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 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>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>Then you were allowing it to be called automatically by PETSc, and it would have been this time as well.</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><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>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><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>[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 class="HOEnZb"><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 class="HOEnZb"><font color="#888888">
</font></span></blockquote></div><span class="HOEnZb"><font color="#888888"><br></font></span></div></div><span class="HOEnZb"><font color="#888888">
</font></span></blockquote></div></div></div><span class="HOEnZb"><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><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>