<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Thu, Oct 31, 2013 at 6:13 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 12: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 Tue, Oct 29, 2013 at 5:37 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 Tue, Oct 29, 2013 at 5:31 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"><div><div>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>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:</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 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:</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 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:</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"><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:</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">
<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><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></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></div></div></div></div></blockquote><div><br></div></div></div><div>Forget about multiple degrees of freedom until your scalar problem works.</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 = 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></div></div></div></div></blockquote><div><br></div></div><div>I can't imagine that this is what you want. You have two non-interacting problems.</div></div></div></div></blockquote>
<div><br></div><div>I just wanted a very basic example to use PetscSection with multiple dofs. In this example I'm solving Lap(u) = f where Lap is a component wise laplacian and u is a vector. So if I want to apply the Lap component-wise independently I could solve a system thrice for fx, fy and fz; but here just to demonstrate, I wanted to have all fx, fy and fz line up in a single vector.<br>
</div><div></div><div><br>I did as you suggeted to write the values in the matrix, thanks. But now is there a way to let Petsc know automatically the size of the matrix it needs to create for the linear system based on the PetsSection I created? When I used only DMDA without Petscsection, this was possible, so I could just use:<br>
ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr);<br> ierr = KSPSetComputeOperators(ksp,computeMatrix2d,(void*)&testPoisson);CHKERRQ(ierr);<br> ierr = KSPSetComputeRHS(ksp,computeRhs2d,(void*)&testPoisson);CHKERRQ(ierr);<br>
</div><div>where in my computeMatrix2d function I would set the non-zero values in the matrix. <br></div></div></div></div></blockquote><div><br></div></div></div><div>DMCreateMatrix(), just as before.</div></div></div></div>
</blockquote><div><br></div><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>
<br></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>
<br>[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"><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>