[petsc-users] MatSetValuesStencil for cols not in the dmda stencil width

[Barry Smith]

   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 <knepley at gmail.com> wrote:
> 
> On Thu, Jan 26, 2017 at 10:20 AM, Xiangdong <epscodes at gmail.com> 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 <knepley at gmail.com> wrote:
> On Wed, Jan 25, 2017 at 9:17 AM, Xiangdong <epscodes at gmail.com> 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