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