[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Matthew Knepley
knepley at gmail.com
Thu Oct 31 12:13:07 CDT 2013
On Thu, Oct 31, 2013 at 12:02 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
> On Tue, Oct 29, 2013 at 5:37 PM, Matthew Knepley <knepley at gmail.com>wrote:
>
>> On Tue, Oct 29, 2013 at 5:31 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>
>>> 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 ?
>>>
>>
>> Forget about multiple degrees of freedom until your scalar problem works.
>>
>>
>>> 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;
>>> }
>>>
>>
>> I can't imagine that this is what you want. You have two non-interacting
>> problems.
>>
>
> I just wanted a very basic example to use PetscSection with multiple dofs.
> In this example I'm solving Lap(u) = f where Lap is a component wise
> laplacian and u is a vector. So if I want to apply the Lap component-wise
> independently I could solve a system thrice for fx, fy and fz; but here
> just to demonstrate, I wanted to have all fx, fy and fz line up in a single
> vector.
>
> I did as you suggeted to write the values in the matrix, thanks. But now
> is there a way to let Petsc know automatically the size of the matrix it
> needs to create for the linear system based on the PetsSection I created?
> When I used only DMDA without Petscsection, this was possible, so I could
> just use:
> ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr);
> ierr =
> KSPSetComputeOperators(ksp,computeMatrix2d,(void*)&testPoisson);CHKERRQ(ierr);
> ierr =
> KSPSetComputeRHS(ksp,computeRhs2d,(void*)&testPoisson);CHKERRQ(ierr);
> where in my computeMatrix2d function I would set the non-zero values in
> the matrix.
>
DMCreateMatrix(), just as before.
Matt
> However, now since I do not want the rows in the matrix corresponding to
> the points in the grid which do not lie in the computational domain. This
> means the size of the matrix will not necessarily equal the total number of
> cells in DMDA. So in the following code:
>
> ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
> ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s);CHKERRQ(ierr);
> ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
> ierr = PetscSectionSetChart(s, 0, nC);CHKERRQ(ierr);
> for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {
> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {
> PetscInt point;
> if(isPosInDomain(&testPoisson,i,j,0)) {
> ierr = DMDAGetCellPoint(dm, i, j, 0, &point);CHKERRQ(ierr);
> ierr = PetscSectionSetDof(s, point, testPoisson.mDof); //
> No. of dofs associated with the point.
> }
>
> }
> }
> ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
> ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
> ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
>
> should I myself compute proper nC (i.e. total no. of points for which I
> want the matrix to have a corresponding row) before calling
> PetscSectionSetChart() or is the masking out of the points inside the loop
> with if(isPosInDomain(&testPoisson,i,j,0)) sufficient ?
> And, when you write:
> PetscSectionGetOffset(gs, p, &ind);
> row = ind < 0 ? -(ind+1) : ind;
>
> it seems you allow the possibility of getting a -ve ind when using
> PetscSectionGetOffset. When would an offset to a particular dof and point?
>
>
>
>
>>
>>> 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;
>>>
>>
>> Here you would call
>>
>> ierr = DMDAGetCellPoint(da, i, j, 0, &cell);CHKERRQ(ierr);
>> ierr = PetscSectionGetOffset(da, cell, &row);CHKERRQ(ierr);
>> ierr = PetscSectionGetDof(da, cell, &dof);CHKERRQ(ierr);
>>
>> This means that this cell has global rows = [row, row+dof).
>>
>>
>>> 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*/
>>> }
>>>
>>
>> You do the same thing for the columns, and then call
>>
>> ierr = MatSetValues(A, dof, rows, num*dof, cols, v,
>> INSERT_VALUES);CHKERRQ(ierr);
>>
>>
>>> 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
>>>>
>>>
>>>
>>
>>
>> --
>> 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/20131031/01b7d425/attachment-0001.html>
More information about the petsc-users
mailing list