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

Mon Oct 28 09:06:06 CDT 2013

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?

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