[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Matthew Knepley
knepley at gmail.com
Mon Oct 28 15:00:17 CDT 2013
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;
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/20131028/1e1d0584/attachment.html>
More information about the petsc-users
mailing list