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

Matthew Knepley knepley at gmail.com
Mon Nov 4 09:54:23 CST 2013


On Mon, Nov 4, 2013 at 4:49 AM, Bishesh Khanal <bisheshkh at gmail.com> wrote:

> On Fri, Nov 1, 2013 at 1:30 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Fri, Nov 1, 2013 at 12:08 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>
>>> On Thu, Oct 31, 2013 at 11:34 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>
>>>> On Thu, Oct 31, 2013 at 5:02 PM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>
>>>>> 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 ?
>>>>>
>>>>
>>>> You are missing something basic. The row is correctly identified as 0,
>>>> and that means the 2 is
>>>> point 10.
>>>>
>>> The problem must be preallocation. First, you should send the result of
>>>> MatView().
>>>>
>>>
>>> MatView shows that the preallocation and the matrix created is for the
>>> full dmda grid which is probably not what we want here, right ? Here is the
>>> matview result:
>>>
>>
>> Thats right. Did you call
>>
>>   DMSetDefaultSection(da, s)
>>
>> after you made the new section and before you called DMGetMatrix()?
>>
> Do you mean before I called DMCreateMatrix() ? I called
> DMSetDefaultSection(da,s) as below:
>
> ...
>  ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
>     ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
>     ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
> //    ierr = DMGetMatrix(dm,&A);CHKERRQ(ierr);  //It seems DMGetMatrix is
> changed to DMCreateMatrix().
>     ierr = DMCreateMatrix(dm,&A);CHKERRQ(ierr);
>     ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
>

Yes, this is correct. I do exactly this in the example code and it works
fine. I suspect
you are not creating the section you think you are. You can send the code
if you want
and I will look at it.

   Matt


> This MatView on A shows (result below) that it still has 16 rows for all
> DMDA points instead of only those petscsection points for which I set a
> non-zero dof. Is it what it should not be doing ?
>
>>  Mat Object: 1 MPI processes
>>>
>>> type: seqaij
>>>
>>> row 0: (0, 0) (1, 0) (4, 0)
>>>
>>> row 1: (0, 0) (1, 0) (2, 0) (5, 0)
>>>
>>> row 2: (1, 0) (2, 0) (3, 0) (6, 0)
>>>
>>> row 3: (2, 0) (3, 0) (7, 0)
>>>
>>> row 4: (0, 0) (4, 0) (5, 0) (8, 0)
>>>
>>> row 5: (1, 0) (4, 0) (5, 0) (6, 0) (9, 0)
>>>
>>> row 6: (2, 0) (5, 0) (6, 0) (7, 0) (10, 0)
>>>
>>> row 7: (3, 0) (6, 0) (7, 0) (11, 0)
>>>
>>> row 8: (4, 0) (8, 0) (9, 0) (12, 0)
>>>
>>> row 9: (5, 0) (8, 0) (9, 0) (10, 0) (13, 0)
>>>
>>> row 10: (6, 0) (9, 0) (10, 0) (11, 0) (14, 0)
>>>
>>> row 11: (7, 0) (10, 0) (11, 0) (15, 0)
>>>
>>> row 12: (8, 0) (12, 0) (13, 0)
>>>
>>> row 13: (9, 0) (12, 0) (13, 0) (14, 0)
>>>
>>> row 14: (10, 0) (13, 0) (14, 0) (15, 0)
>>>
>>> row 15: (11, 0) (14, 0) (15, 0)
>>>
>>> [0]PETSC ERROR: --------------------- Error Message
>>> ------------------------------------
>>>
>>> [0]PETSC ERROR: Argument out of range!
>>>
>>> [0]PETSC ERROR: New nonzero at (0,2) caused a malloc!
>>>
>>>
>>>
> For reference here is the code snippet again:
> ...
> 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);
>     ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
>
> 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);
>
> ....
> The function computeMatrix2dSection:
>
> #undef __FUNCT__
> #define __FUNCT__ "computeMatrix2dSection"
> PetscErrorCode computeMatrix2dSection(KSP ksp, Mat A, Mat B, MatStructure
> *str, void *context) {
>     PoissonModel            *ctx = (PoissonModel*)context;
>     PetscErrorCode          ierr;
>     DM                      da;
>     DMDALocalInfo           info;
>     PetscInt                row, col[5];
>     PetscInt                dof;
>     PetscScalar             v[5]; //array to store 5 point stencil for
> each row
>     PetscInt                num; //non-zero position in the current row
>     //    PetscScalar             dScale = 1;     //to scale Dirichlet
> identity rows if needed
>     PetscSection            gs;             //Global Section that keeps
> the grid info and indices
>     PetscInt                point;          //Current point of the
> petscSection
>
>     PetscFunctionBeginUser;
>     ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
>     ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
>     ierr = DMGetDefaultGlobalSection(da,&gs);CHKERRQ(ierr);
>     ierr = DMCreateMatrix(da,&A);CHKERRQ(ierr);
> //    ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);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 < ctx->mDof; ++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);
>
>
>
>
>>
>>>>    Matt
>>>>
>>>>
>>>>>
>>>>>>   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
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>> --
>>>> 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/20131104/3102f5f9/attachment-0001.html>


More information about the petsc-users mailing list