[petsc-users] using petsc tools to solve isolated irregular domains with finite difference
Bishesh Khanal
bisheshkh at gmail.com
Thu Oct 31 17:02:11 CDT 2013
On Thu, Oct 31, 2013 at 7:43 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Thu, Oct 31, 2013 at 1:25 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>
>> I did not call DMCreateMatrix() before when using just dmda (and it is
>> working). In any case, I added a call to DMCreateMatrix() now but still
>> there are problems. Here are the code snippets and the error I get:
>>
>
> Then you were allowing it to be called automatically by PETSc, and it
> would have been this time as well.
>
>
>> //Setting up the section and linking it to DM:
>> 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);
>> ierr = DMCreateMatrix(dm,&A);CHKERRQ(ierr);
>>
>> //Set up KSP:
>> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
>> ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr);
>> ierr =
>> KSPSetComputeOperators(ksp,computeMatrix2dSection,(void*)&testPoisson);CHKERRQ(ierr);
>> ierr =
>> KSPSetComputeRHS(ksp,computeRhs2dSection,(void*)&testPoisson);CHKERRQ(ierr);
>> ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
>> ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
>>
>> ------------------------------------------------------
>> The computeMatrix2dSection function has:
>>
>> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
>> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
>> ierr = DMGetDefaultGlobalSection(da,&gs);CHKERRQ(ierr);
>> for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {
>> for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {
>> if (isPosInDomain(ctx,i,j)) {
>> ierr = DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);
>> ierr = PetscSectionGetOffset(gs,point,&row);
>> ierr = PetscSectionGetDof(gs,point,&dof);
>> for(PetscInt cDof = 0; cDof < dof; ++cDof) {
>> num = 0;
>> row+=cDof;
>> col[num] = row; //(i,j) position
>> v[num++] = -4;
>> if(isPosInDomain(ctx,i+1,j)) {
>> ierr =
>> DMDAGetCellPoint(da,i+1,j,0,&point);CHKERRQ(ierr);
>> ierr = PetscSectionGetOffset(gs,point,&col[num]);
>> col[num] += cDof;
>> v[num++] = 1;
>> }
>> if(isPosInDomain(ctx,i-1,j)) {
>> ierr =
>> DMDAGetCellPoint(da,i-1,j,0,&point);CHKERRQ(ierr);
>> ierr = PetscSectionGetOffset(gs,point,&col[num]);
>> col[num] += cDof;
>> v[num++] = 1;
>> }
>> if(isPosInDomain(ctx,i,j+1)) {
>> ierr =
>> DMDAGetCellPoint(da,i,j+1,0,&point);CHKERRQ(ierr);
>> ierr = PetscSectionGetOffset(gs,point,&col[num]);
>> col[num] += cDof;
>> v[num++] = 1;
>> }
>> if(isPosInDomain(ctx,i,j-1)) {
>> ierr =
>> DMDAGetCellPoint(da,i,j-1,0,&point);CHKERRQ(ierr);
>> ierr = PetscSectionGetOffset(gs,point,&col[num]);
>> col[num] += cDof;
>> v[num++] = 1;
>> }
>> ierr =
>> MatSetValues(A,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
>> }
>> }
>> }
>> }
>> ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>> ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>>
>> But this is not working. For e.g. for a 6X6 grid size I get the following
>> error:
>>
>
> Okay, here is how to debug this:
>
> 0) Go to a single scalar equations to make things easier
>
> 1) Use MatView(A, PETSC_VIEWER_STDOUT_WORLD), which should output a 36
> row matrix with
> the preallocated nonzero pattern, filled with zeros
>
> 2) Make sure this is the pattern you want
>
> 3) Run in the debugger with -start_in_debugger
>
> 4) When you get the error, see
>
> a) If the (i, j) is one that should be setting a value
>
> b) Why this (i, j) was not preallocated
>
Up to 4 (a), it is correct. There is a problem in the way DMDAGetCellPoint
and PetscSectionGetOffset works in this case scenario. I will try to
explain it below for the case of 4X4 grid.
First case:
If I set the computational domain to be all the points of the dmda grid,
(i.e. isPosInDomain(..,i,j,..) returns true for all i and j in the dmda
grid), the program runs fine and does not give any error.
Second case:
I want the computational domain to be some part of the whole grid. There is
a problem in this case.
The following test is in a case where,
isPosInDomain(..,i,j,..) returns true ONLY for (i,j) pairs of (2,1) (1,2)
(2,2) (3,2) and (3,2). The grid with its corresponding point number in the
petscsection is shown below:
12 13 14 15
8 9 10 11
4 5 6 7
0 1 2 3
of which 6, 9, 10, 11 and 14 correspond to the (i,j) points that returns
true for isPosInDomain(..,i,j,..)
MatView gives me the 16-row matrix with the star stencil non-zero structure
as expected.
The error I get is: new non-zero at (0,2) caused a malloc.
This error is for the value of (i,j) = (2,1), i.e. point 6.
The loop to set values in the matrix starts as:
for(PetscInt j = info.ys; j<info.ys+info.ym; ++j) {
for(PetscInt i = info.xs; i<info.xs+info.xm; ++i) {
if (isPosInDomain(ctx,i,j)) {
ierr = DMDAGetCellPoint(da,i,j,0,&point);CHKERRQ(ierr);
ierr = PetscSectionGetOffset(gs,point,&row);
Now, what I wanted from the beginning is to create the matrix containing
the rows corresponding to only those points (i,j) which has true values for
isPosInDomain(ctx,i,j), that means in this case 5 row matrix. In the
current test, the first (i,j) that reaches DMDAGetCellPoint is (2,1), which
gives sets point variable to 6.
The line PetscSectionGetOffset(gs,point,&row) sets the value of row to 0.
So I think this where the inconsistency lies.
MatSetValues(A,1,&row,num,col,v,INSERT_VALUES) expects row variable and
col[] array to correspond to (i,j) based on DMDA, but PetsSectionGetOffset
gives result based on how we masked away some of the points from dmda grid.
Or am I missing something very basic here ?
> Thanks,
>
> Matt
>
>
>> [0]PETSC ERROR: --------------------- Error Message
>> ------------------------------------
>> [0]PETSC ERROR: Argument out of range!
>> [0]PETSC ERROR: New nonzero at (0,6) caused a malloc!
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: Petsc Development GIT revision:
>> 61e5e40bb2c5bf2423e94b71a15fef47e411b0da GIT Date: 2013-10-25 21:47:45
>> -0500
>> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
>> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
>> [0]PETSC ERROR: See docs/index.html for manual pages.
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: build/poissonIrregular on a arch-linux2-cxx-debug named
>> edwards by bkhanal Thu Oct 31 19:23:58 2013
>> [0]PETSC ERROR: Libraries linked from
>> /home/bkhanal/Documents/softwares/petsc/arch-linux2-cxx-debug/lib
>> [0]PETSC ERROR: Configure run at Sat Oct 26 16:35:15 2013
>> [0]PETSC ERROR: Configure options --download-mpich
>> -download-f-blas-lapack=1 --download-metis --download-parmetis
>> --download-superlu_dist --download-scalapack --download-mumps
>> --download-hypre --with-clanguage=cxx
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: MatSetValues_SeqAIJ() line 413 in
>> src/mat/impls/aij/seq/aij.c
>> [0]PETSC ERROR: MatSetValues() line 1130 in src/mat/interface/matrix.c
>> [0]PETSC ERROR: computeMatrix2dSection() line 318 in
>> /user/bkhanal/home/works/poissonIrregular/poissonIrregular.cxx
>> [0]PETSC ERROR: KSPSetUp() line 228 in src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: KSPSolve() line 399 in src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: main() line 598 in
>> /user/bkhanal/home/works/poissonIrregular/poissonIrregular.cxx
>> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0
>> [unset]: aborting job:
>> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0
>>
>>
>>> 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
>>>
>>
>>
>
>
> --
> 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/5164a25c/attachment-0001.html>
More information about the petsc-users
mailing list