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 ? > 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 >
