[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Bishesh Khanal
bisheshkh at gmail.com
Thu Oct 31 12:02:16 CDT 2013
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.
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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131031/a99ad675/attachment-0001.html>
More information about the petsc-users
mailing list