On Thu, Oct 31, 2013 at 12:02 PM, Bishesh Khanal <[email protected]>wrote:
> On Tue, Oct 29, 2013 at 5:37 PM, Matthew Knepley <[email protected]>wrote: > >> On Tue, Oct 29, 2013 at 5:31 AM, Bishesh Khanal <[email protected]>wrote: >> >>> On Mon, Oct 28, 2013 at 9:00 PM, Matthew Knepley <[email protected]>wrote: >>> >>>> On Mon, Oct 28, 2013 at 2:51 PM, Bishesh Khanal <[email protected]>wrote: >>>> >>>>> On Mon, Oct 28, 2013 at 6:19 PM, Matthew Knepley <[email protected]>wrote: >>>>> >>>>>> On Mon, Oct 28, 2013 at 10:13 AM, Bishesh Khanal <[email protected] >>>>>> > wrote: >>>>>> >>>>>>> On Mon, Oct 28, 2013 at 3:49 PM, Matthew Knepley >>>>>>> <[email protected]>wrote: >>>>>>> >>>>>>>> On Mon, Oct 28, 2013 at 9:06 AM, Bishesh Khanal < >>>>>>>> [email protected]> wrote: >>>>>>>> >>>>>>>>> On Mon, Oct 28, 2013 at 1:40 PM, Matthew Knepley < >>>>>>>>> [email protected]> wrote: >>>>>>>>> >>>>>>>>>> On Mon, Oct 28, 2013 at 5:30 AM, Bishesh Khanal < >>>>>>>>>> [email protected]> wrote: >>>>>>>>>> >>>>>>>>>>> On Fri, Oct 25, 2013 at 10:21 PM, Matthew Knepley < >>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>> >>>>>>>>>>>> On Fri, Oct 25, 2013 at 2:55 PM, Bishesh Khanal < >>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>> >>>>>>>>>>>>> On Fri, Oct 25, 2013 at 8:18 PM, Matthew Knepley < >>>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>>> >>>>>>>>>>>>>> On Fri, Oct 25, 2013 at 12:09 PM, Bishesh Khanal < >>>>>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>>>>> >>>>>>>>>>>>>>> Dear all, >>>>>>>>>>>>>>> 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): >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> Let's say that I have a 3D binary image (a cube). >>>>>>>>>>>>>>> The binary information of the image partitions the cube into >>>>>>>>>>>>>>> a computational domain and non-computational domain. >>>>>>>>>>>>>>> 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 >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> 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). >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> 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. >>>>>>>>>>>>>>> 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. >>>>>>>>>>>>>>> Then, once the system is solved, put the values of the >>>>>>>>>>>>>>> solution back to the corresponding positions in the binary >>>>>>>>>>>>>>> image. >>>>>>>>>>>>>>> Later, I might have to extend this for the staggered grid >>>>>>>>>>>>>>> case too. >>>>>>>>>>>>>>> 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 ? >>>>>>>>>>>>>>> >>>>>>>>>>>>>>> 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!) >>>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> 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 >>>>>>>>>>>>>> alpha code. If you feel comfortable with that level of >>>>>>>>>>>>>> development, we can help you. >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> 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. >>>>>>>>>>>>> Just a note: I'm currently using the petsc3.4.3 and have not >>>>>>>>>>>>> used the development version before. >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> Okay, >>>>>>>>>>>> >>>>>>>>>>>> 1) clone the repository using Git and build the 'next' branch. >>>>>>>>>>>> >>>>>>>>>>>> 2) then we will need to create a PetscSection that puts >>>>>>>>>>>> unknowns where you want them >>>>>>>>>>>> >>>>>>>>>>>> 3) Setup the solver as usual >>>>>>>>>>>> >>>>>>>>>>>> You can do 1) an 3) before we do 2). >>>>>>>>>>>> >>>>>>>>>>>> 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). >>>>>>>>>>> Please let me know what I should do now. >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> Okay, now write a loop to setup the PetscSection. I have the DMDA >>>>>>>>>> partitioning cells, so you would have >>>>>>>>>> unknowns in cells. The code should look like this: >>>>>>>>>> >>>>>>>>>> PetscSectionCreate(comm, &s); >>>>>>>>>> DMDAGetNumCells(dm, NULL, NULL, NULL, &nC); >>>>>>>>>> PetscSectionSetChart(s, 0, nC); >>>>>>>>>> for (k = zs; k < zs+zm; ++k) { >>>>>>>>>> for (j = ys; j < ys+ym; ++j) { >>>>>>>>>> for (i = xs; i < xs+xm; ++i) { >>>>>>>>>> PetscInt point; >>>>>>>>>> >>>>>>>>>> DMDAGetCellPoint(dm, i, j, k, &point); >>>>>>>>>> PetscSectionSetDof(s, point, dof); // You know how many dof >>>>>>>>>> are on each vertex >>>>>>>>>> } >>>>>>>>>> } >>>>>>>>>> } >>>>>>>>>> PetscSectionSetUp(s); >>>>>>>>>> DMSetDefaultSection(dm, s); >>>>>>>>>> PetscSectionDestroy(&s); >>>>>>>>>> >>>>>>>>>> I will merge the necessary stuff into 'next' to make this work. >>>>>>>>>> >>>>>>>>> >>>>>>>>> I have put an example without the PetscSection here: >>>>>>>>> https://github.com/bishesh/petscPoissonIrregular/blob/master/poissonIrregular.cxx >>>>>>>>> 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? >>>>>>>>> >>>>>>>> >>>>>>>> 1) Give me permission to comment on the source (I am 'knepley') >>>>>>>> >>>>>>>> 2) You mask out the (i,j,k) that you do not want in that loop >>>>>>>> >>>>>>> >>>>>>> Done. >>>>>>> I mask it out using isPosInDomain() function:: >>>>>>> if(isPosInDomain(&testPoisson,i,j,k)) { >>>>>>> ierr = DMDAGetCellPoint(dm, i, j, k, &point);CHKERRQ( >>>>>>> ierr); >>>>>>> ierr = PetscSectionSetDof(s, point, testPoisson.mDof); // >>>>>>> You know how many dof are on each vertex >>>>>>> } >>>>>>> >>>>>>> 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. >>>>>>> >>>>>> >>>>>> Pushed. >>>>>> >>>>> >>>>> Pulled it thanks, but compiler was complaining that it didn't find >>>>> DMDAGetCellPoint. >>>>> Doing grep -R 'DMDAGetCellPoint' include/ >>>>> returned nothing although the implementation is there in dalocal.c >>>>> So I added the declaration to petscdmda.h file just above the >>>>> declaration of DMDAGetNumCells. >>>>> >>>> >>>> I will add the declaration. >>>> >>>> >>>>> 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 ? >>>>> >>>> >>>> Instead of MatSetValuesStencil(), you would go back to using >>>> MatSetValues(), and calculate the indices using PetscSection. >>>> You get the section >>>> >>>> DMGetDefaultGlobalSection(dm, &gs) >>>> >>>> and then >>>> >>>> DMDAGetCellPoint(dm, i, j, k, &p); >>>> PetscSectionGetOffset(gs, p, &ind); >>>> row = ind < 0 ? -(ind+1) : ind; >>>> >>>> >>> >>> 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 ? >>> How would you change the following code that sets 2D laplacian stencil >>> for 2 dof to use petscsection and Matsetvalues instead ? >>> >> >> Forget about multiple degrees of freedom until your scalar problem works. >> >> >>> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr); >>> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); >>> ierr = DMGetDefaultGlobalSection(dm,&gs);CHKERRQ(ierr); /*Here I get >>> the PetscSection*/ >>> for(PetscInt dof = 0; dof < 2; ++dof) { >>> row.c = dof; >>> for(PetscInt k = 0; k < 5; ++k) { >>> col[k].c = dof; >>> } >>> >> >> I can't imagine that this is what you want. You have two non-interacting >> problems. >> > > 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. > > 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: > ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); > ierr = > KSPSetComputeOperators(ksp,computeMatrix2d,(void*)&testPoisson);CHKERRQ(ierr); > ierr = > KSPSetComputeRHS(ksp,computeRhs2d,(void*)&testPoisson);CHKERRQ(ierr); > where in my computeMatrix2d function I would set the non-zero values in > the matrix. > DMCreateMatrix(), just as before. Matt > 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: > > ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); > ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s);CHKERRQ(ierr); > ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); > ierr = PetscSectionSetChart(s, 0, nC);CHKERRQ(ierr); > for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) { > for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) { > PetscInt point; > if(isPosInDomain(&testPoisson,i,j,0)) { > ierr = DMDAGetCellPoint(dm, i, j, 0, &point);CHKERRQ(ierr); > ierr = PetscSectionSetDof(s, point, testPoisson.mDof); // > No. of dofs associated with the point. > } > > } > } > ierr = PetscSectionSetUp(s);CHKERRQ(ierr); > ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr); > ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); > > 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 ? > And, when you write: > PetscSectionGetOffset(gs, p, &ind); > row = ind < 0 ? -(ind+1) : ind; > > it seems you allow the possibility of getting a -ve ind when using > PetscSectionGetOffset. When would an offset to a particular dof and point? > > > > >> >>> for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) { >>> for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) { >>> num = 0; >>> /*ierr = >>> DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);*/ /*Get the PetscPoint*/ >>> /*But I did now understand how I would emulate the row.c >>> and col.c info with PetscSection*/ >>> row.i = i; row.j = j; >>> >> >> Here you would call >> >> ierr = DMDAGetCellPoint(da, i, j, 0, &cell);CHKERRQ(ierr); >> ierr = PetscSectionGetOffset(da, cell, &row);CHKERRQ(ierr); >> ierr = PetscSectionGetDof(da, cell, &dof);CHKERRQ(ierr); >> >> This means that this cell has global rows = [row, row+dof). >> >> >>> col[num].i = i; col[num].j = j; >>> if (isPosInDomain(ctx,i,j)) { /*put the operator >>> for only the values for only the points lying in the domain*/ >>> v[num++] = -4; >>> if(isPosInDomain(ctx,i+1,j)) { >>> col[num].i = i+1; col[num].j = j; >>> v[num++] = 1; >>> } >>> if(isPosInDomain(ctx,i-1,j)) { >>> col[num].i = i-1; col[num].j = j; >>> v[num++] = 1; >>> } >>> if(isPosInDomain(ctx,i,j+1)) { >>> col[num].i = i; col[num].j = j+1; >>> v[num++] = 1; >>> } >>> if(isPosInDomain(ctx,i,j-1)) { >>> col[num].i = i; col[num].j = j-1; >>> v[num++] = 1; >>> } >>> } else { >>> v[num++] = dScale; /*set Dirichlet identity rows >>> for points in the rectangle but outside the computational domain*/ >>> } >>> >> >> You do the same thing for the columns, and then call >> >> ierr = MatSetValues(A, dof, rows, num*dof, cols, v, >> INSERT_VALUES);CHKERRQ(ierr); >> >> >>> ierr = >>> MatSetValuesStencil(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr); >>> } >>> } >>> } >>> >>> >>>> Thanks, >>>> >>>> Matt >>>> >>>> >>>>> >>>>>> Matt >>>>>> >>>>>> >>>>>>> >>>>>>>> Matt >>>>>>>> >>>>>>>> >>>>>>>>> >>>>>>>>>> Thanks, >>>>>>>>>> >>>>>>>>>> Matt >>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>>> If not, just put the identity into >>>>>>>>>>>>>> the rows you do not use on the full cube. It will not hurt >>>>>>>>>>>>>> scalability or convergence. >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> 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: >>>>>>>>>>>>> 1) Having read Jed's answer here : >>>>>>>>>>>>> http://scicomp.stackexchange.com/questions/3426/why-is-pinning-a-point-to-remove-a-null-space-bad/3427#3427 >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> Jed is talking about a constraint on a the pressure at a point. >>>>>>>>>>>> This is just decoupling these unknowns from the rest >>>>>>>>>>>> of the problem. >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>>> 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. >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> It can change the exact iteration, but it does not make the >>>>>>>>>>>> matrix conditioning worse. >>>>>>>>>>>> >>>>>>>>>>>> Matt >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>>> 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. >>>>>>>>>>>>> >>>>>>>>>>>>> Thanks, >>>>>>>>>>>>> Bishesh >>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> Matt >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>>> Thanks, >>>>>>>>>>>>>>> Bishesh >>>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> >>>>>>>>>>>>>> -- >>>>>>>>>>>>>> What most experimenters take for granted before they begin >>>>>>>>>>>>>> their experiments is infinitely more interesting than any >>>>>>>>>>>>>> results to which >>>>>>>>>>>>>> their experiments lead. >>>>>>>>>>>>>> -- Norbert Wiener >>>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> -- >>>>>>>>>>>> What most experimenters take for granted before they begin >>>>>>>>>>>> their experiments is infinitely more interesting than any results >>>>>>>>>>>> to which >>>>>>>>>>>> their experiments lead. >>>>>>>>>>>> -- Norbert Wiener >>>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> -- >>>>>>>>>> What most experimenters take for granted before they begin their >>>>>>>>>> experiments is infinitely more interesting than any results to which >>>>>>>>>> their >>>>>>>>>> experiments lead. >>>>>>>>>> -- Norbert Wiener >>>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> -- >>>>>>>> What most experimenters take for granted before they begin their >>>>>>>> experiments is infinitely more interesting than any results to which >>>>>>>> their >>>>>>>> experiments lead. >>>>>>>> -- Norbert Wiener >>>>>>>> >>>>>>> >>>>>>> >>>>>> >>>>>> >>>>>> -- >>>>>> What most experimenters take for granted before they begin their >>>>>> experiments is infinitely more interesting than any results to which >>>>>> their >>>>>> experiments lead. >>>>>> -- Norbert Wiener >>>>>> >>>>> >>>>> >>>> >>>> >>>> -- >>>> What most experimenters take for granted before they begin their >>>> experiments is infinitely more interesting than any results to which their >>>> experiments lead. >>>> -- Norbert Wiener >>>> >>> >>> >> >> >> -- >> What most experimenters take for granted before they begin their >> experiments is infinitely more interesting than any results to which their >> experiments lead. >> -- Norbert Wiener >> > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
