Actual source code: ex6.c

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

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

  7:     Consider the domain (with natural numbering)

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

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

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

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

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

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

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

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

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

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

 50:   /* Create distributed array and get vectors */
 51:   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));
 52:   PetscCall(DMSetUp(da));
 53:   PetscCall(DMCreateGlobalVector(da, &global));
 54:   PetscCall(DMCreateLocalVector(da, &local));

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

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

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

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

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

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

137: /*TEST

139:    test:

141:    test:
142:       suffix: 2
143:       nsize: 2

145:    test:
146:       suffix: 4
147:       nsize: 4

149:    test:
150:       suffix: 9
151:       nsize: 9

153: TEST*/