In the natural numbering on the mesh it is I = ix + iy*Nx + iz*Ny*Nx but PETSc does not use the natural numbering. You can use
DMDAGetAO() then AOApplicationToPetsc() to make the value from the natural numbering to the PETSc number, then this number can be used with MatSetValues(), VecSetValues(), What Matt describes below is a more cumbersome way to achieve the same result. Barry > On Jan 26, 2017, at 10:29 AM, Matthew Knepley <[email protected]> wrote: > > On Thu, Jan 26, 2017 at 10:20 AM, Xiangdong <[email protected]> wrote: > Thanks, Matt. For a cell in DMDA 3d with global id (ix,iy,iz), what is the > global row id of that cell corresponding to the matrix generated by > DMCreateMatrix? It is not always ix + iy*Nx + iz*Ny*Nx for different > da_processors_xyz, is it? How can I obtain that global row id? > > It is not simple: > > 1) Figure out the process division using > > > http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMDA/DMDAGetInfo.html > > 2) Figure out which process p your i,j,k in on > > You have to march through the processes, knowing how big the chunk of > indices is on each, > which you get from 1) and > > > http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMDA/DMDAGetOwnershipRanges.html > > 3) Add up the sizes from all process less than p, and then the offset on > process p. This the global offset > > This is not a scalable thing to do, which is why we do not include it in the > API. > > Matt > > Thanks. > > Xiangdong > > > > On Wed, Jan 25, 2017 at 10:53 AM, Matthew Knepley <[email protected]> wrote: > On Wed, Jan 25, 2017 at 9:17 AM, Xiangdong <[email protected]> wrote: > Hello everyone, > > I have a question on setting matrix entries which are not in the stencil > width. Take ksp ex45.c as an example, > http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex45.c.html > > Instead of setting up the standard 7-point stencil, now for each cell, the > matrix also has a additional dependency on the cell (Mx, My, Mz). Namely, for > each row, the col corresponding to (Mx, My, Mz) is always nonzero. I modify > the example code to add this entries: > > + MatSetOption(B,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE); > + MatSetOption(jac,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE); > > + v[7] = 100; col[7].i = mx-1; col[7].j = my-1; col[7].k = mz-1; > + ierr = > MatSetValuesStencil(B,1,&row,8,col,v,INSERT_VALUES);CHKERRQ(ierr); > > It is okay to for np=1, but crash for np>=2 with the error message: > > You can do this, but > > 1) You cannot use MatSetStencil, since your entry is not actually in your > stencil. You will have to make a separate call to MatSetValues() using the > global index. > > 2) The nonzero pattern we have allocated will be wrong. You will have to set > the MatOption which gives an error on new nonzeros to PETSC_FALSE. > > 3) You will have a dense row in your Jacobian, which is hard to make perform > well, and also stymies most preconditioners. > > Thanks, > > Matt > > [0]PETSC ERROR: --------------------- Error Message > -------------------------------------------------------------- > [0]PETSC ERROR: Argument out of range > [0]PETSC ERROR: Local index 342 too large 244 (max) at 7 > > [0]PETSC ERROR: #1 ISLocalToGlobalMappingApply() line 423 in > petsc-3.7.4/src/vec/is/utils/isltog.c > [0]PETSC ERROR: #2 MatSetValuesLocal() line 2052 in > petsc-3.7.4/src/mat/interface/matrix.c > [0]PETSC ERROR: #3 MatSetValuesStencil() line 1447 in > petsc-3.7.4/src/mat/interface/matrix.c > [0]PETSC ERROR: #4 ComputeMatrix() line 151 in extest.c > > Can I add new entries to the cols not in the stencil width into the dmda > matrix or Jacobian? > > Attached please find the modifed ex45 example, the diff file as well as the > run log. > > Thanks for your help. > > Best, > Xiangdong > > > > > -- > 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
