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