1: static char help[] = "Demonstrates the DMLocalToLocal bug in PETSc 3.6.\n\n";

  3: /*
  4: Use the options
  5:      -da_grid_x <nx> - number of grid points in x direction, if M < 0
  6:      -da_grid_y <ny> - number of grid points in y direction, if N < 0
  7:      -da_processors_x <MX> number of processors in x directio
  8:      -da_processors_y <MY> number of processors in x direction

 10:   Contributed by Constantine Khroulev
 11: */

 13: #include <petscdm.h>
 14: #include <petscdmda.h>

 18: PetscErrorCode PrintVecWithGhosts(DM da, Vec v)
 19: {
 20:   PetscScalar    **p;
 21:   PetscInt       i, j;
 23:   MPI_Comm       com;
 24:   PetscMPIInt    rank;
 25:   DMDALocalInfo  info;

 27:   com = PetscObjectComm((PetscObject)da);
 28:   MPI_Comm_rank(com, &rank);

 30:   DMDAGetLocalInfo(da, &info);

 32:   PetscSynchronizedPrintf(com, "begin rank %d portion (with ghosts, %D x %D)\n",rank, info.gxm, info.gym);

 34:   DMDAVecGetArray(da, v, &p);
 35:   for (i = info.gxs; i < info.gxs + info.gxm; i++) {
 36:     for (j = info.gys; j < info.gys + info.gym; j++) {
 37:       PetscSynchronizedPrintf(com, "%g, ", (double) PetscRealPart(p[j][i]));
 38:     }
 39:     PetscSynchronizedPrintf(com, "\n");
 40:   }
 41:   DMDAVecRestoreArray(da, v, &p);

 43:   PetscSynchronizedPrintf(com, "end rank %d portion\n", rank);
 44:   PetscSynchronizedFlush(com, PETSC_STDOUT);
 45:   return 0;
 46: }

 48: /* Set a Vec v to value, but do not touch ghosts. */
 51: PetscErrorCode VecSetOwned(DM da, Vec v, PetscScalar value)
 52: {
 53:   PetscScalar    **p;
 54:   PetscInt         i, j, xs, xm, ys, ym;
 55:   PetscErrorCode   ierr;

 57:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
 58:   DMDAVecGetArray(da, v, &p);
 59:   for (i = xs; i < xs + xm; i++) {
 60:     for (j = ys; j < ys + ym; j++) {
 61:       p[j][i] = value;
 62:     }
 63:   }
 64:   DMDAVecRestoreArray(da, v, &p);

 66:   return 0;
 67: }

 71: int main(int argc, char **argv)
 72: {
 73:   PetscInt         M = -4, N = -3;
 74:   PetscErrorCode   ierr;
 75:   DM               da;
 76:   Vec              local;
 77:   PetscScalar      value = 0.0;
 78:   DMBoundaryType   bx    = DM_BOUNDARY_PERIODIC, by = DM_BOUNDARY_PERIODIC;
 79:   DMDAStencilType  stype = DMDA_STENCIL_BOX;

 81:   PetscInitialize(&argc, &argv, (char*)0, help);

 83:   /* Create distributed array and get vectors */
 84:   DMDACreate2d(PETSC_COMM_WORLD,
 85:                       bx, by,
 86:                       stype,
 87:                       M, N,
 88:                       PETSC_DECIDE, PETSC_DECIDE,
 89:                       1, 1,
 90:                       NULL, NULL,
 91:                       &da);

 93:   DMCreateLocalVector(da, &local);

 95:   VecSet(local, value);
 96:   PetscPrintf(PETSC_COMM_WORLD, "\nAfter setting all values to %d:\n", (int)PetscRealPart(value));
 97:   PrintVecWithGhosts(da, local);
 98:   PetscPrintf(PETSC_COMM_WORLD, "done\n");

100:   value += 1.0;
101:   /* set values owned by a process, leaving ghosts alone */
102:   VecSetOwned(da, local, value);

104:   /* print after re-setting interior values again */
105:   PetscPrintf(PETSC_COMM_WORLD, "\nAfter setting interior values to %d:\n", (int)PetscRealPart(value));
106:   PrintVecWithGhosts(da, local);
107:   PetscPrintf(PETSC_COMM_WORLD, "done\n");

109:   /* communicate ghosts */
110:   DMLocalToLocalBegin(da, local, INSERT_VALUES, local);
111:   DMLocalToLocalEnd(da, local, INSERT_VALUES, local);

113:   /* print again */
114:   PetscPrintf(PETSC_COMM_WORLD, "\nAfter local-to-local communication:\n");
115:   PrintVecWithGhosts(da, local);
116:   PetscPrintf(PETSC_COMM_WORLD, "done\n");

118:   /* Free memory */
119:   VecDestroy(&local);
120:   DMDestroy(&da);
121:   PetscFinalize();
122:   return 0;
123: }