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