1
2 static char help[] = "\n\n";
3
4 /*
5 Demonstrates using DM_BOUNDARY_GHOSTED how to handle a rotated boundary conditions where one edge
6 is connected to its immediate neighbor
7
8 Consider the domain (with natural numbering)
9
10 6 7 8
11 3 4 5
12 0 1 2
13
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
16
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)
19
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
25
26 where x indicates a location that is not updated by the communication and should be used.
27
28 For this to make sense the number of grid points in the x and y directions must be the same
29
30 This ghost point mapping was suggested by: Wenbo Zhao <zhaowenbo.npic@gmail.com>
31 */
32
33 #include <petscdm.h>
34 #include <petscdmda.h>
35
main(int argc,char ** argv)36 int main(int argc,char **argv)
37 {
38 PetscInt M = 6;
39 PetscErrorCode ierr;
40 DM da;
41 Vec local,global,natural;
42 PetscInt i,start,end,*ifrom,x,y,xm,ym;
43 PetscScalar *xnatural;
44 IS from,to;
45 AO ao;
46 VecScatter scatter1,scatter2;
47 PetscViewer subviewer;
48
49 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
50
51 /* Create distributed array and get vectors */
52 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,DM_BOUNDARY_GHOSTED,DMDA_STENCIL_STAR,M,M,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr);
53 ierr = DMSetUp(da);CHKERRQ(ierr);
54 ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);
55 ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);
56
57 /* construct global to local scatter for the left side of the domain to the ghost on the bottom */
58 ierr = DMDAGetCorners(da,&x,&y,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
59 if (!y) { /* only processes on the bottom of the domain fill up the ghost locations */
60 ierr = ISCreateStride(PETSC_COMM_SELF,xm,1,1,&to);CHKERRQ(ierr);
61 } else {
62 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);CHKERRQ(ierr);
63 }
64 ierr = PetscMalloc1(xm,&ifrom);CHKERRQ(ierr);
65 for (i=x;i<x+xm;i++) {
66 ifrom[i-x] = M*i;
67 }
68 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr);
69 ierr = AOApplicationToPetsc(ao,xm,ifrom);CHKERRQ(ierr);
70 if (!y) {
71 ierr = ISCreateGeneral(PETSC_COMM_WORLD,xm,ifrom,PETSC_OWN_POINTER,&from);CHKERRQ(ierr);
72 } else {
73 ierr = PetscFree(ifrom);CHKERRQ(ierr);
74 ierr = ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
75 }
76 ierr = VecScatterCreate(global,from,local,to,&scatter1);CHKERRQ(ierr);
77 ierr = ISDestroy(&to);CHKERRQ(ierr);
78 ierr = ISDestroy(&from);CHKERRQ(ierr);
79
80 /* construct global to local scatter for the bottom side of the domain to the ghost on the right */
81 if (!x) { /* only processes on the left side of the domain fill up the ghost locations */
82 ierr = ISCreateStride(PETSC_COMM_SELF,ym,xm+2,xm+2,&to);CHKERRQ(ierr);
83 } else {
84 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);CHKERRQ(ierr);
85 }
86 ierr = PetscMalloc1(ym,&ifrom);CHKERRQ(ierr);
87 for (i=y;i<y+ym;i++) {
88 ifrom[i-y] = i;
89 }
90 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr);
91 ierr = AOApplicationToPetsc(ao,ym,ifrom);CHKERRQ(ierr);
92 if (!x) {
93 ierr = ISCreateGeneral(PETSC_COMM_WORLD,ym,ifrom,PETSC_OWN_POINTER,&from);CHKERRQ(ierr);
94 } else {
95 ierr = PetscFree(ifrom);CHKERRQ(ierr);
96 ierr = ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
97 }
98 ierr = VecScatterCreate(global,from,local,to,&scatter2);CHKERRQ(ierr);
99 ierr = ISDestroy(&to);CHKERRQ(ierr);
100 ierr = ISDestroy(&from);CHKERRQ(ierr);
101
102 /*
103 fill the global vector with the natural global numbering for each local entry
104 this is only done for testing purposes since it is easy to see if the scatter worked correctly
105 */
106 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
107 ierr = VecGetOwnershipRange(natural,&start,&end);CHKERRQ(ierr);
108 ierr = VecGetArray(natural,&xnatural);CHKERRQ(ierr);
109 for (i=start; i<end; i++) {
110 xnatural[i-start] = i;
111 }
112 ierr = VecRestoreArray(natural,&xnatural);CHKERRQ(ierr);
113 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,global);CHKERRQ(ierr);
114 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,global);CHKERRQ(ierr);
115 ierr = VecDestroy(&natural);CHKERRQ(ierr);
116
117 /* scatter from global to local */
118 ierr = VecScatterBegin(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
119 ierr = VecScatterEnd(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
120 ierr = VecScatterBegin(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
121 ierr = VecScatterEnd(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
122 /*
123 normally here you would also call
124 ierr = DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
125 ierr = DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);CHKERRQ(ierr);
126 to update all the interior ghost cells between neighboring processes.
127 We don't do it here since this is only a test of "special" ghost points.
128 */
129
130 /* view each local ghosted vector */
131 ierr = PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);CHKERRQ(ierr);
132 ierr = VecView(local,subviewer);CHKERRQ(ierr);
133 ierr = PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);CHKERRQ(ierr);
134
135 ierr = VecScatterDestroy(&scatter1);CHKERRQ(ierr);
136 ierr = VecScatterDestroy(&scatter2);CHKERRQ(ierr);
137 ierr = VecDestroy(&local);CHKERRQ(ierr);
138 ierr = VecDestroy(&global);CHKERRQ(ierr);
139 ierr = DMDestroy(&da);CHKERRQ(ierr);
140 ierr = PetscFinalize();
141 return ierr;
142 }
143
144
145
146 /*TEST
147
148 test:
149
150 test:
151 suffix: 2
152 nsize: 2
153
154 test:
155 suffix: 4
156 nsize: 4
157
158 test:
159 suffix: 9
160 nsize: 9
161
162 TEST*/
163