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