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