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: if (remove_sent_points) {
260: for (p = 0; p < npoints_curr; p++) {
261: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
262: /* kill point */
263: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
264: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL)); /* you need to update npoints as the list size decreases! */
265: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
266: p--; /* check replacement point */
267: }
268: }
269: }
270: PetscCall(DMSwarmGetLocalSize(dm, &npoints_prior_migration));
271: }
273: /* locate points newly received */
274: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
276: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
277: Vec npos;
278: IS nis;
279: PetscInt npoints_from_neighbours, bs;
280: DMSwarmDataField PField;
282: npoints_from_neighbours = npoints2 - npoints_prior_migration;
284: PetscCall(DMSwarmCreateLocalVectorFromFields(dm, Nfc, coordFields, &pos));
285: PetscCall(VecGetBlockSize(pos, &bs));
286: PetscCall(ISCreateStride(PETSC_COMM_SELF, npoints_from_neighbours * bs, npoints_prior_migration * bs, 1, &nis));
287: /*
288: Device VecGetSubVector to zero sized subvector triggers
289: debug error with mismatching memory types due to the device
290: pointer being host memtype without anything to point to in
291: Vec"Type"GetArrays(...), and we still need to pass something to
292: DMLocatePoints to avoid deadlock
293: */
294: #if defined(PETSC_HAVE_DEVICE)
295: if (npoints_from_neighbours > 0) {
296: PetscCall(VecGetSubVector(pos, nis, &npos));
297: PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
298: PetscCall(VecRestoreSubVector(pos, nis, &npos));
299: } else {
300: PetscCall(VecCreate(PETSC_COMM_SELF, &npos));
301: PetscCall(VecSetSizes(npos, 0, PETSC_DETERMINE));
302: PetscCall(VecSetBlockSize(npos, bs));
303: PetscCall(VecSetType(npos, dm->vectype));
304: PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
305: PetscCall(VecDestroy(&npos));
306: }
307: #else
308: PetscCall(VecGetSubVector(pos, nis, &npos));
309: PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
310: PetscCall(VecRestoreSubVector(pos, nis, &npos));
311: #endif
312: PetscCall(ISDestroy(&nis));
313: PetscCall(DMSwarmDestroyLocalVectorFromFields(dm, Nfc, coordFields, &pos));
315: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
316: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
317: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&p_cellid));
318: for (p = 0; p < npoints_from_neighbours; p++) {
319: rankval[npoints_prior_migration + p] = rank;
320: p_cellid[npoints_prior_migration + p] = LA_sfcell[p].index;
321: }
323: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));
324: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
326: PetscCall(PetscSFDestroy(&sfcell));
328: /* remove points which left processor */
329: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, cellid, &PField));
330: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
332: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
333: if (remove_sent_points) {
334: for (p = npoints_prior_migration; p < npoints2; p++) {
335: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
336: /* kill point */
337: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
338: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL)); /* you need to update npoints as the list size decreases! */
339: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
340: p--; /* check replacement point */
341: }
342: }
343: }
344: }
346: /* check for error on removed points */
347: if (error_check) {
348: PetscCall(DMSwarmGetSize(dm, &npoints2g));
349: 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);
350: }
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm, PetscBool remove_sent_points)
355: {
356: PetscFunctionBegin;
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: /*
361: Redundant as this assumes points can only be sent to a single rank
362: */
363: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize)
364: {
365: DM_Swarm *swarm = (DM_Swarm *)dm->data;
366: DMSwarmDataEx de;
367: PetscInt p, npoints, *rankval, n_points_recv;
368: PetscMPIInt rank, nrank, negrank;
369: void *point_buffer, *recv_points;
370: size_t sizeof_dmswarm_point;
372: PetscFunctionBegin;
373: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
374: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
375: *globalsize = npoints;
376: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
377: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
378: PetscCall(DMSwarmDataExTopologyInitialize(de));
379: for (p = 0; p < npoints; p++) {
380: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
381: if (negrank < 0) {
382: nrank = -negrank - 1;
383: PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
384: }
385: }
386: PetscCall(DMSwarmDataExTopologyFinalize(de));
387: PetscCall(DMSwarmDataExInitializeSendCount(de));
388: for (p = 0; p < npoints; p++) {
389: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
390: if (negrank < 0) {
391: nrank = -negrank - 1;
392: PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
393: }
394: }
395: PetscCall(DMSwarmDataExFinalizeSendCount(de));
396: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
397: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
398: for (p = 0; p < npoints; p++) {
399: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
400: if (negrank < 0) {
401: nrank = -negrank - 1;
402: rankval[p] = nrank;
403: /* copy point into buffer */
404: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
405: /* insert point buffer into DMSwarmDataExchanger */
406: PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
407: rankval[p] = negrank;
408: }
409: }
410: PetscCall(DMSwarmDataExPackFinalize(de));
411: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
412: PetscCall(DMSwarmDataExBegin(de));
413: PetscCall(DMSwarmDataExEnd(de));
414: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
415: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
416: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
417: for (p = 0; p < n_points_recv; p++) {
418: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
420: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
421: }
422: PetscCall(DMSwarmDataExView(de));
423: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
424: PetscCall(DMSwarmDataExDestroy(de));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: typedef struct {
429: PetscMPIInt owner_rank;
430: PetscReal min[3], max[3];
431: } CollectBBox;
433: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm, PetscInt *globalsize)
434: {
435: DM_Swarm *swarm = (DM_Swarm *)dm->data;
436: DMSwarmDataEx de;
437: PetscInt p, pk, npoints, *rankval, n_points_recv, n_bbox_recv, dim, neighbour_cells;
438: PetscMPIInt rank, nrank;
439: void *point_buffer, *recv_points;
440: size_t sizeof_dmswarm_point, sizeof_bbox_ctx;
441: PetscBool isdmda;
442: CollectBBox *bbox, *recv_bbox;
443: const PetscMPIInt *dmneighborranks;
444: DM dmcell;
446: PetscFunctionBegin;
447: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
449: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
450: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
451: isdmda = PETSC_FALSE;
452: PetscCall(PetscObjectTypeCompare((PetscObject)dmcell, DMDA, &isdmda));
453: PetscCheck(isdmda, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only DMDA support for CollectBoundingBox");
455: PetscCall(DMGetDimension(dm, &dim));
456: sizeof_bbox_ctx = sizeof(CollectBBox);
457: PetscCall(PetscMalloc1(1, &bbox));
458: bbox->owner_rank = rank;
460: /* compute the bounding box based on the overlapping / stenctil size */
461: {
462: Vec lcoor;
464: PetscCall(DMGetCoordinatesLocal(dmcell, &lcoor));
465: if (dim >= 1) {
466: PetscCall(VecStrideMin(lcoor, 0, NULL, &bbox->min[0]));
467: PetscCall(VecStrideMax(lcoor, 0, NULL, &bbox->max[0]));
468: }
469: if (dim >= 2) {
470: PetscCall(VecStrideMin(lcoor, 1, NULL, &bbox->min[1]));
471: PetscCall(VecStrideMax(lcoor, 1, NULL, &bbox->max[1]));
472: }
473: if (dim == 3) {
474: PetscCall(VecStrideMin(lcoor, 2, NULL, &bbox->min[2]));
475: PetscCall(VecStrideMax(lcoor, 2, NULL, &bbox->max[2]));
476: }
477: }
478: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
479: *globalsize = npoints;
480: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
481: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
482: /* use DMDA neighbours */
483: PetscCall(DMDAGetNeighbors(dmcell, &dmneighborranks));
484: if (dim == 1) {
485: neighbour_cells = 3;
486: } else if (dim == 2) {
487: neighbour_cells = 9;
488: } else {
489: neighbour_cells = 27;
490: }
491: PetscCall(DMSwarmDataExTopologyInitialize(de));
492: for (p = 0; p < neighbour_cells; p++) {
493: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, dmneighborranks[p]));
494: }
495: PetscCall(DMSwarmDataExTopologyFinalize(de));
496: PetscCall(DMSwarmDataExInitializeSendCount(de));
497: for (p = 0; p < neighbour_cells; p++) {
498: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExAddToSendCount(de, dmneighborranks[p], 1));
499: }
500: PetscCall(DMSwarmDataExFinalizeSendCount(de));
501: /* send bounding boxes */
502: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_bbox_ctx));
503: for (p = 0; p < neighbour_cells; p++) {
504: nrank = dmneighborranks[p];
505: if ((nrank >= 0) && (nrank != rank)) {
506: /* insert bbox buffer into DMSwarmDataExchanger */
507: PetscCall(DMSwarmDataExPackData(de, nrank, 1, bbox));
508: }
509: }
510: PetscCall(DMSwarmDataExPackFinalize(de));
511: /* recv bounding boxes */
512: PetscCall(DMSwarmDataExBegin(de));
513: PetscCall(DMSwarmDataExEnd(de));
514: PetscCall(DMSwarmDataExGetRecvData(de, &n_bbox_recv, (void **)&recv_bbox));
515: /* Wrong, should not be using PETSC_COMM_WORLD */
516: for (p = 0; p < n_bbox_recv; p++) {
517: 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],
518: (double)recv_bbox[p].max[1]));
519: }
520: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
521: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
522: PetscCall(DMSwarmDataExInitializeSendCount(de));
523: for (pk = 0; pk < n_bbox_recv; pk++) {
524: PetscReal *array_x, *array_y;
526: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
527: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
528: for (p = 0; p < npoints; p++) {
529: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
530: 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));
531: }
532: }
533: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
534: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
535: }
536: PetscCall(DMSwarmDataExFinalizeSendCount(de));
537: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
538: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
539: for (pk = 0; pk < n_bbox_recv; pk++) {
540: PetscReal *array_x, *array_y;
542: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
543: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
544: for (p = 0; p < npoints; p++) {
545: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
546: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
547: /* copy point into buffer */
548: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
549: /* insert point buffer into DMSwarmDataExchanger */
550: PetscCall(DMSwarmDataExPackData(de, recv_bbox[pk].owner_rank, 1, point_buffer));
551: }
552: }
553: }
554: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
555: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
556: }
557: PetscCall(DMSwarmDataExPackFinalize(de));
558: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
559: PetscCall(DMSwarmDataExBegin(de));
560: PetscCall(DMSwarmDataExEnd(de));
561: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
562: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
563: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
564: for (p = 0; p < n_points_recv; p++) {
565: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
567: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
568: }
569: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
570: PetscCall(PetscFree(bbox));
571: PetscCall(DMSwarmDataExView(de));
572: PetscCall(DMSwarmDataExDestroy(de));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /* General collection when no order, or neighbour information is provided */
577: /*
578: User provides context and collect() method
579: Broadcast user context
581: for each context / rank {
582: collect(swarm,context,n,list)
583: }
584: */
585: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm, PetscErrorCode (*collect)(DM, void *, PetscInt *, PetscInt **), size_t ctx_size, void *ctx, PetscInt *globalsize)
586: {
587: DM_Swarm *swarm = (DM_Swarm *)dm->data;
588: DMSwarmDataEx de;
589: PetscInt p, npoints, n_points_recv;
590: PetscMPIInt size, rank, len;
591: void *point_buffer, *recv_points;
592: void *ctxlist;
593: PetscInt *n2collect, **collectlist;
594: size_t sizeof_dmswarm_point;
596: PetscFunctionBegin;
597: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
598: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
599: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
600: *globalsize = npoints;
601: /* Broadcast user context */
602: PetscCall(PetscMPIIntCast(ctx_size, &len));
603: PetscCall(PetscMalloc(ctx_size * size, &ctxlist));
604: PetscCallMPI(MPI_Allgather(ctx, len, MPI_CHAR, ctxlist, len, MPI_CHAR, PetscObjectComm((PetscObject)dm)));
605: PetscCall(PetscMalloc1(size, &n2collect));
606: PetscCall(PetscMalloc1(size, &collectlist));
607: for (PetscMPIInt r = 0; r < size; r++) {
608: PetscInt _n2collect;
609: PetscInt *_collectlist;
610: void *_ctx_r;
612: _n2collect = 0;
613: _collectlist = NULL;
614: if (r != rank) { /* don't collect data from yourself */
615: _ctx_r = (void *)((char *)ctxlist + r * ctx_size);
616: PetscCall(collect(dm, _ctx_r, &_n2collect, &_collectlist));
617: }
618: n2collect[r] = _n2collect;
619: collectlist[r] = _collectlist;
620: }
621: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
622: /* Define topology */
623: PetscCall(DMSwarmDataExTopologyInitialize(de));
624: for (PetscMPIInt r = 0; r < size; r++) {
625: if (n2collect[r] > 0) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, r));
626: }
627: PetscCall(DMSwarmDataExTopologyFinalize(de));
628: /* Define send counts */
629: PetscCall(DMSwarmDataExInitializeSendCount(de));
630: for (PetscMPIInt r = 0; r < size; r++) {
631: if (n2collect[r] > 0) PetscCall(DMSwarmDataExAddToSendCount(de, r, n2collect[r]));
632: }
633: PetscCall(DMSwarmDataExFinalizeSendCount(de));
634: /* Pack data */
635: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
636: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
637: for (PetscMPIInt r = 0; r < size; r++) {
638: for (p = 0; p < n2collect[r]; p++) {
639: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, collectlist[r][p], point_buffer));
640: /* insert point buffer into the data exchanger */
641: PetscCall(DMSwarmDataExPackData(de, r, 1, point_buffer));
642: }
643: }
644: PetscCall(DMSwarmDataExPackFinalize(de));
645: /* Scatter */
646: PetscCall(DMSwarmDataExBegin(de));
647: PetscCall(DMSwarmDataExEnd(de));
648: /* Collect data in DMSwarm container */
649: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
650: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
651: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
652: for (p = 0; p < n_points_recv; p++) {
653: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
655: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
656: }
657: /* Release memory */
658: for (PetscMPIInt r = 0; r < size; r++) {
659: if (collectlist[r]) PetscCall(PetscFree(collectlist[r]));
660: }
661: PetscCall(PetscFree(collectlist));
662: PetscCall(PetscFree(n2collect));
663: PetscCall(PetscFree(ctxlist));
664: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
665: PetscCall(DMSwarmDataExView(de));
666: PetscCall(DMSwarmDataExDestroy(de));
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: /*@
671: DMSwarmGetMigrateType - Get the style of point migration
673: Logically Collective
675: Input Parameter:
676: . dm - the `DMSWARM`
678: Output Parameter:
679: . mtype - The migration type, see `DMSwarmMigrateType`
681: Level: intermediate
683: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmMigrate()`
684: @*/
685: PetscErrorCode DMSwarmGetMigrateType(DM dm, DMSwarmMigrateType *mtype)
686: {
687: PetscFunctionBegin;
689: PetscAssertPointer(mtype, 2);
690: *mtype = ((DM_Swarm *)dm->data)->migrate_type;
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: /*@
695: DMSwarmSetMigrateType - Set the style of point migration
697: Logically Collective
699: Input Parameters:
700: + dm - the `DMSWARM`
701: - mtype - The migration type, see `DMSwarmMigrateType`
703: Level: intermediate
705: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmGetMigrateType()`, `DMSwarmMigrate()`
706: @*/
707: PetscErrorCode DMSwarmSetMigrateType(DM dm, DMSwarmMigrateType mtype)
708: {
709: PetscFunctionBegin;
712: ((DM_Swarm *)dm->data)->migrate_type = mtype;
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }