1 
2 static char help[] = "Tests DMCreateDomainDecomposition.\n\n";
3 
4 /*
5 Use the options
6      -da_grid_x <nx> - number of grid points in x direction, if M < 0
7      -da_grid_y <ny> - number of grid points in y direction, if N < 0
8      -da_processors_x <MX> number of processors in x directio
9      -da_processors_y <MY> number of processors in x direction
10 */
11 
12 #include <petscdm.h>
13 #include <petscdmda.h>
14 
FillLocalSubdomain(DM da,Vec gvec)15 PetscErrorCode FillLocalSubdomain(DM da, Vec gvec)
16 {
17   DMDALocalInfo  info;
18   PetscMPIInt    rank;
19   PetscInt       i,j,k,l;
20   PetscErrorCode ierr;
21 
22   PetscFunctionBeginUser;
23   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
24   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
25 
26   if (info.dim == 3) {
27     PetscScalar    ***g;
28     ierr = DMDAVecGetArray(da,gvec,&g);CHKERRQ(ierr);
29     /* loop over ghosts */
30     for (k=info.zs; k<info.zs+info.zm; k++) {
31       for (j=info.ys; j<info.ys+info.ym; j++) {
32         for (i=info.xs; i<info.xs+info.xm; i++) {
33           g[k][j][info.dof*i+0]   = i;
34           g[k][j][info.dof*i+1]   = j;
35           g[k][j][info.dof*i+2]   = k;
36         }
37       }
38     }
39     ierr = DMDAVecRestoreArray(da,gvec,&g);CHKERRQ(ierr);
40   }
41   if (info.dim == 2) {
42     PetscScalar    **g;
43     ierr = DMDAVecGetArray(da,gvec,&g);CHKERRQ(ierr);
44     /* loop over ghosts */
45     for (j=info.ys; j<info.ys+info.ym; j++) {
46       for (i=info.xs; i<info.xs+info.xm; i++) {
47         for (l = 0;l<info.dof;l++) {
48           g[j][info.dof*i+0]   = i;
49           g[j][info.dof*i+1]   = j;
50           g[j][info.dof*i+2]   = rank;
51         }
52       }
53     }
54     ierr = DMDAVecRestoreArray(da,gvec,&g);CHKERRQ(ierr);
55   }
56   PetscFunctionReturn(0);
57 }
58 
main(int argc,char ** argv)59 int main(int argc,char **argv)
60 {
61   PetscErrorCode ierr;
62   DM             da,*subda;
63   PetscInt       i,dim=3;
64   PetscMPIInt    size,rank;
65   Vec            v;
66   Vec            slvec,sgvec;
67   IS             *ois,*iis;
68   VecScatter     oscata;
69   VecScatter     *iscat,*oscat,*gscat;
70   DMDALocalInfo  info;
71 
72   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
73   ierr = PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);CHKERRQ(ierr);
74 
75   /* Create distributed array and get vectors */
76   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
77   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
78   if (dim == 2) {
79     ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);CHKERRQ(ierr);
80   } else if (dim == 3) {
81     ierr = DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,-4,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,NULL,&da);CHKERRQ(ierr);
82   }
83   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
84   ierr = DMSetUp(da);CHKERRQ(ierr);
85   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
86 
87   ierr = DMCreateDomainDecomposition(da,NULL,NULL,&iis,&ois,&subda);CHKERRQ(ierr);
88   ierr = DMCreateDomainDecompositionScatters(da,1,subda,&iscat,&oscat,&gscat);CHKERRQ(ierr);
89 
90   {
91     DMDALocalInfo subinfo;
92     MatStencil    lower,upper;
93     IS            patchis,subpatchis;
94     Vec           smallvec;
95     Vec           largevec;
96     VecScatter    patchscat;
97 
98     ierr = DMDAGetLocalInfo(subda[0],&subinfo);CHKERRQ(ierr);
99 
100     lower.i = info.xs;
101     lower.j = info.ys;
102     lower.k = info.zs;
103     upper.i = info.xs+info.xm;
104     upper.j = info.ys+info.ym;
105     upper.k = info.zs+info.zm;
106 
107     /* test the patch IS as a thing to scatter to/from */
108     ierr = DMDACreatePatchIS(da,&lower,&upper,&patchis);CHKERRQ(ierr);
109     ierr = DMGetGlobalVector(da,&largevec);CHKERRQ(ierr);
110 
111     ierr = VecCreate(PETSC_COMM_SELF,&smallvec);CHKERRQ(ierr);
112     ierr = VecSetSizes(smallvec,info.dof*(upper.i - lower.i)*(upper.j - lower.j)*(upper.k - lower.k),PETSC_DECIDE);CHKERRQ(ierr);
113     ierr = VecSetFromOptions(smallvec);CHKERRQ(ierr);
114     ierr = VecScatterCreate(smallvec,NULL,largevec,patchis,&patchscat);CHKERRQ(ierr);
115 
116     ierr = FillLocalSubdomain(subda[0],smallvec);CHKERRQ(ierr);
117     ierr = VecSet(largevec,0);CHKERRQ(ierr);
118 
119     ierr = VecScatterBegin(patchscat,smallvec,largevec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
120     ierr = VecScatterEnd(patchscat,smallvec,largevec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
121 
122     for (i = 0; i < size; i++) {
123       if (i == rank) {
124         ierr = ISView(patchis,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
125         ierr = VecScatterView(patchscat,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
126         ierr = VecView(smallvec,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
127       }
128       ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
129     }
130 
131     ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
132     ierr = VecView(largevec,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
133 
134     ierr = VecDestroy(&smallvec);CHKERRQ(ierr);
135     ierr = DMRestoreGlobalVector(da,&largevec);CHKERRQ(ierr);
136     ierr = ISDestroy(&patchis);CHKERRQ(ierr);
137     ierr = VecScatterDestroy(&patchscat);CHKERRQ(ierr);
138   }
139 
140   /* view the various parts */
141   {
142     for (i = 0; i < size; i++) {
143       if (i == rank) {
144         ierr = PetscPrintf(PETSC_COMM_SELF,"Processor %d: \n",i);CHKERRQ(ierr);
145         ierr = DMView(subda[0],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
146       }
147       ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
148     }
149 
150     ierr = DMGetLocalVector(subda[0],&slvec);CHKERRQ(ierr);
151     ierr = DMGetGlobalVector(subda[0],&sgvec);CHKERRQ(ierr);
152     ierr = DMGetGlobalVector(da,&v);CHKERRQ(ierr);
153 
154     /* test filling outer between the big DM and the small ones with the IS scatter*/
155     ierr = VecScatterCreate(v,ois[0],sgvec,NULL,&oscata);CHKERRQ(ierr);
156 
157     ierr = FillLocalSubdomain(subda[0],sgvec);CHKERRQ(ierr);
158 
159     ierr = VecScatterBegin(oscata,sgvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
160     ierr = VecScatterEnd(oscata,sgvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
161 
162     /* test the local-to-local scatter */
163 
164     /* fill up the local subdomain and then add them together */
165     ierr = FillLocalSubdomain(da,v);CHKERRQ(ierr);
166 
167     ierr = VecScatterBegin(gscat[0],v,slvec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
168     ierr = VecScatterEnd(gscat[0],v,slvec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
169 
170     ierr = VecView(v,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
171 
172     /* test ghost scattering backwards */
173 
174     ierr = VecSet(v,0);CHKERRQ(ierr);
175 
176     ierr = VecScatterBegin(gscat[0],slvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
177     ierr = VecScatterEnd(gscat[0],slvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
178 
179     ierr = VecView(v,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
180 
181     /* test overlap scattering backwards */
182 
183     ierr = DMLocalToGlobalBegin(subda[0],slvec,ADD_VALUES,sgvec);CHKERRQ(ierr);
184     ierr = DMLocalToGlobalEnd(subda[0],slvec,ADD_VALUES,sgvec);CHKERRQ(ierr);
185 
186     ierr = VecSet(v,0);CHKERRQ(ierr);
187 
188     ierr = VecScatterBegin(oscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
189     ierr = VecScatterEnd(oscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
190 
191     ierr = VecView(v,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
192 
193     /* test interior scattering backwards */
194 
195     ierr = VecSet(v,0);CHKERRQ(ierr);
196 
197     ierr = VecScatterBegin(iscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
198     ierr = VecScatterEnd(iscat[0],sgvec,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
199 
200     ierr = VecView(v,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
201 
202     /* test matrix allocation */
203     for (i = 0; i < size; i++) {
204       if (i == rank) {
205         Mat m;
206         ierr = PetscPrintf(PETSC_COMM_SELF,"Processor %d: \n",i);CHKERRQ(ierr);
207         ierr = DMSetMatType(subda[0],MATAIJ);CHKERRQ(ierr);
208         ierr = DMCreateMatrix(subda[0],&m);CHKERRQ(ierr);
209         ierr = MatView(m,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
210         ierr = MatDestroy(&m);CHKERRQ(ierr);
211       }
212       ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRQ(ierr);
213     }
214     ierr = DMRestoreLocalVector(subda[0],&slvec);CHKERRQ(ierr);
215     ierr = DMRestoreGlobalVector(subda[0],&sgvec);CHKERRQ(ierr);
216     ierr = DMRestoreGlobalVector(da,&v);CHKERRQ(ierr);
217   }
218 
219 
220   ierr = DMDestroy(&subda[0]);CHKERRQ(ierr);
221   ierr = ISDestroy(&ois[0]);CHKERRQ(ierr);
222   ierr = ISDestroy(&iis[0]);CHKERRQ(ierr);
223 
224   ierr = VecScatterDestroy(&iscat[0]);CHKERRQ(ierr);
225   ierr = VecScatterDestroy(&oscat[0]);CHKERRQ(ierr);
226   ierr = VecScatterDestroy(&gscat[0]);CHKERRQ(ierr);
227   ierr = VecScatterDestroy(&oscata);CHKERRQ(ierr);
228 
229   ierr = PetscFree(iscat);CHKERRQ(ierr);
230   ierr = PetscFree(oscat);CHKERRQ(ierr);
231   ierr = PetscFree(gscat);CHKERRQ(ierr);
232   ierr = PetscFree(oscata);CHKERRQ(ierr);
233 
234   ierr = PetscFree(subda);CHKERRQ(ierr);
235   ierr = PetscFree(ois);CHKERRQ(ierr);
236   ierr = PetscFree(iis);CHKERRQ(ierr);
237 
238   ierr = DMDestroy(&da);CHKERRQ(ierr);
239   ierr = PetscFinalize();
240   return ierr;
241 }
242 
243