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