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; 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
