[petsc-users] using petsc tools to solve isolated irregular domains with finite difference

Bishesh Khanal bisheshkh at gmail.com
Mon Oct 28 14:51:30 CDT 2013 

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


>    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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131028/ac85a9aa/attachment-0001.html>

More information about the petsc-users mailing list