Actual source code: ex6.c


  2: static char help[] = "\n\n";

  4: /*
  5:      Demonstrates using DM_BOUNDARY_GHOSTED how to handle a rotated boundary conditions where one edge
  6:     is connected to its immediate neighbor

  8:     Consider the domain (with natural numbering)

 10:      6   7   8
 11:      3   4   5
 12:      0   1   2

 14:     The ghost points along the bottom (directly below the three columns above) should be 0 3 and 6
 15:     while the ghost points along the left side should be 0 1 2

 17:     Note that the ghosted local vectors extend in both the x and y directions so, for example if we have a
 18:     single MPI process the ghosted vector has (in the original natural numbering)

 20:      x  x  x  x  x
 21:      2  6  7  8  x
 22:      1  3  4  5  x
 23:      0  0  1  2  x
 24:      x  0  3  6  x

 26:     where x indicates a location that is not updated by the communication and should be used.

 28:     For this to make sense the number of grid points in the x and y directions must be the same

 30:     This ghost point mapping was suggested by: Wenbo Zhao <zhaowenbo.npic@gmail.com>
 31: */

 33: #include <petscdm.h>
 34: #include <petscdmda.h>

 36: int main(int argc, char **argv)
 37: {
 38:   PetscInt     M = 6;
 39:   DM           da;
 40:   Vec          local, global, natural;
 41:   PetscInt     i, start, end, *ifrom, x, y, xm, ym;
 42:   PetscScalar *xnatural;
 43:   IS           from, to;
 44:   AO           ao;
 45:   VecScatter   scatter1, scatter2;
 46:   PetscViewer  subviewer;

 48:   PetscFunctionBeginUser;
 49:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 51:   /* Create distributed array and get vectors */
 52:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DMDA_STENCIL_STAR, M, M, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
 53:   PetscCall(DMSetUp(da));
 54:   PetscCall(DMCreateGlobalVector(da, &global));
 55:   PetscCall(DMCreateLocalVector(da, &local));

 57:   /* construct global to local scatter for the left side of the domain to the ghost on the bottom */
 58:   PetscCall(DMDAGetCorners(da, &x, &y, NULL, &xm, &ym, NULL));
 59:   if (!y) { /* only processes on the bottom of the domain fill up the ghost locations */
 60:     PetscCall(ISCreateStride(PETSC_COMM_SELF, xm, 1, 1, &to));
 61:   } else {
 62:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 0, &to));
 63:   }
 64:   PetscCall(PetscMalloc1(xm, &ifrom));
 65:   for (i = x; i < x + xm; i++) ifrom[i - x] = M * i;
 66:   PetscCall(DMDAGetAO(da, &ao));
 67:   PetscCall(AOApplicationToPetsc(ao, xm, ifrom));
 68:   if (!y) {
 69:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, xm, ifrom, PETSC_OWN_POINTER, &from));
 70:   } else {
 71:     PetscCall(PetscFree(ifrom));
 72:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_COPY_VALUES, &from));
 73:   }
 74:   PetscCall(VecScatterCreate(global, from, local, to, &scatter1));
 75:   PetscCall(ISDestroy(&to));
 76:   PetscCall(ISDestroy(&from));

 78:   /* construct global to local scatter for the bottom side of the domain to the ghost on the right */
 79:   if (!x) { /* only processes on the left side of the domain fill up the ghost locations */
 80:     PetscCall(ISCreateStride(PETSC_COMM_SELF, ym, xm + 2, xm + 2, &to));
 81:   } else {
 82:     PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 0, &to));
 83:   }
 84:   PetscCall(PetscMalloc1(ym, &ifrom));
 85:   for (i = y; i < y + ym; i++) ifrom[i - y] = i;
 86:   PetscCall(DMDAGetAO(da, &ao));
 87:   PetscCall(AOApplicationToPetsc(ao, ym, ifrom));
 88:   if (!x) {
 89:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ym, ifrom, PETSC_OWN_POINTER, &from));
 90:   } else {
 91:     PetscCall(PetscFree(ifrom));
 92:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_COPY_VALUES, &from));
 93:   }
 94:   PetscCall(VecScatterCreate(global, from, local, to, &scatter2));
 95:   PetscCall(ISDestroy(&to));
 96:   PetscCall(ISDestroy(&from));

 98:   /*
 99:      fill the global vector with the natural global numbering for each local entry
100:      this is only done for testing purposes since it is easy to see if the scatter worked correctly
101:   */
102:   PetscCall(DMDACreateNaturalVector(da, &natural));
103:   PetscCall(VecGetOwnershipRange(natural, &start, &end));
104:   PetscCall(VecGetArray(natural, &xnatural));
105:   for (i = start; i < end; i++) xnatural[i - start] = i;
106:   PetscCall(VecRestoreArray(natural, &xnatural));
107:   PetscCall(DMDANaturalToGlobalBegin(da, natural, INSERT_VALUES, global));
108:   PetscCall(DMDANaturalToGlobalEnd(da, natural, INSERT_VALUES, global));
109:   PetscCall(VecDestroy(&natural));

111:   /* scatter from global to local */
112:   PetscCall(VecScatterBegin(scatter1, global, local, INSERT_VALUES, SCATTER_FORWARD));
113:   PetscCall(VecScatterEnd(scatter1, global, local, INSERT_VALUES, SCATTER_FORWARD));
114:   PetscCall(VecScatterBegin(scatter2, global, local, INSERT_VALUES, SCATTER_FORWARD));
115:   PetscCall(VecScatterEnd(scatter2, global, local, INSERT_VALUES, SCATTER_FORWARD));
116:   /*
117:      normally here you would also call
118:   PetscCall(DMGlobalToLocalBegin(da,global,INSERT_VALUES,local));
119:   PetscCall(DMGlobalToLocalEnd(da,global,INSERT_VALUES,local));
120:     to update all the interior ghost cells between neighboring processes.
121:     We don't do it here since this is only a test of "special" ghost points.
122:   */

124:   /* view each local ghosted vector */
125:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer));
126:   PetscCall(VecView(local, subviewer));
127:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &subviewer));

129:   PetscCall(VecScatterDestroy(&scatter1));
130:   PetscCall(VecScatterDestroy(&scatter2));
131:   PetscCall(VecDestroy(&local));
132:   PetscCall(VecDestroy(&global));
133:   PetscCall(DMDestroy(&da));
134:   PetscCall(PetscFinalize());
135:   return 0;
136: }

138: /*TEST

140:    test:

142:    test:
143:       suffix: 2
144:       nsize: 2

146:    test:
147:       suffix: 4
148:       nsize: 4

150:    test:
151:       suffix: 9
152:       nsize: 9

154: TEST*/