1
2 static char help[] = "Tests DMSwarm\n\n";
3
4 #include <petscdm.h>
5 #include <petscdmda.h>
6 #include <petscdmswarm.h>
7
8 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM,PetscErrorCode (*)(DM,void*,PetscInt*,PetscInt**),size_t,void*,PetscInt*);
9 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM,PetscInt*);
10
ex1_1(void)11 PetscErrorCode ex1_1(void)
12 {
13 DM dms;
14 PetscErrorCode ierr;
15 Vec x;
16 PetscMPIInt rank,size;
17 PetscInt p;
18
19 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
20 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
21 if ((size > 1) && (size != 4)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Must be run wuth 4 MPI ranks");
22
23 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr);
24 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr);
25
26 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr);
27 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr);
28 ierr = DMSwarmRegisterPetscDatatypeField(dms,"strain",1,PETSC_REAL);CHKERRQ(ierr);
29 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr);
30 ierr = DMSwarmSetLocalSizes(dms,5+rank,4);CHKERRQ(ierr);
31 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
32
33 {
34 PetscReal *array;
35 ierr = DMSwarmGetField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr);
36 for (p=0; p<5+rank; p++) {
37 array[p] = 11.1 + p*0.1 + rank*100.0;
38 }
39 ierr = DMSwarmRestoreField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr);
40 }
41
42 {
43 PetscReal *array;
44 ierr = DMSwarmGetField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr);
45 for (p=0; p<5+rank; p++) {
46 array[p] = 2.0e-2 + p*0.001 + rank*1.0;
47 }
48 ierr = DMSwarmRestoreField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr);
49 }
50
51 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
52 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
53
54 ierr = DMSwarmVectorDefineField(dms,"strain");CHKERRQ(ierr);
55 ierr = DMCreateGlobalVector(dms,&x);CHKERRQ(ierr);
56 ierr = VecDestroy(&x);CHKERRQ(ierr);
57
58 {
59 PetscInt *rankval;
60 PetscInt npoints[2],npoints_orig[2];
61
62 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr);
63 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr);
64 ierr = DMSwarmGetField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
65 if ((rank == 0) && (size > 1)) {
66 rankval[0] = 1;
67 rankval[3] = 1;
68 }
69 if (rank == 3) {
70 rankval[2] = 1;
71 }
72 ierr = DMSwarmRestoreField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
73 ierr = DMSwarmMigrate(dms,PETSC_TRUE);CHKERRQ(ierr);
74 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr);
75 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr);
76 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
77 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D) after(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1],npoints[0],npoints[1]);CHKERRQ(ierr);
78 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
79 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
80 }
81 {
82 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
83 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
84 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
85 }
86 {
87 ierr = DMSwarmCreateGlobalVectorFromField(dms,"strain",&x);CHKERRQ(ierr);
88 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
89 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"strain",&x);CHKERRQ(ierr);
90 }
91
92 ierr = DMDestroy(&dms);CHKERRQ(ierr);
93 PetscFunctionReturn(0);
94 }
95
ex1_2(void)96 PetscErrorCode ex1_2(void)
97 {
98 DM dms;
99 PetscErrorCode ierr;
100 Vec x;
101 PetscMPIInt rank,size;
102 PetscInt p;
103
104 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
105 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
106 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr);
107 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr);
108 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr);
109
110 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr);
111 ierr = DMSwarmRegisterPetscDatatypeField(dms,"strain",1,PETSC_REAL);CHKERRQ(ierr);
112 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr);
113 ierr = DMSwarmSetLocalSizes(dms,5+rank,4);CHKERRQ(ierr);
114 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
115 {
116 PetscReal *array;
117 ierr = DMSwarmGetField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr);
118 for (p=0; p<5+rank; p++) {
119 array[p] = 11.1 + p*0.1 + rank*100.0;
120 }
121 ierr = DMSwarmRestoreField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr);
122 }
123 {
124 PetscReal *array;
125 ierr = DMSwarmGetField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr);
126 for (p=0; p<5+rank; p++) {
127 array[p] = 2.0e-2 + p*0.001 + rank*1.0;
128 }
129 ierr = DMSwarmRestoreField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr);
130 }
131 {
132 PetscInt *rankval;
133 PetscInt npoints[2],npoints_orig[2];
134
135 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr);
136 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr);
137 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
138 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1]);CHKERRQ(ierr);
139 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
140 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
141
142 ierr = DMSwarmGetField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
143
144 if (rank == 1) {
145 rankval[0] = -1;
146 }
147 if (rank == 2) {
148 rankval[1] = -1;
149 }
150 if (rank == 3) {
151 rankval[3] = -1;
152 rankval[4] = -1;
153 }
154 ierr = DMSwarmRestoreField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
155 ierr = DMSwarmCollectViewCreate(dms);CHKERRQ(ierr);
156 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr);
157 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr);
158 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
159 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr);
160 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
161 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
162
163 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
164 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
165 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
166
167 ierr = DMSwarmCollectViewDestroy(dms);CHKERRQ(ierr);
168 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr);
169 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr);
170 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
171 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after_v(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr);
172 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
173 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
174
175 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
176 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
177 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr);
178 }
179 ierr = DMDestroy(&dms);CHKERRQ(ierr);
180 PetscFunctionReturn(0);
181 }
182
183 /*
184 splot "c-rank0.gp","c-rank1.gp","c-rank2.gp","c-rank3.gp"
185 */
ex1_3(void)186 PetscErrorCode ex1_3(void)
187 {
188 DM dms;
189 PetscErrorCode ierr;
190 PetscMPIInt rank,size;
191 PetscInt is,js,ni,nj,overlap;
192 DM dmcell;
193
194 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
195 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
196 overlap = 2;
197 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,13,13,PETSC_DECIDE,PETSC_DECIDE,1,overlap,NULL,NULL,&dmcell);CHKERRQ(ierr);
198 ierr = DMSetFromOptions(dmcell);CHKERRQ(ierr);
199 ierr = DMSetUp(dmcell);CHKERRQ(ierr);
200 ierr = DMDASetUniformCoordinates(dmcell,-1.0,1.0,-1.0,1.0,-1.0,1.0);CHKERRQ(ierr);
201 ierr = DMDAGetCorners(dmcell,&is,&js,NULL,&ni,&nj,NULL);CHKERRQ(ierr);
202 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr);
203 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr);
204 ierr = DMSwarmSetCellDM(dms,dmcell);CHKERRQ(ierr);
205
206 /* load in data types */
207 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr);
208 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr);
209 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coorx",1,PETSC_REAL);CHKERRQ(ierr);
210 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coory",1,PETSC_REAL);CHKERRQ(ierr);
211 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr);
212 ierr = DMSwarmSetLocalSizes(dms,ni*nj*4,4);CHKERRQ(ierr);
213 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
214
215 /* set values within the swarm */
216 {
217 PetscReal *array_x,*array_y;
218 PetscInt npoints,i,j,cnt;
219 DMDACoor2d **LA_coor;
220 Vec coor;
221 DM dmcellcdm;
222
223 ierr = DMGetCoordinateDM(dmcell,&dmcellcdm);CHKERRQ(ierr);
224 ierr = DMGetCoordinates(dmcell,&coor);CHKERRQ(ierr);
225 ierr = DMDAVecGetArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr);
226 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr);
227 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
228 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
229 cnt = 0;
230 for (j=js; j<js+nj; j++) {
231 for (i=is; i<is+ni; i++) {
232 PetscReal xp,yp;
233 xp = PetscRealPart(LA_coor[j][i].x);
234 yp = PetscRealPart(LA_coor[j][i].y);
235 array_x[4*cnt+0] = xp - 0.05; if (array_x[4*cnt+0] < -1.0) { array_x[4*cnt+0] = -1.0+1.0e-12; }
236 array_x[4*cnt+1] = xp + 0.05; if (array_x[4*cnt+1] > 1.0) { array_x[4*cnt+1] = 1.0-1.0e-12; }
237 array_x[4*cnt+2] = xp - 0.05; if (array_x[4*cnt+2] < -1.0) { array_x[4*cnt+2] = -1.0+1.0e-12; }
238 array_x[4*cnt+3] = xp + 0.05; if (array_x[4*cnt+3] > 1.0) { array_x[4*cnt+3] = 1.0-1.0e-12; }
239
240 array_y[4*cnt+0] = yp - 0.05; if (array_y[4*cnt+0] < -1.0) { array_y[4*cnt+0] = -1.0+1.0e-12; }
241 array_y[4*cnt+1] = yp - 0.05; if (array_y[4*cnt+1] < -1.0) { array_y[4*cnt+1] = -1.0+1.0e-12; }
242 array_y[4*cnt+2] = yp + 0.05; if (array_y[4*cnt+2] > 1.0) { array_y[4*cnt+2] = 1.0-1.0e-12; }
243 array_y[4*cnt+3] = yp + 0.05; if (array_y[4*cnt+3] > 1.0) { array_y[4*cnt+3] = 1.0-1.0e-12; }
244 cnt++;
245 }
246 }
247 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
248 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
249 ierr = DMDAVecRestoreArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr);
250 }
251 {
252 PetscInt npoints[2],npoints_orig[2],ng;
253
254 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr);
255 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr);
256 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
257 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1]);CHKERRQ(ierr);
258 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
259 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
260 ierr = DMSwarmCollect_DMDABoundingBox(dms,&ng);CHKERRQ(ierr);
261
262 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr);
263 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr);
264 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
265 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr);
266 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
267 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
268 }
269 {
270 PetscReal *array_x,*array_y;
271 PetscInt npoints,p;
272 FILE *fp = NULL;
273 char name[PETSC_MAX_PATH_LEN];
274
275 ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"c-rank%d.gp",rank);CHKERRQ(ierr);
276 fp = fopen(name,"w");
277 if (!fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name);
278 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr);
279 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
280 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
281 for (p=0; p<npoints; p++) {
282 fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array_x[p],array_y[p],(double)rank);
283 }
284 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
285 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
286 fclose(fp);
287 }
288 ierr = DMDestroy(&dmcell);CHKERRQ(ierr);
289 ierr = DMDestroy(&dms);CHKERRQ(ierr);
290 PetscFunctionReturn(0);
291 }
292
293 typedef struct {
294 PetscReal cx[2];
295 PetscReal radius;
296 } CollectZoneCtx;
297
collect_zone(DM dm,void * ctx,PetscInt * nfound,PetscInt ** foundlist)298 PetscErrorCode collect_zone(DM dm,void *ctx,PetscInt *nfound,PetscInt **foundlist)
299 {
300 CollectZoneCtx *zone = (CollectZoneCtx*)ctx;
301 PetscInt p,npoints;
302 PetscReal *array_x,*array_y,r2;
303 PetscInt p2collect,*plist;
304 PetscMPIInt rank;
305 PetscErrorCode ierr;
306
307 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
308 ierr = DMSwarmGetLocalSize(dm,&npoints);CHKERRQ(ierr);
309 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
310 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
311
312 r2 = zone->radius * zone->radius;
313 p2collect = 0;
314 for (p=0; p<npoints; p++) {
315 PetscReal sep2;
316
317 sep2 = (array_x[p] - zone->cx[0])*(array_x[p] - zone->cx[0]);
318 sep2 += (array_y[p] - zone->cx[1])*(array_y[p] - zone->cx[1]);
319 if (sep2 < r2) {
320 p2collect++;
321 }
322 }
323
324 ierr = PetscMalloc1(p2collect+1,&plist);CHKERRQ(ierr);
325 p2collect = 0;
326 for (p=0; p<npoints; p++) {
327 PetscReal sep2;
328
329 sep2 = (array_x[p] - zone->cx[0])*(array_x[p] - zone->cx[0]);
330 sep2 += (array_y[p] - zone->cx[1])*(array_y[p] - zone->cx[1]);
331 if (sep2 < r2) {
332 plist[p2collect] = p;
333 p2collect++;
334 }
335 }
336 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
337 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
338
339 *nfound = p2collect;
340 *foundlist = plist;
341 PetscFunctionReturn(0);
342 }
343
ex1_4(void)344 PetscErrorCode ex1_4(void)
345 {
346 DM dms;
347 PetscErrorCode ierr;
348 PetscMPIInt rank,size;
349 PetscInt is,js,ni,nj,overlap,nn;
350 DM dmcell;
351 CollectZoneCtx *zone;
352 PetscReal dx;
353
354 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
355 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
356 nn = 101;
357 dx = 2.0/ (PetscReal)(nn-1);
358 overlap = 0;
359 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nn,nn,PETSC_DECIDE,PETSC_DECIDE,1,overlap,NULL,NULL,&dmcell);CHKERRQ(ierr);
360 ierr = DMSetFromOptions(dmcell);CHKERRQ(ierr);
361 ierr = DMSetUp(dmcell);CHKERRQ(ierr);
362 ierr = DMDASetUniformCoordinates(dmcell,-1.0,1.0,-1.0,1.0,-1.0,1.0);CHKERRQ(ierr);
363 ierr = DMDAGetCorners(dmcell,&is,&js,NULL,&ni,&nj,NULL);CHKERRQ(ierr);
364 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr);
365 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr);
366
367 /* load in data types */
368 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr);
369 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr);
370 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coorx",1,PETSC_REAL);CHKERRQ(ierr);
371 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coory",1,PETSC_REAL);CHKERRQ(ierr);
372 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr);
373 ierr = DMSwarmSetLocalSizes(dms,ni*nj*4,4);CHKERRQ(ierr);
374 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
375
376 /* set values within the swarm */
377 {
378 PetscReal *array_x,*array_y;
379 PetscInt npoints,i,j,cnt;
380 DMDACoor2d **LA_coor;
381 Vec coor;
382 DM dmcellcdm;
383
384 ierr = DMGetCoordinateDM(dmcell,&dmcellcdm);CHKERRQ(ierr);
385 ierr = DMGetCoordinates(dmcell,&coor);CHKERRQ(ierr);
386 ierr = DMDAVecGetArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr);
387 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr);
388 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
389 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
390 cnt = 0;
391 for (j=js; j<js+nj; j++) {
392 for (i=is; i<is+ni; i++) {
393 PetscReal xp,yp;
394
395 xp = PetscRealPart(LA_coor[j][i].x);
396 yp = PetscRealPart(LA_coor[j][i].y);
397 array_x[4*cnt+0] = xp - dx*0.1; /*if (array_x[4*cnt+0] < -1.0) { array_x[4*cnt+0] = -1.0+1.0e-12; }*/
398 array_x[4*cnt+1] = xp + dx*0.1; /*if (array_x[4*cnt+1] > 1.0) { array_x[4*cnt+1] = 1.0-1.0e-12; }*/
399 array_x[4*cnt+2] = xp - dx*0.1; /*if (array_x[4*cnt+2] < -1.0) { array_x[4*cnt+2] = -1.0+1.0e-12; }*/
400 array_x[4*cnt+3] = xp + dx*0.1; /*if (array_x[4*cnt+3] > 1.0) { array_x[4*cnt+3] = 1.0-1.0e-12; }*/
401 array_y[4*cnt+0] = yp - dx*0.1; /*if (array_y[4*cnt+0] < -1.0) { array_y[4*cnt+0] = -1.0+1.0e-12; }*/
402 array_y[4*cnt+1] = yp - dx*0.1; /*if (array_y[4*cnt+1] < -1.0) { array_y[4*cnt+1] = -1.0+1.0e-12; }*/
403 array_y[4*cnt+2] = yp + dx*0.1; /*if (array_y[4*cnt+2] > 1.0) { array_y[4*cnt+2] = 1.0-1.0e-12; }*/
404 array_y[4*cnt+3] = yp + dx*0.1; /*if (array_y[4*cnt+3] > 1.0) { array_y[4*cnt+3] = 1.0-1.0e-12; }*/
405 cnt++;
406 }
407 }
408 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
409 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
410 ierr = DMDAVecRestoreArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr);
411 }
412 ierr = PetscMalloc1(1,&zone);CHKERRQ(ierr);
413 if (size == 4) {
414 if (rank == 0) {
415 zone->cx[0] = 0.5;
416 zone->cx[1] = 0.5;
417 zone->radius = 0.3;
418 }
419 if (rank == 1) {
420 zone->cx[0] = -0.5;
421 zone->cx[1] = 0.5;
422 zone->radius = 0.25;
423 }
424 if (rank == 2) {
425 zone->cx[0] = 0.5;
426 zone->cx[1] = -0.5;
427 zone->radius = 0.2;
428 }
429 if (rank == 3) {
430 zone->cx[0] = -0.5;
431 zone->cx[1] = -0.5;
432 zone->radius = 0.1;
433 }
434 } else {
435 if (rank == 0) {
436 zone->cx[0] = 0.5;
437 zone->cx[1] = 0.5;
438 zone->radius = 0.8;
439 } else {
440 zone->cx[0] = 10.0;
441 zone->cx[1] = 10.0;
442 zone->radius = 0.0;
443 }
444 }
445 {
446 PetscInt npoints[2],npoints_orig[2],ng;
447
448 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr);
449 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr);
450 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
451 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1]);CHKERRQ(ierr);
452 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
453 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
454 ierr = DMSwarmCollect_General(dms,collect_zone,sizeof(CollectZoneCtx),zone,&ng);CHKERRQ(ierr);
455 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr);
456 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr);
457 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
458 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr);
459 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
460 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
461 }
462 {
463 PetscReal *array_x,*array_y;
464 PetscInt npoints,p;
465 FILE *fp = NULL;
466 char name[PETSC_MAX_PATH_LEN];
467
468 ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"c-rank%d.gp",rank);CHKERRQ(ierr);
469 fp = fopen(name,"w");
470 if (!fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name);
471 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr);
472 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
473 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
474 for (p=0; p<npoints; p++) {
475 fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array_x[p],array_y[p],(double)rank);
476 }
477 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
478 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
479 fclose(fp);
480 }
481 ierr = DMDestroy(&dmcell);CHKERRQ(ierr);
482 ierr = DMDestroy(&dms);CHKERRQ(ierr);
483 ierr = PetscFree(zone);CHKERRQ(ierr);
484 PetscFunctionReturn(0);
485 }
486
main(int argc,char ** argv)487 int main(int argc,char **argv)
488 {
489 PetscErrorCode ierr;
490 PetscInt test_mode = 4;
491
492 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
493 ierr = PetscOptionsGetInt(NULL,NULL,"-test_mode",&test_mode,NULL);CHKERRQ(ierr);
494 if (test_mode == 1) {
495 ierr = ex1_1();CHKERRQ(ierr);
496 } else if (test_mode == 2) {
497 ierr = ex1_2();CHKERRQ(ierr);
498 } else if (test_mode == 3) {
499 ierr = ex1_3();CHKERRQ(ierr);
500 } else if (test_mode == 4) {
501 ierr = ex1_4();CHKERRQ(ierr);
502 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown test_mode value, should be 1,2,3,4");
503 ierr = PetscFinalize();
504 return ierr;
505 }
506
507 /*TEST
508
509 build:
510 requires: !complex double
511
512 test:
513 args: -test_mode 1
514 filter: grep -v atomic
515 filter_output: grep -v atomic
516
517 test:
518 suffix: 2
519 args: -test_mode 2
520 filter: grep -v atomic
521 filter_output: grep -v atomic
522
523 test:
524 suffix: 3
525 args: -test_mode 3
526 filter: grep -v atomic
527 filter_output: grep -v atomic
528 TODO: broken
529
530 test:
531 suffix: 4
532 args: -test_mode 4
533 filter: grep -v atomic
534 filter_output: grep -v atomic
535
536 test:
537 suffix: 5
538 nsize: 4
539 args: -test_mode 1
540 filter: grep -v atomic
541 filter_output: grep -v atomic
542
543 test:
544 suffix: 6
545 nsize: 4
546 args: -test_mode 2
547 filter: grep -v atomic
548 filter_output: grep -v atomic
549
550 test:
551 suffix: 7
552 nsize: 4
553 args: -test_mode 3
554 filter: grep -v atomic
555 filter_output: grep -v atomic
556 TODO: broken
557
558 test:
559 suffix: 8
560 nsize: 4
561 args: -test_mode 4
562 filter: grep -v atomic
563 filter_output: grep -v atomic
564
565 TEST*/
566