1 static char help[] = "Demonstrates the DMLocalToLocal bug in PETSc 3.6.\n\n"; 2 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 9 10 Contributed by Constantine Khroulev 11 */ 12 13 #include <petscdm.h> 14 #include <petscdmda.h> 15 16 PetscErrorCode PrintVecWithGhosts(DM da, Vec v) 17 { 18 PetscScalar **p; 19 PetscInt i, j; 20 PetscErrorCode ierr; 21 MPI_Comm com; 22 PetscMPIInt rank; 23 DMDALocalInfo info; 24 25 com = PetscObjectComm((PetscObject)da); 26 ierr = MPI_Comm_rank(com, &rank);CHKERRMPI(ierr); 27 ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr); 28 ierr = PetscSynchronizedPrintf(com, "begin rank %d portion (with ghosts, %D x %D)\n",rank, info.gxm, info.gym);CHKERRQ(ierr); 29 ierr = DMDAVecGetArray(da, v, &p);CHKERRQ(ierr); 30 for (i = info.gxs; i < info.gxs + info.gxm; i++) { 31 for (j = info.gys; j < info.gys + info.gym; j++) { 32 ierr = PetscSynchronizedPrintf(com, "%g, ", (double) PetscRealPart(p[j][i]));CHKERRQ(ierr); 33 } 34 ierr = PetscSynchronizedPrintf(com, "\n");CHKERRQ(ierr); 35 } 36 ierr = DMDAVecRestoreArray(da, v, &p);CHKERRQ(ierr); 37 ierr = PetscSynchronizedPrintf(com, "end rank %d portion\n", rank);CHKERRQ(ierr); 38 ierr = PetscSynchronizedFlush(com, PETSC_STDOUT);CHKERRQ(ierr); 39 return 0; 40 } 41 42 /* Set a Vec v to value, but do not touch ghosts. */ 43 PetscErrorCode VecSetOwned(DM da, Vec v, PetscScalar value) 44 { 45 PetscScalar **p; 46 PetscInt i, j, xs, xm, ys, ym; 47 PetscErrorCode ierr; 48 49 ierr = DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);CHKERRQ(ierr); 50 ierr = DMDAVecGetArray(da, v, &p);CHKERRQ(ierr); 51 for (i = xs; i < xs + xm; i++) { 52 for (j = ys; j < ys + ym; j++) { 53 p[j][i] = value; 54 } 55 } 56 ierr = DMDAVecRestoreArray(da, v, &p);CHKERRQ(ierr); 57 return 0; 58 } 59 60 int main(int argc, char **argv) 61 { 62 PetscInt M = 4, N = 3; 63 PetscErrorCode ierr; 64 DM da; 65 Vec local; 66 PetscScalar value = 0.0; 67 DMBoundaryType bx = DM_BOUNDARY_PERIODIC, by = DM_BOUNDARY_PERIODIC; 68 DMDAStencilType stype = DMDA_STENCIL_BOX; 69 70 ierr = PetscInitialize(&argc, &argv, (char*)0, help);if (ierr) return ierr; 71 ierr = DMDACreate2d(PETSC_COMM_WORLD,bx,by,stype,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr); 72 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 73 ierr = DMSetUp(da);CHKERRQ(ierr); 74 ierr = DMCreateLocalVector(da, &local);CHKERRQ(ierr); 75 76 ierr = VecSet(local, value);CHKERRQ(ierr); 77 ierr = PetscPrintf(PETSC_COMM_WORLD, "\nAfter setting all values to %d:\n", (int)PetscRealPart(value));CHKERRQ(ierr); 78 ierr = PrintVecWithGhosts(da, local);CHKERRQ(ierr); 79 ierr = PetscPrintf(PETSC_COMM_WORLD, "done\n");CHKERRQ(ierr); 80 81 value += 1.0; 82 /* set values owned by a process, leaving ghosts alone */ 83 ierr = VecSetOwned(da, local, value);CHKERRQ(ierr); 84 85 /* print after re-setting interior values again */ 86 ierr = PetscPrintf(PETSC_COMM_WORLD, "\nAfter setting interior values to %d:\n", (int)PetscRealPart(value));CHKERRQ(ierr); 87 ierr = PrintVecWithGhosts(da, local);CHKERRQ(ierr); 88 ierr = PetscPrintf(PETSC_COMM_WORLD, "done\n");CHKERRQ(ierr); 89 90 /* communicate ghosts */ 91 ierr = DMLocalToLocalBegin(da, local, INSERT_VALUES, local);CHKERRQ(ierr); 92 ierr = DMLocalToLocalEnd(da, local, INSERT_VALUES, local);CHKERRQ(ierr); 93 94 /* print again */ 95 ierr = PetscPrintf(PETSC_COMM_WORLD, "\nAfter local-to-local communication:\n");CHKERRQ(ierr); 96 ierr = PrintVecWithGhosts(da, local);CHKERRQ(ierr); 97 ierr = PetscPrintf(PETSC_COMM_WORLD, "done\n");CHKERRQ(ierr); 98 99 /* Free memory */ 100 ierr = VecDestroy(&local);CHKERRQ(ierr); 101 ierr = DMDestroy(&da);CHKERRQ(ierr); 102 ierr = PetscFinalize(); 103 return ierr; 104 } 105 106 /*TEST 107 108 test: 109 nsize: 2 110 111 TEST*/ 112