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

[Xiangdong]

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:

[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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170125/aaacc7d1/attachment.html>
-------------- next part --------------
diff --git a/ex45.c b/ex45.c
index b4bb565..6bcb074 100644
--- a/ex45.c
+++ b/ex45.c
@@ -117,8 +117,8 @@ PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void *ctx)
   DM             da;
   PetscErrorCode ierr;
   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
-  PetscScalar    v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
-  MatStencil     row,col[7];
+  PetscScalar    v[8],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
+  MatStencil     row,col[8];
 
   PetscFunctionBeginUser;
   ierr    = KSPGetDM(ksp,&da);CHKERRQ(ierr);
@@ -127,6 +127,9 @@ PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void *ctx)
   HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
   ierr    = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
 
+  MatSetOption(B,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
+  MatSetOption(jac,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
+
   for (k=zs; k<zs+zm; k++) {
     for (j=ys; j<ys+ym; j++) {
       for (i=xs; i<xs+xm; i++) {
@@ -142,7 +145,10 @@ PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void *ctx)
           v[4] = -HyHzdHx;col[4].i = i+1; col[4].j = j; col[4].k = k;
           v[5] = -HxHzdHy;col[5].i = i; col[5].j = j+1; col[5].k = k;
           v[6] = -HxHydHz;col[6].i = i; col[6].j = j; col[6].k = k+1;
-          ierr = MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
+          //ierr = MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
+
+	  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);
         }
       }
     }
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex45modify.c
Type: text/x-csrc
Size: 5473 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170125/aaacc7d1/attachment.c>
-------------- next part --------------
DM Object: 2 MPI processes
  type: da
Processor [0] M 7 N 7 P 7 m 1 n 1 p 2 w 1 s 1
X range of indices: 0 7, Y range of indices: 0 7, Z range of indices: 0 4
Processor [1] M 7 N 7 P 7 m 1 n 1 p 2 w 1 s 1
X range of indices: 0 7, Y range of indices: 0 7, Z range of indices: 4 7
[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: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.7.4, Oct, 02, 2016 
[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
[0]PETSC ERROR: #5 KSPSetUp() line 341 in petsc-3.7.4/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #6 KSPSolve() line 599 in petsc-3.7.4/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #7 main() line 51 in extest.c
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -dm_view
[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint at mcs.anl.gov----------