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
11: It should be storing the rank information as MPIInt not Int
12: */
13: PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm, PetscBool remove_sent_points)
14: {
15: DM_Swarm *swarm = (DM_Swarm *)dm->data;
16: DMSwarmDataEx de;
17: PetscInt p, npoints, *rankval, n_points_recv;
18: PetscMPIInt rank, nrank;
19: void *point_buffer, *recv_points;
20: size_t sizeof_dmswarm_point;
21: PetscBool debug = PETSC_FALSE;
23: PetscFunctionBegin;
24: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
26: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
27: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
28: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
29: PetscCall(DMSwarmDataExTopologyInitialize(de));
30: for (p = 0; p < npoints; ++p) {
31: PetscCall(PetscMPIIntCast(rankval[p], &nrank));
32: if (nrank != rank) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
33: }
34: PetscCall(DMSwarmDataExTopologyFinalize(de));
35: PetscCall(DMSwarmDataExInitializeSendCount(de));
36: for (p = 0; p < npoints; p++) {
37: PetscCall(PetscMPIIntCast(rankval[p], &nrank));
38: if (nrank != rank) PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
39: }
40: PetscCall(DMSwarmDataExFinalizeSendCount(de));
41: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
42: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
43: for (p = 0; p < npoints; p++) {
44: PetscCall(PetscMPIIntCast(rankval[p], &nrank));
45: if (nrank != rank) {
46: /* copy point into buffer */
47: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
48: /* insert point buffer into DMSwarmDataExchanger */
49: PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
50: }
51: }
52: PetscCall(DMSwarmDataExPackFinalize(de));
53: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
55: if (remove_sent_points) {
56: DMSwarmDataField gfield;
58: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &gfield));
59: PetscCall(DMSwarmDataFieldGetAccess(gfield));
60: PetscCall(DMSwarmDataFieldGetEntries(gfield, (void **)&rankval));
62: /* remove points which left processor */
63: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
64: for (p = 0; p < npoints; p++) {
65: PetscCall(PetscMPIIntCast(rankval[p], &nrank));
66: if (nrank != rank) {
67: /* kill point */
68: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
70: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
72: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL)); /* you need to update npoints as the list size decreases! */
73: PetscCall(DMSwarmDataFieldGetAccess(gfield));
74: PetscCall(DMSwarmDataFieldGetEntries(gfield, (void **)&rankval));
75: p--; /* check replacement point */
76: }
77: }
78: PetscCall(DMSwarmDataFieldRestoreEntries(gfield, (void **)&rankval));
79: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
80: }
81: PetscCall(DMSwarmDataExBegin(de));
82: PetscCall(DMSwarmDataExEnd(de));
83: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
84: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
85: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
86: for (p = 0; p < n_points_recv; p++) {
87: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
89: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
90: }
91: if (debug) PetscCall(DMSwarmDataExView(de));
92: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
93: PetscCall(DMSwarmDataExDestroy(de));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm, DM dmcell, PetscBool remove_sent_points, PetscInt *npoints_prior_migration)
98: {
99: DM_Swarm *swarm = (DM_Swarm *)dm->data;
100: DMSwarmCellDM celldm;
101: DMSwarmDataEx de;
102: PetscInt r, p, npoints, *p_cellid, n_points_recv;
103: PetscMPIInt rank, _rank;
104: const PetscMPIInt *neighbourranks;
105: void *point_buffer, *recv_points;
106: size_t sizeof_dmswarm_point;
107: PetscInt nneighbors;
108: PetscMPIInt mynneigh, *myneigh;
109: const char *cellid;
111: PetscFunctionBegin;
112: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
113: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
114: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
115: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
116: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&p_cellid));
117: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
118: PetscCall(DMGetNeighbors(dmcell, &nneighbors, &neighbourranks));
119: PetscCall(DMSwarmDataExTopologyInitialize(de));
120: for (r = 0; r < nneighbors; r++) {
121: _rank = neighbourranks[r];
122: if ((_rank != rank) && (_rank > 0)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, _rank));
123: }
124: PetscCall(DMSwarmDataExTopologyFinalize(de));
125: PetscCall(DMSwarmDataExTopologyGetNeighbours(de, &mynneigh, &myneigh));
126: PetscCall(DMSwarmDataExInitializeSendCount(de));
127: for (p = 0; p < npoints; p++) {
128: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
129: for (r = 0; r < mynneigh; r++) {
130: _rank = myneigh[r];
131: PetscCall(DMSwarmDataExAddToSendCount(de, _rank, 1));
132: }
133: }
134: }
135: PetscCall(DMSwarmDataExFinalizeSendCount(de));
136: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
137: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
138: for (p = 0; p < npoints; p++) {
139: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
140: for (r = 0; r < mynneigh; r++) {
141: _rank = myneigh[r];
142: /* copy point into buffer */
143: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
144: /* insert point buffer into DMSwarmDataExchanger */
145: PetscCall(DMSwarmDataExPackData(de, _rank, 1, point_buffer));
146: }
147: }
148: }
149: PetscCall(DMSwarmDataExPackFinalize(de));
150: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));
151: if (remove_sent_points) {
152: DMSwarmDataField PField;
154: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, cellid, &PField));
155: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
156: /* remove points which left processor */
157: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
158: for (p = 0; p < npoints; p++) {
159: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
160: /* kill point */
161: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
162: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL)); /* you need to update npoints as the list size decreases! */
163: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point increase realloc performed */
164: p--; /* check replacement point */
165: }
166: }
167: }
168: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, npoints_prior_migration, NULL, NULL));
169: PetscCall(DMSwarmDataExBegin(de));
170: PetscCall(DMSwarmDataExEnd(de));
171: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
172: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
173: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
174: for (p = 0; p < n_points_recv; p++) {
175: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
177: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
178: }
179: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
180: PetscCall(DMSwarmDataExDestroy(de));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm, PetscBool remove_sent_points)
185: {
186: DM_Swarm *swarm = (DM_Swarm *)dm->data;
187: DMSwarmCellDM celldm;
188: PetscInt p, npoints, npointsg = 0, npoints2, npoints2g, *rankval, *p_cellid, npoints_prior_migration, Nfc;
189: PetscSF sfcell = NULL;
190: const PetscSFNode *LA_sfcell;
191: DM dmcell;
192: Vec pos;
193: PetscBool error_check = swarm->migrate_error_on_missing_point;
194: const char **coordFields;
195: PetscMPIInt size, rank;
196: const char *cellid;
198: PetscFunctionBegin;
199: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
200: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
201: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
202: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
204: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
205: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
207: #if 1
208: {
209: PetscInt npoints_curr, range = 0;
210: PetscSFNode *sf_cells;
212: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
213: PetscCall(PetscMalloc1(npoints_curr, &sf_cells));
215: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&p_cellid));
216: for (p = 0; p < npoints_curr; p++) {
217: sf_cells[p].rank = 0;
218: sf_cells[p].index = p_cellid[p];
219: if (p_cellid[p] > range) range = p_cellid[p];
220: }
221: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));
223: /* PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell)); */
224: PetscCall(PetscSFCreate(PETSC_COMM_SELF, &sfcell));
225: PetscCall(PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER));
226: }
227: #endif
229: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
230: PetscCall(DMSwarmCreateLocalVectorFromFields(dm, Nfc, coordFields, &pos));
231: PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
232: PetscCall(DMSwarmDestroyLocalVectorFromFields(dm, Nfc, coordFields, &pos));
234: if (error_check) PetscCall(DMSwarmGetSize(dm, &npointsg));
235: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
236: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&p_cellid));
237: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
238: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
240: for (p = 0; p < npoints; p++) {
241: p_cellid[p] = LA_sfcell[p].index;
242: rankval[p] = rank;
243: }
244: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
245: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));
246: PetscCall(PetscSFDestroy(&sfcell));
248: if (size > 1) {
249: PetscCall(DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration));
250: } else {
251: DMSwarmDataField PField;
252: PetscInt npoints_curr;
254: /* remove points which the domain */
255: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, cellid, &PField));
256: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
258: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
259: for (p = 0; p < npoints_curr; p++) {
260: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
261: /* kill point */
262: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
263: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL)); /* you need to update npoints as the list size decreases! */
264: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
265: p--; /* check replacement point */
266: }
267: }
268: PetscCall(DMSwarmGetLocalSize(dm, &npoints_prior_migration));
269: }
271: /* locate points newly received */
272: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
274: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
275: Vec npos;
276: IS nis;
277: PetscInt npoints_from_neighbours, bs;
278: DMSwarmDataField PField;
280: npoints_from_neighbours = npoints2 - npoints_prior_migration;
282: PetscCall(DMSwarmCreateLocalVectorFromFields(dm, Nfc, coordFields, &pos));
283: PetscCall(VecGetBlockSize(pos, &bs));
284: PetscCall(ISCreateStride(PETSC_COMM_SELF, npoints_from_neighbours * bs, npoints_prior_migration * bs, 1, &nis));
285: /*
286: Device VecGetSubVector to zero sized subvector triggers
287: debug error with mismatching memory types due to the device
288: pointer being host memtype without anything to point to in
289: Vec"Type"GetArrays(...), and we still need to pass something to
290: DMLocatePoints to avoid deadlock
291: */
292: #if defined(PETSC_HAVE_DEVICE)
293: if (npoints_from_neighbours > 0) {
294: PetscCall(VecGetSubVector(pos, nis, &npos));
295: PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
296: PetscCall(VecRestoreSubVector(pos, nis, &npos));
297: } else {
298: PetscCall(VecCreate(PETSC_COMM_SELF, &npos));
299: PetscCall(VecSetSizes(npos, 0, PETSC_DETERMINE));
300: PetscCall(VecSetBlockSize(npos, bs));
301: PetscCall(VecSetType(npos, dm->vectype));
302: PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
303: PetscCall(VecDestroy(&npos));
304: }
305: #else
306: PetscCall(VecGetSubVector(pos, nis, &npos));
307: PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
308: PetscCall(VecRestoreSubVector(pos, nis, &npos));
309: #endif
310: PetscCall(ISDestroy(&nis));
311: PetscCall(DMSwarmDestroyLocalVectorFromFields(dm, Nfc, coordFields, &pos));
313: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
314: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
315: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&p_cellid));
316: for (p = 0; p < npoints_from_neighbours; p++) {
317: rankval[npoints_prior_migration + p] = rank;
318: p_cellid[npoints_prior_migration + p] = LA_sfcell[p].index;
319: }
321: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));
322: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
324: PetscCall(PetscSFDestroy(&sfcell));
326: /* remove points which left processor */
327: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, cellid, &PField));
328: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
330: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
331: for (p = npoints_prior_migration; p < npoints2; p++) {
332: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
333: /* kill point */
334: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
335: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL)); /* you need to update npoints as the list size decreases! */
336: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
337: p--; /* check replacement point */
338: }
339: }
340: }
342: /* check for error on removed points */
343: if (error_check) {
344: PetscCall(DMSwarmGetSize(dm, &npoints2g));
345: PetscCheck(npointsg == npoints2g, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points from the DMSwarm must remain constant during migration (initial %" PetscInt_FMT " - final %" PetscInt_FMT ")", npointsg, npoints2g);
346: }
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm, PetscBool remove_sent_points)
351: {
352: PetscFunctionBegin;
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*
357: Redundant as this assumes points can only be sent to a single rank
358: */
359: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize)
360: {
361: DM_Swarm *swarm = (DM_Swarm *)dm->data;
362: DMSwarmDataEx de;
363: PetscInt p, npoints, *rankval, n_points_recv;
364: PetscMPIInt rank, nrank, negrank;
365: void *point_buffer, *recv_points;
366: size_t sizeof_dmswarm_point;
368: PetscFunctionBegin;
369: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
370: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
371: *globalsize = npoints;
372: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
373: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
374: PetscCall(DMSwarmDataExTopologyInitialize(de));
375: for (p = 0; p < npoints; p++) {
376: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
377: if (negrank < 0) {
378: nrank = -negrank - 1;
379: PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
380: }
381: }
382: PetscCall(DMSwarmDataExTopologyFinalize(de));
383: PetscCall(DMSwarmDataExInitializeSendCount(de));
384: for (p = 0; p < npoints; p++) {
385: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
386: if (negrank < 0) {
387: nrank = -negrank - 1;
388: PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
389: }
390: }
391: PetscCall(DMSwarmDataExFinalizeSendCount(de));
392: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
393: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
394: for (p = 0; p < npoints; p++) {
395: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
396: if (negrank < 0) {
397: nrank = -negrank - 1;
398: rankval[p] = nrank;
399: /* copy point into buffer */
400: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
401: /* insert point buffer into DMSwarmDataExchanger */
402: PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
403: rankval[p] = negrank;
404: }
405: }
406: PetscCall(DMSwarmDataExPackFinalize(de));
407: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
408: PetscCall(DMSwarmDataExBegin(de));
409: PetscCall(DMSwarmDataExEnd(de));
410: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
411: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
412: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
413: for (p = 0; p < n_points_recv; p++) {
414: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
416: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
417: }
418: PetscCall(DMSwarmDataExView(de));
419: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
420: PetscCall(DMSwarmDataExDestroy(de));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: typedef struct {
425: PetscMPIInt owner_rank;
426: PetscReal min[3], max[3];
427: } CollectBBox;
429: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm, PetscInt *globalsize)
430: {
431: DM_Swarm *swarm = (DM_Swarm *)dm->data;
432: DMSwarmDataEx de;
433: PetscInt p, pk, npoints, *rankval, n_points_recv, n_bbox_recv, dim, neighbour_cells;
434: PetscMPIInt rank, nrank;
435: void *point_buffer, *recv_points;
436: size_t sizeof_dmswarm_point, sizeof_bbox_ctx;
437: PetscBool isdmda;
438: CollectBBox *bbox, *recv_bbox;
439: const PetscMPIInt *dmneighborranks;
440: DM dmcell;
442: PetscFunctionBegin;
443: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
445: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
446: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
447: isdmda = PETSC_FALSE;
448: PetscCall(PetscObjectTypeCompare((PetscObject)dmcell, DMDA, &isdmda));
449: PetscCheck(isdmda, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only DMDA support for CollectBoundingBox");
451: PetscCall(DMGetDimension(dm, &dim));
452: sizeof_bbox_ctx = sizeof(CollectBBox);
453: PetscCall(PetscMalloc1(1, &bbox));
454: bbox->owner_rank = rank;
456: /* compute the bounding box based on the overlapping / stenctil size */
457: {
458: Vec lcoor;
460: PetscCall(DMGetCoordinatesLocal(dmcell, &lcoor));
461: if (dim >= 1) {
462: PetscCall(VecStrideMin(lcoor, 0, NULL, &bbox->min[0]));
463: PetscCall(VecStrideMax(lcoor, 0, NULL, &bbox->max[0]));
464: }
465: if (dim >= 2) {
466: PetscCall(VecStrideMin(lcoor, 1, NULL, &bbox->min[1]));
467: PetscCall(VecStrideMax(lcoor, 1, NULL, &bbox->max[1]));
468: }
469: if (dim == 3) {
470: PetscCall(VecStrideMin(lcoor, 2, NULL, &bbox->min[2]));
471: PetscCall(VecStrideMax(lcoor, 2, NULL, &bbox->max[2]));
472: }
473: }
474: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
475: *globalsize = npoints;
476: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
477: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
478: /* use DMDA neighbours */
479: PetscCall(DMDAGetNeighbors(dmcell, &dmneighborranks));
480: if (dim == 1) {
481: neighbour_cells = 3;
482: } else if (dim == 2) {
483: neighbour_cells = 9;
484: } else {
485: neighbour_cells = 27;
486: }
487: PetscCall(DMSwarmDataExTopologyInitialize(de));
488: for (p = 0; p < neighbour_cells; p++) {
489: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, dmneighborranks[p]));
490: }
491: PetscCall(DMSwarmDataExTopologyFinalize(de));
492: PetscCall(DMSwarmDataExInitializeSendCount(de));
493: for (p = 0; p < neighbour_cells; p++) {
494: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExAddToSendCount(de, dmneighborranks[p], 1));
495: }
496: PetscCall(DMSwarmDataExFinalizeSendCount(de));
497: /* send bounding boxes */
498: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_bbox_ctx));
499: for (p = 0; p < neighbour_cells; p++) {
500: nrank = dmneighborranks[p];
501: if ((nrank >= 0) && (nrank != rank)) {
502: /* insert bbox buffer into DMSwarmDataExchanger */
503: PetscCall(DMSwarmDataExPackData(de, nrank, 1, bbox));
504: }
505: }
506: PetscCall(DMSwarmDataExPackFinalize(de));
507: /* recv bounding boxes */
508: PetscCall(DMSwarmDataExBegin(de));
509: PetscCall(DMSwarmDataExEnd(de));
510: PetscCall(DMSwarmDataExGetRecvData(de, &n_bbox_recv, (void **)&recv_bbox));
511: /* Wrong, should not be using PETSC_COMM_WORLD */
512: for (p = 0; p < n_bbox_recv; p++) {
513: 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],
514: (double)recv_bbox[p].max[1]));
515: }
516: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
517: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
518: PetscCall(DMSwarmDataExInitializeSendCount(de));
519: for (pk = 0; pk < n_bbox_recv; pk++) {
520: PetscReal *array_x, *array_y;
522: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
523: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
524: for (p = 0; p < npoints; p++) {
525: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
526: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) PetscCall(DMSwarmDataExAddToSendCount(de, recv_bbox[pk].owner_rank, 1));
527: }
528: }
529: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
530: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
531: }
532: PetscCall(DMSwarmDataExFinalizeSendCount(de));
533: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
534: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
535: for (pk = 0; pk < n_bbox_recv; pk++) {
536: PetscReal *array_x, *array_y;
538: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
539: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
540: for (p = 0; p < npoints; p++) {
541: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
542: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
543: /* copy point into buffer */
544: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
545: /* insert point buffer into DMSwarmDataExchanger */
546: PetscCall(DMSwarmDataExPackData(de, recv_bbox[pk].owner_rank, 1, point_buffer));
547: }
548: }
549: }
550: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
551: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
552: }
553: PetscCall(DMSwarmDataExPackFinalize(de));
554: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
555: PetscCall(DMSwarmDataExBegin(de));
556: PetscCall(DMSwarmDataExEnd(de));
557: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
558: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
559: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
560: for (p = 0; p < n_points_recv; p++) {
561: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
563: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
564: }
565: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
566: PetscCall(PetscFree(bbox));
567: PetscCall(DMSwarmDataExView(de));
568: PetscCall(DMSwarmDataExDestroy(de));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /* General collection when no order, or neighbour information is provided */
573: /*
574: User provides context and collect() method
575: Broadcast user context
577: for each context / rank {
578: collect(swarm,context,n,list)
579: }
580: */
581: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm, PetscErrorCode (*collect)(DM, void *, PetscInt *, PetscInt **), size_t ctx_size, void *ctx, PetscInt *globalsize)
582: {
583: DM_Swarm *swarm = (DM_Swarm *)dm->data;
584: DMSwarmDataEx de;
585: PetscInt p, npoints, n_points_recv;
586: PetscMPIInt size, rank, len;
587: void *point_buffer, *recv_points;
588: void *ctxlist;
589: PetscInt *n2collect, **collectlist;
590: size_t sizeof_dmswarm_point;
592: PetscFunctionBegin;
593: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
594: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
595: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
596: *globalsize = npoints;
597: /* Broadcast user context */
598: PetscCall(PetscMPIIntCast(ctx_size, &len));
599: PetscCall(PetscMalloc(ctx_size * size, &ctxlist));
600: PetscCallMPI(MPI_Allgather(ctx, len, MPI_CHAR, ctxlist, len, MPI_CHAR, PetscObjectComm((PetscObject)dm)));
601: PetscCall(PetscMalloc1(size, &n2collect));
602: PetscCall(PetscMalloc1(size, &collectlist));
603: for (PetscMPIInt r = 0; r < size; r++) {
604: PetscInt _n2collect;
605: PetscInt *_collectlist;
606: void *_ctx_r;
608: _n2collect = 0;
609: _collectlist = NULL;
610: if (r != rank) { /* don't collect data from yourself */
611: _ctx_r = (void *)((char *)ctxlist + r * ctx_size);
612: PetscCall(collect(dm, _ctx_r, &_n2collect, &_collectlist));
613: }
614: n2collect[r] = _n2collect;
615: collectlist[r] = _collectlist;
616: }
617: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
618: /* Define topology */
619: PetscCall(DMSwarmDataExTopologyInitialize(de));
620: for (PetscMPIInt r = 0; r < size; r++) {
621: if (n2collect[r] > 0) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, r));
622: }
623: PetscCall(DMSwarmDataExTopologyFinalize(de));
624: /* Define send counts */
625: PetscCall(DMSwarmDataExInitializeSendCount(de));
626: for (PetscMPIInt r = 0; r < size; r++) {
627: if (n2collect[r] > 0) PetscCall(DMSwarmDataExAddToSendCount(de, r, n2collect[r]));
628: }
629: PetscCall(DMSwarmDataExFinalizeSendCount(de));
630: /* Pack data */
631: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
632: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
633: for (PetscMPIInt r = 0; r < size; r++) {
634: for (p = 0; p < n2collect[r]; p++) {
635: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, collectlist[r][p], point_buffer));
636: /* insert point buffer into the data exchanger */
637: PetscCall(DMSwarmDataExPackData(de, r, 1, point_buffer));
638: }
639: }
640: PetscCall(DMSwarmDataExPackFinalize(de));
641: /* Scatter */
642: PetscCall(DMSwarmDataExBegin(de));
643: PetscCall(DMSwarmDataExEnd(de));
644: /* Collect data in DMSwarm container */
645: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
646: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
647: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
648: for (p = 0; p < n_points_recv; p++) {
649: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
651: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
652: }
653: /* Release memory */
654: for (PetscMPIInt r = 0; r < size; r++) {
655: if (collectlist[r]) PetscCall(PetscFree(collectlist[r]));
656: }
657: PetscCall(PetscFree(collectlist));
658: PetscCall(PetscFree(n2collect));
659: PetscCall(PetscFree(ctxlist));
660: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
661: PetscCall(DMSwarmDataExView(de));
662: PetscCall(DMSwarmDataExDestroy(de));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: /*@
667: DMSwarmGetMigrateType - Get the style of point migration
669: Logically Collective
671: Input Parameter:
672: . dm - the `DMSWARM`
674: Output Parameter:
675: . mtype - The migration type, see `DMSwarmMigrateType`
677: Level: intermediate
679: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmMigrate()`
680: @*/
681: PetscErrorCode DMSwarmGetMigrateType(DM dm, DMSwarmMigrateType *mtype)
682: {
683: PetscFunctionBegin;
685: PetscAssertPointer(mtype, 2);
686: *mtype = ((DM_Swarm *)dm->data)->migrate_type;
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: /*@
691: DMSwarmSetMigrateType - Set the style of point migration
693: Logically Collective
695: Input Parameters:
696: + dm - the `DMSWARM`
697: - mtype - The migration type, see `DMSwarmMigrateType`
699: Level: intermediate
701: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmGetMigrateType()`, `DMSwarmMigrate()`
702: @*/
703: PetscErrorCode DMSwarmSetMigrateType(DM dm, DMSwarmMigrateType mtype)
704: {
705: PetscFunctionBegin;
708: ((DM_Swarm *)dm->data)->migrate_type = mtype;
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }