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*/