On Thu, Oct 31, 2013 at 5:02 PM, Bishesh Khanal <[email protected]> wrote:
> > > > On Thu, Oct 31, 2013 at 7:43 PM, Matthew Knepley <[email protected]>wrote: > >> On Thu, Oct 31, 2013 at 1:25 PM, Bishesh Khanal <[email protected]>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(). 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
