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 ? 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; } 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; 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*/ } 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 >
