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