Actual source code: swarm_migrate.c
1: #include <petscsf.h>
2: #include <petscdmswarm.h>
3: #include <petscdmda.h>
4: #include <petsc/private/dmswarmimpl.h>
5: #include "../src/dm/impls/swarm/data_bucket.h"
6: #include "../src/dm/impls/swarm/data_ex.h"
8: /*
9: User loads desired location (MPI rank) into field DMSwarm_rank
10: */
11: PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm, PetscBool remove_sent_points)
12: {
13: DM_Swarm *swarm = (DM_Swarm *)dm->data;
14: DMSwarmDataEx de;
15: PetscInt p, npoints, *rankval, n_points_recv;
16: PetscMPIInt rank, nrank;
17: void *point_buffer, *recv_points;
18: size_t sizeof_dmswarm_point;
19: PetscBool debug = PETSC_FALSE;
21: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
23: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
24: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
25: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
26: DMSwarmDataExTopologyInitialize(de);
27: for (p = 0; p < npoints; ++p) {
28: nrank = rankval[p];
29: if (nrank != rank) DMSwarmDataExTopologyAddNeighbour(de, nrank);
30: }
31: DMSwarmDataExTopologyFinalize(de);
32: DMSwarmDataExInitializeSendCount(de);
33: for (p = 0; p < npoints; p++) {
34: nrank = rankval[p];
35: if (nrank != rank) DMSwarmDataExAddToSendCount(de, nrank, 1);
36: }
37: DMSwarmDataExFinalizeSendCount(de);
38: DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
39: DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
40: for (p = 0; p < npoints; p++) {
41: nrank = rankval[p];
42: if (nrank != rank) {
43: /* copy point into buffer */
44: DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
45: /* insert point buffer into DMSwarmDataExchanger */
46: DMSwarmDataExPackData(de, nrank, 1, point_buffer);
47: }
48: }
49: DMSwarmDataExPackFinalize(de);
50: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
52: if (remove_sent_points) {
53: DMSwarmDataField gfield;
55: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &gfield);
56: DMSwarmDataFieldGetAccess(gfield);
57: DMSwarmDataFieldGetEntries(gfield, (void **)&rankval);
59: /* remove points which left processor */
60: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
61: for (p = 0; p < npoints; p++) {
62: nrank = rankval[p];
63: if (nrank != rank) {
64: /* kill point */
65: DMSwarmDataFieldRestoreAccess(gfield);
67: DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);
69: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL); /* you need to update npoints as the list size decreases! */
70: DMSwarmDataFieldGetAccess(gfield);
71: DMSwarmDataFieldGetEntries(gfield, (void **)&rankval);
72: p--; /* check replacement point */
73: }
74: }
75: DMSwarmDataFieldRestoreEntries(gfield, (void **)&rankval);
76: DMSwarmDataFieldRestoreAccess(gfield);
77: }
78: DMSwarmDataExBegin(de);
79: DMSwarmDataExEnd(de);
80: DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
81: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
82: DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
83: for (p = 0; p < n_points_recv; p++) {
84: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
86: DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
87: }
88: if (debug) DMSwarmDataExView(de);
89: DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
90: DMSwarmDataExDestroy(de);
91: return 0;
92: }
94: PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm, DM dmcell, PetscBool remove_sent_points, PetscInt *npoints_prior_migration)
95: {
96: DM_Swarm *swarm = (DM_Swarm *)dm->data;
97: DMSwarmDataEx de;
98: PetscInt r, p, npoints, *rankval, n_points_recv;
99: PetscMPIInt rank, _rank;
100: const PetscMPIInt *neighbourranks;
101: void *point_buffer, *recv_points;
102: size_t sizeof_dmswarm_point;
103: PetscInt nneighbors;
104: PetscMPIInt mynneigh, *myneigh;
106: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
107: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
108: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
109: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
110: DMGetNeighbors(dmcell, &nneighbors, &neighbourranks);
111: DMSwarmDataExTopologyInitialize(de);
112: for (r = 0; r < nneighbors; r++) {
113: _rank = neighbourranks[r];
114: if ((_rank != rank) && (_rank > 0)) DMSwarmDataExTopologyAddNeighbour(de, _rank);
115: }
116: DMSwarmDataExTopologyFinalize(de);
117: DMSwarmDataExTopologyGetNeighbours(de, &mynneigh, &myneigh);
118: DMSwarmDataExInitializeSendCount(de);
119: for (p = 0; p < npoints; p++) {
120: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
121: for (r = 0; r < mynneigh; r++) {
122: _rank = myneigh[r];
123: DMSwarmDataExAddToSendCount(de, _rank, 1);
124: }
125: }
126: }
127: DMSwarmDataExFinalizeSendCount(de);
128: DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
129: DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
130: for (p = 0; p < npoints; p++) {
131: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
132: for (r = 0; r < mynneigh; r++) {
133: _rank = myneigh[r];
134: /* copy point into buffer */
135: DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
136: /* insert point buffer into DMSwarmDataExchanger */
137: DMSwarmDataExPackData(de, _rank, 1, point_buffer);
138: }
139: }
140: }
141: DMSwarmDataExPackFinalize(de);
142: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
143: if (remove_sent_points) {
144: DMSwarmDataField PField;
146: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &PField);
147: DMSwarmDataFieldGetEntries(PField, (void **)&rankval);
148: /* remove points which left processor */
149: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
150: for (p = 0; p < npoints; p++) {
151: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
152: /* kill point */
153: DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);
154: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL); /* you need to update npoints as the list size decreases! */
155: DMSwarmDataFieldGetEntries(PField, (void **)&rankval); /* update date point increase realloc performed */
156: p--; /* check replacement point */
157: }
158: }
159: }
160: DMSwarmDataBucketGetSizes(swarm->db, npoints_prior_migration, NULL, NULL);
161: DMSwarmDataExBegin(de);
162: DMSwarmDataExEnd(de);
163: DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
164: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
165: DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
166: for (p = 0; p < n_points_recv; p++) {
167: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
169: DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
170: }
171: DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
172: DMSwarmDataExDestroy(de);
173: return 0;
174: }
176: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm, PetscBool remove_sent_points)
177: {
178: DM_Swarm *swarm = (DM_Swarm *)dm->data;
179: PetscInt p, npoints, npointsg = 0, npoints2, npoints2g, *rankval, npoints_prior_migration;
180: PetscSF sfcell = NULL;
181: const PetscSFNode *LA_sfcell;
182: DM dmcell;
183: Vec pos;
184: PetscBool error_check = swarm->migrate_error_on_missing_point;
185: PetscMPIInt size, rank;
187: DMSwarmGetCellDM(dm, &dmcell);
190: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
191: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
193: #if 1
194: {
195: PetscInt *p_cellid;
196: PetscInt npoints_curr, range = 0;
197: PetscSFNode *sf_cells;
199: DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL);
200: PetscMalloc1(npoints_curr, &sf_cells);
202: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
203: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
204: for (p = 0; p < npoints_curr; p++) {
205: sf_cells[p].rank = 0;
206: sf_cells[p].index = p_cellid[p];
207: if (p_cellid[p] > range) range = p_cellid[p];
208: }
209: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
210: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
212: /* PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell); */
213: PetscSFCreate(PETSC_COMM_SELF, &sfcell);
214: PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
215: }
216: #endif
218: DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
219: DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
220: DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
222: if (error_check) DMSwarmGetSize(dm, &npointsg);
223: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
224: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
225: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
226: for (p = 0; p < npoints; p++) rankval[p] = LA_sfcell[p].index;
227: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
228: PetscSFDestroy(&sfcell);
230: if (size > 1) {
231: DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration);
232: } else {
233: DMSwarmDataField PField;
234: PetscInt npoints_curr;
236: /* remove points which the domain */
237: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &PField);
238: DMSwarmDataFieldGetEntries(PField, (void **)&rankval);
240: DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL);
241: for (p = 0; p < npoints_curr; p++) {
242: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
243: /* kill point */
244: DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);
245: DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL); /* you need to update npoints as the list size decreases! */
246: DMSwarmDataFieldGetEntries(PField, (void **)&rankval); /* update date point increase realloc performed */
247: p--; /* check replacement point */
248: }
249: }
250: DMSwarmGetSize(dm, &npoints_prior_migration);
251: }
253: /* locate points newly received */
254: DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL);
256: #if 0
257: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
258: PetscScalar *LA_coor;
259: PetscInt bs;
260: DMSwarmDataField PField;
262: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
263: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);
264: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
266: VecDestroy(&pos);
267: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
269: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
270: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
271: for (p=0; p<npoints2; p++) {
272: rankval[p] = LA_sfcell[p].index;
273: }
274: PetscSFDestroy(&sfcell);
275: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
277: /* remove points which left processor */
278: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
279: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
281: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
282: for (p=0; p<npoints2; p++) {
283: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
284: /* kill point */
285: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
286: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
287: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
288: p--; /* check replacement point */
289: }
290: }
291: }
292: #endif
294: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
295: PetscScalar *LA_coor;
296: PetscInt npoints_from_neighbours, bs;
297: DMSwarmDataField PField;
299: npoints_from_neighbours = npoints2 - npoints_prior_migration;
301: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor);
302: VecCreateSeqWithArray(PETSC_COMM_SELF, bs, bs * npoints_from_neighbours, (const PetscScalar *)&LA_coor[bs * npoints_prior_migration], &pos);
304: DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
306: VecDestroy(&pos);
307: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor);
309: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
310: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
311: for (p = 0; p < npoints_from_neighbours; p++) rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
312: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
313: PetscSFDestroy(&sfcell);
315: /* remove points which left processor */
316: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &PField);
317: DMSwarmDataFieldGetEntries(PField, (void **)&rankval);
319: DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL);
320: for (p = npoints_prior_migration; p < npoints2; p++) {
321: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
322: /* kill point */
323: DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);
324: DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL); /* you need to update npoints as the list size decreases! */
325: DMSwarmDataFieldGetEntries(PField, (void **)&rankval); /* update date point increase realloc performed */
326: p--; /* check replacement point */
327: }
328: }
329: }
331: {
332: PetscInt *p_cellid;
334: DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL);
335: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
336: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
337: for (p = 0; p < npoints2; p++) p_cellid[p] = rankval[p];
338: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
339: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
340: }
342: /* check for error on removed points */
343: if (error_check) {
344: DMSwarmGetSize(dm, &npoints2g);
346: }
347: return 0;
348: }
350: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm, PetscBool remove_sent_points)
351: {
352: return 0;
353: }
355: /*
356: Redundant as this assumes points can only be sent to a single rank
357: */
358: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize)
359: {
360: DM_Swarm *swarm = (DM_Swarm *)dm->data;
361: DMSwarmDataEx de;
362: PetscInt p, npoints, *rankval, n_points_recv;
363: PetscMPIInt rank, nrank, negrank;
364: void *point_buffer, *recv_points;
365: size_t sizeof_dmswarm_point;
367: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
368: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
369: *globalsize = npoints;
370: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
371: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
372: DMSwarmDataExTopologyInitialize(de);
373: for (p = 0; p < npoints; p++) {
374: negrank = rankval[p];
375: if (negrank < 0) {
376: nrank = -negrank - 1;
377: DMSwarmDataExTopologyAddNeighbour(de, nrank);
378: }
379: }
380: DMSwarmDataExTopologyFinalize(de);
381: DMSwarmDataExInitializeSendCount(de);
382: for (p = 0; p < npoints; p++) {
383: negrank = rankval[p];
384: if (negrank < 0) {
385: nrank = -negrank - 1;
386: DMSwarmDataExAddToSendCount(de, nrank, 1);
387: }
388: }
389: DMSwarmDataExFinalizeSendCount(de);
390: DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
391: DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
392: for (p = 0; p < npoints; p++) {
393: negrank = rankval[p];
394: if (negrank < 0) {
395: nrank = -negrank - 1;
396: rankval[p] = nrank;
397: /* copy point into buffer */
398: DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
399: /* insert point buffer into DMSwarmDataExchanger */
400: DMSwarmDataExPackData(de, nrank, 1, point_buffer);
401: rankval[p] = negrank;
402: }
403: }
404: DMSwarmDataExPackFinalize(de);
405: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
406: DMSwarmDataExBegin(de);
407: DMSwarmDataExEnd(de);
408: DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
409: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
410: DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
411: for (p = 0; p < n_points_recv; p++) {
412: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
414: DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
415: }
416: DMSwarmDataExView(de);
417: DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
418: DMSwarmDataExDestroy(de);
419: return 0;
420: }
422: typedef struct {
423: PetscMPIInt owner_rank;
424: PetscReal min[3], max[3];
425: } CollectBBox;
427: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm, PetscInt *globalsize)
428: {
429: DM_Swarm *swarm = (DM_Swarm *)dm->data;
430: DMSwarmDataEx de;
431: PetscInt p, pk, npoints, *rankval, n_points_recv, n_bbox_recv, dim, neighbour_cells;
432: PetscMPIInt rank, nrank;
433: void *point_buffer, *recv_points;
434: size_t sizeof_dmswarm_point, sizeof_bbox_ctx;
435: PetscBool isdmda;
436: CollectBBox *bbox, *recv_bbox;
437: const PetscMPIInt *dmneighborranks;
438: DM dmcell;
440: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
442: DMSwarmGetCellDM(dm, &dmcell);
444: isdmda = PETSC_FALSE;
445: PetscObjectTypeCompare((PetscObject)dmcell, DMDA, &isdmda);
448: DMGetDimension(dm, &dim);
449: sizeof_bbox_ctx = sizeof(CollectBBox);
450: PetscMalloc1(1, &bbox);
451: bbox->owner_rank = rank;
453: /* compute the bounding box based on the overlapping / stenctil size */
454: {
455: Vec lcoor;
457: DMGetCoordinatesLocal(dmcell, &lcoor);
458: if (dim >= 1) {
459: VecStrideMin(lcoor, 0, NULL, &bbox->min[0]);
460: VecStrideMax(lcoor, 0, NULL, &bbox->max[0]);
461: }
462: if (dim >= 2) {
463: VecStrideMin(lcoor, 1, NULL, &bbox->min[1]);
464: VecStrideMax(lcoor, 1, NULL, &bbox->max[1]);
465: }
466: if (dim == 3) {
467: VecStrideMin(lcoor, 2, NULL, &bbox->min[2]);
468: VecStrideMax(lcoor, 2, NULL, &bbox->max[2]);
469: }
470: }
471: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
472: *globalsize = npoints;
473: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
474: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
475: /* use DMDA neighbours */
476: DMDAGetNeighbors(dmcell, &dmneighborranks);
477: if (dim == 1) {
478: neighbour_cells = 3;
479: } else if (dim == 2) {
480: neighbour_cells = 9;
481: } else {
482: neighbour_cells = 27;
483: }
484: DMSwarmDataExTopologyInitialize(de);
485: for (p = 0; p < neighbour_cells; p++) {
486: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) DMSwarmDataExTopologyAddNeighbour(de, dmneighborranks[p]);
487: }
488: DMSwarmDataExTopologyFinalize(de);
489: DMSwarmDataExInitializeSendCount(de);
490: for (p = 0; p < neighbour_cells; p++) {
491: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) DMSwarmDataExAddToSendCount(de, dmneighborranks[p], 1);
492: }
493: DMSwarmDataExFinalizeSendCount(de);
494: /* send bounding boxes */
495: DMSwarmDataExPackInitialize(de, sizeof_bbox_ctx);
496: for (p = 0; p < neighbour_cells; p++) {
497: nrank = dmneighborranks[p];
498: if ((nrank >= 0) && (nrank != rank)) {
499: /* insert bbox buffer into DMSwarmDataExchanger */
500: DMSwarmDataExPackData(de, nrank, 1, bbox);
501: }
502: }
503: DMSwarmDataExPackFinalize(de);
504: /* recv bounding boxes */
505: DMSwarmDataExBegin(de);
506: DMSwarmDataExEnd(de);
507: DMSwarmDataExGetRecvData(de, &n_bbox_recv, (void **)&recv_bbox);
508: /* Wrong, should not be using PETSC_COMM_WORLD */
509: for (p = 0; p < n_bbox_recv; p++) {
510: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n", rank, recv_bbox[p].owner_rank, (double)recv_bbox[p].min[0], (double)recv_bbox[p].max[0], (double)recv_bbox[p].min[1],
511: (double)recv_bbox[p].max[1]));
512: }
513: PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout);
514: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
515: DMSwarmDataExInitializeSendCount(de);
516: for (pk = 0; pk < n_bbox_recv; pk++) {
517: PetscReal *array_x, *array_y;
519: DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x);
520: DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y);
521: for (p = 0; p < npoints; p++) {
522: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
523: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) DMSwarmDataExAddToSendCount(de, recv_bbox[pk].owner_rank, 1);
524: }
525: }
526: DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y);
527: DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x);
528: }
529: DMSwarmDataExFinalizeSendCount(de);
530: DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
531: DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
532: for (pk = 0; pk < n_bbox_recv; pk++) {
533: PetscReal *array_x, *array_y;
535: DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x);
536: DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y);
537: for (p = 0; p < npoints; p++) {
538: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
539: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
540: /* copy point into buffer */
541: DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
542: /* insert point buffer into DMSwarmDataExchanger */
543: DMSwarmDataExPackData(de, recv_bbox[pk].owner_rank, 1, point_buffer);
544: }
545: }
546: }
547: DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y);
548: DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x);
549: }
550: DMSwarmDataExPackFinalize(de);
551: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
552: DMSwarmDataExBegin(de);
553: DMSwarmDataExEnd(de);
554: DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
555: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
556: DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
557: for (p = 0; p < n_points_recv; p++) {
558: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
560: DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
561: }
562: DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
563: PetscFree(bbox);
564: DMSwarmDataExView(de);
565: DMSwarmDataExDestroy(de);
566: return 0;
567: }
569: /* General collection when no order, or neighbour information is provided */
570: /*
571: User provides context and collect() method
572: Broadcast user context
574: for each context / rank {
575: collect(swarm,context,n,list)
576: }
577: */
578: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm, PetscErrorCode (*collect)(DM, void *, PetscInt *, PetscInt **), size_t ctx_size, void *ctx, PetscInt *globalsize)
579: {
580: DM_Swarm *swarm = (DM_Swarm *)dm->data;
581: DMSwarmDataEx de;
582: PetscInt p, r, npoints, n_points_recv;
583: PetscMPIInt size, rank;
584: void *point_buffer, *recv_points;
585: void *ctxlist;
586: PetscInt *n2collect, **collectlist;
587: size_t sizeof_dmswarm_point;
589: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
590: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
591: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
592: *globalsize = npoints;
593: /* Broadcast user context */
594: PetscMalloc(ctx_size * size, &ctxlist);
595: MPI_Allgather(ctx, ctx_size, MPI_CHAR, ctxlist, ctx_size, MPI_CHAR, PetscObjectComm((PetscObject)dm));
596: PetscMalloc1(size, &n2collect);
597: PetscMalloc1(size, &collectlist);
598: for (r = 0; r < size; r++) {
599: PetscInt _n2collect;
600: PetscInt *_collectlist;
601: void *_ctx_r;
603: _n2collect = 0;
604: _collectlist = NULL;
605: if (r != rank) { /* don't collect data from yourself */
606: _ctx_r = (void *)((char *)ctxlist + r * ctx_size);
607: collect(dm, _ctx_r, &_n2collect, &_collectlist);
608: }
609: n2collect[r] = _n2collect;
610: collectlist[r] = _collectlist;
611: }
612: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
613: /* Define topology */
614: DMSwarmDataExTopologyInitialize(de);
615: for (r = 0; r < size; r++) {
616: if (n2collect[r] > 0) DMSwarmDataExTopologyAddNeighbour(de, (PetscMPIInt)r);
617: }
618: DMSwarmDataExTopologyFinalize(de);
619: /* Define send counts */
620: DMSwarmDataExInitializeSendCount(de);
621: for (r = 0; r < size; r++) {
622: if (n2collect[r] > 0) DMSwarmDataExAddToSendCount(de, r, n2collect[r]);
623: }
624: DMSwarmDataExFinalizeSendCount(de);
625: /* Pack data */
626: DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
627: DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
628: for (r = 0; r < size; r++) {
629: for (p = 0; p < n2collect[r]; p++) {
630: DMSwarmDataBucketFillPackedArray(swarm->db, collectlist[r][p], point_buffer);
631: /* insert point buffer into the data exchanger */
632: DMSwarmDataExPackData(de, r, 1, point_buffer);
633: }
634: }
635: DMSwarmDataExPackFinalize(de);
636: /* Scatter */
637: DMSwarmDataExBegin(de);
638: DMSwarmDataExEnd(de);
639: /* Collect data in DMSwarm container */
640: DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
641: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
642: DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
643: for (p = 0; p < n_points_recv; p++) {
644: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
646: DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
647: }
648: /* Release memory */
649: for (r = 0; r < size; r++) {
650: if (collectlist[r]) PetscFree(collectlist[r]);
651: }
652: PetscFree(collectlist);
653: PetscFree(n2collect);
654: PetscFree(ctxlist);
655: DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
656: DMSwarmDataExView(de);
657: DMSwarmDataExDestroy(de);
658: return 0;
659: }