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: PetscFunctionBegin;
22: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
24: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
25: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
26: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
27: PetscCall(DMSwarmDataExTopologyInitialize(de));
28: for (p = 0; p < npoints; ++p) {
29: nrank = rankval[p];
30: if (nrank != rank) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
31: }
32: PetscCall(DMSwarmDataExTopologyFinalize(de));
33: PetscCall(DMSwarmDataExInitializeSendCount(de));
34: for (p = 0; p < npoints; p++) {
35: nrank = rankval[p];
36: if (nrank != rank) PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
37: }
38: PetscCall(DMSwarmDataExFinalizeSendCount(de));
39: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
40: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
41: for (p = 0; p < npoints; p++) {
42: nrank = rankval[p];
43: if (nrank != rank) {
44: /* copy point into buffer */
45: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
46: /* insert point buffer into DMSwarmDataExchanger */
47: PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
48: }
49: }
50: PetscCall(DMSwarmDataExPackFinalize(de));
51: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
53: if (remove_sent_points) {
54: DMSwarmDataField gfield;
56: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &gfield));
57: PetscCall(DMSwarmDataFieldGetAccess(gfield));
58: PetscCall(DMSwarmDataFieldGetEntries(gfield, (void **)&rankval));
60: /* remove points which left processor */
61: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
62: for (p = 0; p < npoints; p++) {
63: nrank = rankval[p];
64: if (nrank != rank) {
65: /* kill point */
66: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
68: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
70: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL)); /* you need to update npoints as the list size decreases! */
71: PetscCall(DMSwarmDataFieldGetAccess(gfield));
72: PetscCall(DMSwarmDataFieldGetEntries(gfield, (void **)&rankval));
73: p--; /* check replacement point */
74: }
75: }
76: PetscCall(DMSwarmDataFieldRestoreEntries(gfield, (void **)&rankval));
77: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
78: }
79: PetscCall(DMSwarmDataExBegin(de));
80: PetscCall(DMSwarmDataExEnd(de));
81: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
82: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
83: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
84: for (p = 0; p < n_points_recv; p++) {
85: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
87: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
88: }
89: if (debug) PetscCall(DMSwarmDataExView(de));
90: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
91: PetscCall(DMSwarmDataExDestroy(de));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: static PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm, DM dmcell, PetscBool remove_sent_points, PetscInt *npoints_prior_migration)
96: {
97: DM_Swarm *swarm = (DM_Swarm *)dm->data;
98: DMSwarmDataEx de;
99: PetscInt r, p, npoints, *p_cellid, n_points_recv;
100: PetscMPIInt rank, _rank;
101: const PetscMPIInt *neighbourranks;
102: void *point_buffer, *recv_points;
103: size_t sizeof_dmswarm_point;
104: PetscInt nneighbors;
105: PetscMPIInt mynneigh, *myneigh;
107: PetscFunctionBegin;
108: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
109: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
110: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
111: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
112: PetscCall(DMGetNeighbors(dmcell, &nneighbors, &neighbourranks));
113: PetscCall(DMSwarmDataExTopologyInitialize(de));
114: for (r = 0; r < nneighbors; r++) {
115: _rank = neighbourranks[r];
116: if ((_rank != rank) && (_rank > 0)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, _rank));
117: }
118: PetscCall(DMSwarmDataExTopologyFinalize(de));
119: PetscCall(DMSwarmDataExTopologyGetNeighbours(de, &mynneigh, &myneigh));
120: PetscCall(DMSwarmDataExInitializeSendCount(de));
121: for (p = 0; p < npoints; p++) {
122: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
123: for (r = 0; r < mynneigh; r++) {
124: _rank = myneigh[r];
125: PetscCall(DMSwarmDataExAddToSendCount(de, _rank, 1));
126: }
127: }
128: }
129: PetscCall(DMSwarmDataExFinalizeSendCount(de));
130: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
131: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
132: for (p = 0; p < npoints; p++) {
133: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
134: for (r = 0; r < mynneigh; r++) {
135: _rank = myneigh[r];
136: /* copy point into buffer */
137: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
138: /* insert point buffer into DMSwarmDataExchanger */
139: PetscCall(DMSwarmDataExPackData(de, _rank, 1, point_buffer));
140: }
141: }
142: }
143: PetscCall(DMSwarmDataExPackFinalize(de));
144: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
145: if (remove_sent_points) {
146: DMSwarmDataField PField;
148: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
149: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
150: /* remove points which left processor */
151: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
152: for (p = 0; p < npoints; p++) {
153: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
154: /* kill point */
155: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
156: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL)); /* you need to update npoints as the list size decreases! */
157: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point increase realloc performed */
158: p--; /* check replacement point */
159: }
160: }
161: }
162: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, npoints_prior_migration, NULL, NULL));
163: PetscCall(DMSwarmDataExBegin(de));
164: PetscCall(DMSwarmDataExEnd(de));
165: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
166: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
167: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
168: for (p = 0; p < n_points_recv; p++) {
169: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
171: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
172: }
173: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
174: PetscCall(DMSwarmDataExDestroy(de));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm, PetscBool remove_sent_points)
179: {
180: DM_Swarm *swarm = (DM_Swarm *)dm->data;
181: PetscInt p, npoints, npointsg = 0, npoints2, npoints2g, *rankval, *p_cellid, npoints_prior_migration;
182: PetscSF sfcell = NULL;
183: const PetscSFNode *LA_sfcell;
184: DM dmcell;
185: Vec pos;
186: PetscBool error_check = swarm->migrate_error_on_missing_point;
187: PetscMPIInt size, rank;
189: PetscFunctionBegin;
190: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
191: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
193: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
194: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
196: #if 1
197: {
198: PetscInt npoints_curr, range = 0;
199: PetscSFNode *sf_cells;
201: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
202: PetscCall(PetscMalloc1(npoints_curr, &sf_cells));
204: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
205: for (p = 0; p < npoints_curr; p++) {
206: sf_cells[p].rank = 0;
207: sf_cells[p].index = p_cellid[p];
208: if (p_cellid[p] > range) range = p_cellid[p];
209: }
210: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
212: /* PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell)); */
213: PetscCall(PetscSFCreate(PETSC_COMM_SELF, &sfcell));
214: PetscCall(PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER));
215: }
216: #endif
218: PetscCall(DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos));
219: PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
220: PetscCall(DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos));
222: if (error_check) PetscCall(DMSwarmGetSize(dm, &npointsg));
223: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
224: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
225: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
226: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
228: for (p = 0; p < npoints; p++) {
229: p_cellid[p] = LA_sfcell[p].index;
230: rankval[p] = rank;
231: }
232: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
233: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
234: PetscCall(PetscSFDestroy(&sfcell));
236: if (size > 1) {
237: PetscCall(DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration));
238: } else {
239: DMSwarmDataField PField;
240: PetscInt npoints_curr;
242: /* remove points which the domain */
243: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
244: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
246: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
247: for (p = 0; p < npoints_curr; p++) {
248: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
249: /* kill point */
250: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
251: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL)); /* you need to update npoints as the list size decreases! */
252: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
253: p--; /* check replacement point */
254: }
255: }
256: PetscCall(DMSwarmGetLocalSize(dm, &npoints_prior_migration));
257: }
259: /* locate points newly received */
260: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
262: #if 0
263: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
264: PetscScalar *LA_coor;
265: PetscInt bs;
266: DMSwarmDataField PField;
268: PetscCall(DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor));
269: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos));
270: PetscCall(DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell));
272: PetscCall(VecDestroy(&pos));
273: PetscCall(DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor));
275: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
276: PetscCall(DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));
277: for (p=0; p<npoints2; p++) {
278: rankval[p] = LA_sfcell[p].index;
279: }
280: PetscCall(PetscSFDestroy(&sfcell));
281: PetscCall(DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));
283: /* remove points which left processor */
284: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField));
285: PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval));
287: PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL));
288: for (p=0; p<npoints2; p++) {
289: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
290: /* kill point */
291: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db,p));
292: PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL)); /* you need to update npoints as the list size decreases! */
293: PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval)); /* update date point increase realloc performed */
294: p--; /* check replacement point */
295: }
296: }
297: }
298: #endif
300: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
301: PetscScalar *LA_coor;
302: PetscInt npoints_from_neighbours, bs;
303: DMSwarmDataField PField;
305: npoints_from_neighbours = npoints2 - npoints_prior_migration;
307: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor));
308: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, bs, bs * npoints_from_neighbours, (const PetscScalar *)&LA_coor[bs * npoints_prior_migration], &pos));
310: PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
312: PetscCall(VecDestroy(&pos));
313: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor));
315: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
316: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
317: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_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, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
324: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
325: PetscCall(PetscSFDestroy(&sfcell));
327: /* remove points which left processor */
328: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
329: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
331: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
332: for (p = npoints_prior_migration; p < npoints2; p++) {
333: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
334: /* kill point */
335: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
336: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL)); /* you need to update npoints as the list size decreases! */
337: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
338: p--; /* check replacement point */
339: }
340: }
341: }
343: /* check for error on removed points */
344: if (error_check) {
345: PetscCall(DMSwarmGetSize(dm, &npoints2g));
346: 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);
347: }
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm, PetscBool remove_sent_points)
352: {
353: PetscFunctionBegin;
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*
358: Redundant as this assumes points can only be sent to a single rank
359: */
360: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize)
361: {
362: DM_Swarm *swarm = (DM_Swarm *)dm->data;
363: DMSwarmDataEx de;
364: PetscInt p, npoints, *rankval, n_points_recv;
365: PetscMPIInt rank, nrank, negrank;
366: void *point_buffer, *recv_points;
367: size_t sizeof_dmswarm_point;
369: PetscFunctionBegin;
370: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
371: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
372: *globalsize = npoints;
373: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
374: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
375: PetscCall(DMSwarmDataExTopologyInitialize(de));
376: for (p = 0; p < npoints; p++) {
377: negrank = rankval[p];
378: if (negrank < 0) {
379: nrank = -negrank - 1;
380: PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
381: }
382: }
383: PetscCall(DMSwarmDataExTopologyFinalize(de));
384: PetscCall(DMSwarmDataExInitializeSendCount(de));
385: for (p = 0; p < npoints; p++) {
386: negrank = rankval[p];
387: if (negrank < 0) {
388: nrank = -negrank - 1;
389: PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
390: }
391: }
392: PetscCall(DMSwarmDataExFinalizeSendCount(de));
393: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
394: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
395: for (p = 0; p < npoints; p++) {
396: negrank = rankval[p];
397: if (negrank < 0) {
398: nrank = -negrank - 1;
399: rankval[p] = nrank;
400: /* copy point into buffer */
401: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
402: /* insert point buffer into DMSwarmDataExchanger */
403: PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
404: rankval[p] = negrank;
405: }
406: }
407: PetscCall(DMSwarmDataExPackFinalize(de));
408: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
409: PetscCall(DMSwarmDataExBegin(de));
410: PetscCall(DMSwarmDataExEnd(de));
411: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
412: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
413: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
414: for (p = 0; p < n_points_recv; p++) {
415: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
417: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
418: }
419: PetscCall(DMSwarmDataExView(de));
420: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
421: PetscCall(DMSwarmDataExDestroy(de));
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: typedef struct {
426: PetscMPIInt owner_rank;
427: PetscReal min[3], max[3];
428: } CollectBBox;
430: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm, PetscInt *globalsize)
431: {
432: DM_Swarm *swarm = (DM_Swarm *)dm->data;
433: DMSwarmDataEx de;
434: PetscInt p, pk, npoints, *rankval, n_points_recv, n_bbox_recv, dim, neighbour_cells;
435: PetscMPIInt rank, nrank;
436: void *point_buffer, *recv_points;
437: size_t sizeof_dmswarm_point, sizeof_bbox_ctx;
438: PetscBool isdmda;
439: CollectBBox *bbox, *recv_bbox;
440: const PetscMPIInt *dmneighborranks;
441: DM dmcell;
443: PetscFunctionBegin;
444: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
446: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
447: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
448: isdmda = PETSC_FALSE;
449: PetscCall(PetscObjectTypeCompare((PetscObject)dmcell, DMDA, &isdmda));
450: PetscCheck(isdmda, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only DMDA support for CollectBoundingBox");
452: PetscCall(DMGetDimension(dm, &dim));
453: sizeof_bbox_ctx = sizeof(CollectBBox);
454: PetscCall(PetscMalloc1(1, &bbox));
455: bbox->owner_rank = rank;
457: /* compute the bounding box based on the overlapping / stenctil size */
458: {
459: Vec lcoor;
461: PetscCall(DMGetCoordinatesLocal(dmcell, &lcoor));
462: if (dim >= 1) {
463: PetscCall(VecStrideMin(lcoor, 0, NULL, &bbox->min[0]));
464: PetscCall(VecStrideMax(lcoor, 0, NULL, &bbox->max[0]));
465: }
466: if (dim >= 2) {
467: PetscCall(VecStrideMin(lcoor, 1, NULL, &bbox->min[1]));
468: PetscCall(VecStrideMax(lcoor, 1, NULL, &bbox->max[1]));
469: }
470: if (dim == 3) {
471: PetscCall(VecStrideMin(lcoor, 2, NULL, &bbox->min[2]));
472: PetscCall(VecStrideMax(lcoor, 2, NULL, &bbox->max[2]));
473: }
474: }
475: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
476: *globalsize = npoints;
477: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
478: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
479: /* use DMDA neighbours */
480: PetscCall(DMDAGetNeighbors(dmcell, &dmneighborranks));
481: if (dim == 1) {
482: neighbour_cells = 3;
483: } else if (dim == 2) {
484: neighbour_cells = 9;
485: } else {
486: neighbour_cells = 27;
487: }
488: PetscCall(DMSwarmDataExTopologyInitialize(de));
489: for (p = 0; p < neighbour_cells; p++) {
490: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, dmneighborranks[p]));
491: }
492: PetscCall(DMSwarmDataExTopologyFinalize(de));
493: PetscCall(DMSwarmDataExInitializeSendCount(de));
494: for (p = 0; p < neighbour_cells; p++) {
495: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExAddToSendCount(de, dmneighborranks[p], 1));
496: }
497: PetscCall(DMSwarmDataExFinalizeSendCount(de));
498: /* send bounding boxes */
499: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_bbox_ctx));
500: for (p = 0; p < neighbour_cells; p++) {
501: nrank = dmneighborranks[p];
502: if ((nrank >= 0) && (nrank != rank)) {
503: /* insert bbox buffer into DMSwarmDataExchanger */
504: PetscCall(DMSwarmDataExPackData(de, nrank, 1, bbox));
505: }
506: }
507: PetscCall(DMSwarmDataExPackFinalize(de));
508: /* recv bounding boxes */
509: PetscCall(DMSwarmDataExBegin(de));
510: PetscCall(DMSwarmDataExEnd(de));
511: PetscCall(DMSwarmDataExGetRecvData(de, &n_bbox_recv, (void **)&recv_bbox));
512: /* Wrong, should not be using PETSC_COMM_WORLD */
513: for (p = 0; p < n_bbox_recv; p++) {
514: 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],
515: (double)recv_bbox[p].max[1]));
516: }
517: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
518: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
519: PetscCall(DMSwarmDataExInitializeSendCount(de));
520: for (pk = 0; pk < n_bbox_recv; pk++) {
521: PetscReal *array_x, *array_y;
523: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
524: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
525: for (p = 0; p < npoints; p++) {
526: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
527: 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));
528: }
529: }
530: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
531: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
532: }
533: PetscCall(DMSwarmDataExFinalizeSendCount(de));
534: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
535: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
536: for (pk = 0; pk < n_bbox_recv; pk++) {
537: PetscReal *array_x, *array_y;
539: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
540: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
541: for (p = 0; p < npoints; p++) {
542: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
543: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
544: /* copy point into buffer */
545: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
546: /* insert point buffer into DMSwarmDataExchanger */
547: PetscCall(DMSwarmDataExPackData(de, recv_bbox[pk].owner_rank, 1, point_buffer));
548: }
549: }
550: }
551: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
552: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
553: }
554: PetscCall(DMSwarmDataExPackFinalize(de));
555: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
556: PetscCall(DMSwarmDataExBegin(de));
557: PetscCall(DMSwarmDataExEnd(de));
558: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
559: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
560: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
561: for (p = 0; p < n_points_recv; p++) {
562: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
564: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
565: }
566: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
567: PetscCall(PetscFree(bbox));
568: PetscCall(DMSwarmDataExView(de));
569: PetscCall(DMSwarmDataExDestroy(de));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: /* General collection when no order, or neighbour information is provided */
574: /*
575: User provides context and collect() method
576: Broadcast user context
578: for each context / rank {
579: collect(swarm,context,n,list)
580: }
581: */
582: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm, PetscErrorCode (*collect)(DM, void *, PetscInt *, PetscInt **), size_t ctx_size, void *ctx, PetscInt *globalsize)
583: {
584: DM_Swarm *swarm = (DM_Swarm *)dm->data;
585: DMSwarmDataEx de;
586: PetscInt p, r, npoints, n_points_recv;
587: PetscMPIInt size, rank;
588: void *point_buffer, *recv_points;
589: void *ctxlist;
590: PetscInt *n2collect, **collectlist;
591: size_t sizeof_dmswarm_point;
593: PetscFunctionBegin;
594: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
595: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
596: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
597: *globalsize = npoints;
598: /* Broadcast user context */
599: PetscCall(PetscMalloc(ctx_size * size, &ctxlist));
600: PetscCallMPI(MPI_Allgather(ctx, ctx_size, MPI_CHAR, ctxlist, ctx_size, MPI_CHAR, PetscObjectComm((PetscObject)dm)));
601: PetscCall(PetscMalloc1(size, &n2collect));
602: PetscCall(PetscMalloc1(size, &collectlist));
603: for (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 (r = 0; r < size; r++) {
621: if (n2collect[r] > 0) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, (PetscMPIInt)r));
622: }
623: PetscCall(DMSwarmDataExTopologyFinalize(de));
624: /* Define send counts */
625: PetscCall(DMSwarmDataExInitializeSendCount(de));
626: for (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 (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, (void **)&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 (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: }