Actual source code: ex6.c
petsc-3.10.5 2019-03-28
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: 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;
49: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
51: /* Create distributed array and get vectors */
52: 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: DMSetUp(da);
54: DMCreateGlobalVector(da,&global);
55: DMCreateLocalVector(da,&local);
57: /* construct global to local scatter for the left side of the domain to the ghost on the bottom */
58: 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: ISCreateStride(PETSC_COMM_SELF,xm,1,1,&to);
61: } else {
62: ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);
63: }
64: PetscMalloc1(xm,&ifrom);
65: for (i=x;i<x+xm;i++) {
66: ifrom[i-x] = M*i;
67: }
68: DMDAGetAO(da,&ao);
69: AOApplicationToPetsc(ao,xm,ifrom);
70: if (!y) {
71: ISCreateGeneral(PETSC_COMM_WORLD,xm,ifrom,PETSC_OWN_POINTER,&from);
72: } else {
73: PetscFree(ifrom);
74: ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);
75: }
76: VecScatterCreate(global,from,local,to,&scatter1);
77: ISDestroy(&to);
78: ISDestroy(&from);
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: ISCreateStride(PETSC_COMM_SELF,ym,xm+2,xm+2,&to);
83: } else {
84: ISCreateStride(PETSC_COMM_SELF,0,0,0,&to);
85: }
86: PetscMalloc1(ym,&ifrom);
87: for (i=y;i<y+ym;i++) {
88: ifrom[i-y] = i;
89: }
90: DMDAGetAO(da,&ao);
91: AOApplicationToPetsc(ao,ym,ifrom);
92: if (!x) {
93: ISCreateGeneral(PETSC_COMM_WORLD,ym,ifrom,PETSC_OWN_POINTER,&from);
94: } else {
95: PetscFree(ifrom);
96: ISCreateGeneral(PETSC_COMM_WORLD,0,NULL,PETSC_COPY_VALUES,&from);
97: }
98: VecScatterCreate(global,from,local,to,&scatter2);
99: ISDestroy(&to);
100: ISDestroy(&from);
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: DMDACreateNaturalVector(da,&natural);
107: VecGetOwnershipRange(natural,&start,&end);
108: VecGetArray(natural,&xnatural);
109: for (i=start; i<end; i++) {
110: xnatural[i-start] = i;
111: }
112: VecRestoreArray(natural,&xnatural);
113: DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,global);
114: DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,global);
115: VecDestroy(&natural);
117: /* scatter from global to local */
118: VecScatterBegin(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);
119: VecScatterEnd(scatter1,global,local,INSERT_VALUES,SCATTER_FORWARD);
120: VecScatterBegin(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);
121: VecScatterEnd(scatter2,global,local,INSERT_VALUES,SCATTER_FORWARD);
122: /*
123: normally here you would also call
124: DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);
125: DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);
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: */
130: /* view each local ghosted vector */
131: PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);
132: VecView(local,subviewer);
133: PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&subviewer);
135: VecScatterDestroy(&scatter1);
136: VecScatterDestroy(&scatter2);
137: VecDestroy(&local);
138: VecDestroy(&global);
139: DMDestroy(&da);
140: PetscFinalize();
141: return ierr;
142: }
146: /*TEST
148: test:
150: test:
151: suffix: 2
152: nsize: 2
154: test:
155: suffix: 4
156: nsize: 4
158: test:
159: suffix: 9
160: nsize: 9
162: TEST*/