<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Oct 28, 2013 at 9:00 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 class="h5">On Mon, Oct 28, 2013 at 2:51 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</div></div><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"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Oct 28, 2013 at 6:19 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, Oct 28, 2013 at 10:13 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</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"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Oct 28, 2013 at 3:49 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, Oct 28, 2013 at 9:06 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</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"><br><div class="gmail_extra"><div class="gmail_quote"><div><div>On Mon, Oct 28, 2013 at 1:40 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, Oct 28, 2013 at 5:30 AM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</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"><br><div class="gmail_extra"><br><br><div class="gmail_quote">
<div>On Fri, Oct 25, 2013 at 10:21 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<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><div>On Fri, Oct 25, 2013 at 2:55 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 Fri, Oct 25, 2013 at 8:18 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, Oct 25, 2013 at 12:09 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
</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><div><div>Dear all,<br></div>I would like to know if some of the petsc objects that I have not used so far (IS, DMPlex, PetscSection) could be useful in the following case (of irregular domains): <br>
<br>
</div>Let's say that I have a 3D binary image (a cube).<br>The binary information of the image partitions the cube into a computational domain and non-computational domain. <br>I must solve a pde (say a Poisson equation) only on the computational domains (e.g: two isolated spheres within the cube). I'm using finite difference and say a dirichlet boundary condition<br>
</div><div><br>I know that I can create a dmda that will let me access the information from this 3D binary image, get all the coefficients, rhs values etc using the natural indexing (i,j,k).<br><br>Now, I would like to create a matrix corresponding to the laplace
operator (e.g. with standard 7 pt. stencil), and the corresponding RHS
that takes care of the dirchlet values too. <br></div><div>But in this matrix it should have the rows corresponding to the nodes only on the computational domain. It would be nice if I can easily (using (i,j,k) indexing) put on the rhs dirichlet values corresponding to the boundary points. <br>
Then, once the system is solved, put the values of the solution back to the corresponding positions in the binary image.<br></div><div>Later, I might have to extend this for the staggered grid case too.<br></div><div>So is petscsection or dmplex suitable for this so that I can set up the matrix with something like DMCreateMatrix ? Or what would you suggest as a suitable approach to this problem ?<br>
</div><div><br></div><div>I have looked at the manual and that led me to search for a simpler examples in petsc src directories. But most of the ones I encountered are with FEM (and I'm not familiar at all with FEM, so these examples serve more as a distraction with FEM jargon!) <br>
</div></div></blockquote><div><br></div></div></div><div>It sounds like the right solution for this is to use PetscSection on top of DMDA. I am working on this, but it is really</div><div>alpha code. If you feel comfortable with that level of development, we can help you.</div>
</div></div></div></blockquote><div> </div><div>Thanks, with the (short) experience of using Petsc so far and being familiar with the awesomeness (quick and helpful replies) of this mailing list, I would like to give it a try. Please give me some pointers to get going for the example case I mentioned above. A simple example of using PetscSection along with DMDA for finite volume (No FEM) would be great I think.<br>
Just a note: I'm currently using the petsc3.4.3 and have not used the development version before. </div></div></div></div></blockquote><div><br></div></div><div>Okay,</div><div><br></div><div>1) clone the repository using Git and build the 'next' branch.</div>
<div><br></div></div><div><div>2) then we will need to create a PetscSection that puts unknowns where you want them</div><div><br></div><div>3) Setup the solver as usual</div><div><br></div><div>You can do 1) an 3) before we do 2).</div>
<div>
<div><br></div></div></div></div></div></div></blockquote><div>I've done 1) and 3). I have one .cxx file that solves the system using DMDA (putting identity into the rows corresponding to the cells that are not used). <br>
Please let me know what I should do now.</div></div></div></div></blockquote><div><br></div></div></div><div>Okay, now write a loop to setup the PetscSection. I have the DMDA partitioning cells, so you would have</div><div>
unknowns in cells. The code should look like this:</div>
<div><br></div><div>PetscSectionCreate(comm, &s);</div><div><div>DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);<br></div></div><div>PetscSectionSetChart(s, 0, nC);</div><div>for (k = zs; k < zs+zm; ++k) {</div><div>
for (j = ys; j < ys+ym; ++j) {</div><div> for (i = xs; i < xs+xm; ++i) {</div><div> PetscInt point;</div><div><br></div><div> DMDAGetCellPoint(dm, i, j, k, &point);</div><div> PetscSectionSetDof(s, point, dof); // You know how many dof are on each vertex</div>
<div> }</div><div> }</div><div>}</div><div>PetscSectionSetUp(s);</div><div>DMSetDefaultSection(dm, s);</div><div>PetscSectionDestroy(&s);</div><div><br></div><div>I will merge the necessary stuff into 'next' to make this work.</div>
</div></div></div></blockquote><div><br></div></div></div><div>I have put an example without the PetscSection here: <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>From the code you have written above, how do we let PetscSection select only those cells that lie in the computational domain ? Is it that for every "point" inside the above loop, we check whether it lies in the domain or not, e.g using the function isPosInDomain(...) in the .cxx file I linked to?<br>
</div></div></div></div></blockquote><div><br></div></div></div><div>1) Give me permission to comment on the source (I am 'knepley')</div><div><br></div><div>2) You mask out the (i,j,k) that you do not want in that loop</div>
</div></div></div></blockquote><div><br></div><div>Done. <br>I mask it out using isPosInDomain() function::<br><div> <span>if</span><span>(</span><span>isPosInDomain</span><span>(</span><span>&</span><span>testPoisson</span><span>,</span><span>i</span><span>,</span><span>j</span><span>,</span><span>k</span><span>))</span> <span>{</span></div>
<div> <span>ierr</span> <span>=</span> <span>DMDAGetCellPoint</span><span>(</span><span>dm</span><span>,</span> <span>i</span><span>,</span> <span>j</span><span>,</span> <span>k</span><span>,</span> <span>&</span><span>point</span><span>);</span><span>CHKERRQ</span><span>(</span><span>ierr</span><span>);</span></div>
<div> <span>ierr</span> <span>=</span> <span>PetscSectionSetDof</span><span>(</span><span>s</span><span>,</span> <span>point</span><span>,</span> <span>testPoisson</span><span>.</span><span>mDof</span><span>);</span> <span>// You know how many dof are on each vertex </span></div>
<div> <span>}</span></div> <br></div><div>And please let me know when I can rebuild the 'next' branch of petsc so that DMDAGetCellPoint can be used. Currently compiler complains as it cannot find it. <br>
</div></div></div></div></blockquote><div><br></div></div></div><div>Pushed.</div></div></div></div></blockquote><div><br></div><div>Pulled it thanks, but compiler was complaining that it didn't find DMDAGetCellPoint. </div>
Doing grep -R 'DMDAGetCellPoint' include/<br><div>returned nothing although the implementation is there in dalocal.c <br> So I added the declaration to petscdmda.h file just above the declaration of DMDAGetNumCells. <br>
</div></div></div></div></blockquote><div><br></div></div></div><div>I will add the declaration.</div><div class="im"><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>I have also added you as a collaborator in the github project repository so that you can edit and add comments to the file I linked above. My computeMatrix still uses the DMDA instead of the PetscSection, so is it where the changes need to be made ?<br>
</div></div></div></div></blockquote><div><br></div></div><div>Instead of MatSetValuesStencil(), you would go back to using MatSetValues(), and calculate the indices using PetscSection.</div><div>You get the section</div>
<div><br>
</div><div> DMGetDefaultGlobalSection(dm, &gs)</div><div><br></div><div>and then</div><div><br></div><div> DMDAGetCellPoint(dm, i, j, k, &p);</div><div> PetscSectionGetOffset(gs, p, &ind);</div><div> row = ind < 0 ? -(ind+1) : ind;</div>
<div><br></div></div></div></div></blockquote><div><br><br></div><div>I'm sorry but did not understood. For e.g., in 2D, for each point p, I have 2 dof. Does PetsSectionGetOffset give me the two offsets for the dof sth like ind[0] = 0 and ind[1] = 1 ?<br>
</div><div>How would you change the following code that sets 2D laplacian stencil for 2 dof to use petscsection and Matsetvalues instead ?<br><br> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);<br> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);<br>
ierr = DMGetDefaultGlobalSection(dm,&gs);CHKERRQ(ierr); /*Here I get the PetscSection*/<br> for(PetscInt dof = 0; dof < 2; ++dof) {<br> row.c = dof;<br> for(PetscInt k = 0; k < 5; ++k) {<br>
col[k].c = dof;<br> }<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> 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> 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> ierr = MatSetValuesStencil(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);<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> Thanks,</div><div><br></div><div> Matt</div><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 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 class="h5"><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>