[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Bishesh Khanal
bisheshkh at gmail.com
Tue Oct 29 05:31:37 CDT 2013
On Mon, Oct 28, 2013 at 9:00 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Mon, Oct 28, 2013 at 2:51 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>
>>
>>
>>
>> On Mon, Oct 28, 2013 at 6:19 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>
>>> On Mon, Oct 28, 2013 at 10:13 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>
>>>>
>>>>
>>>>
>>>> On Mon, Oct 28, 2013 at 3:49 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>
>>>>> On Mon, Oct 28, 2013 at 9:06 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>>
>>>>>>
>>>>>> On Mon, Oct 28, 2013 at 1:40 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>>>
>>>>>>> On Mon, Oct 28, 2013 at 5:30 AM, Bishesh Khanal <bisheshkh at gmail.com
>>>>>>> > wrote:
>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On Fri, Oct 25, 2013 at 10:21 PM, Matthew Knepley <
>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>
>>>>>>>>> On Fri, Oct 25, 2013 at 2:55 PM, Bishesh Khanal <
>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>> On Fri, Oct 25, 2013 at 8:18 PM, Matthew Knepley <
>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>
>>>>>>>>>>> On Fri, Oct 25, 2013 at 12:09 PM, Bishesh Khanal <
>>>>>>>>>>> bisheshkh at gmail.com> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131029/44c74ce8/attachment-0001.html>
More information about the petsc-users
mailing list