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: DMSwarmDataEx de;
101: PetscInt r, p, npoints, *p_cellid, n_points_recv;
102: PetscMPIInt rank, _rank;
103: const PetscMPIInt *neighbourranks;
104: void *point_buffer, *recv_points;
105: size_t sizeof_dmswarm_point;
106: PetscInt nneighbors;
107: PetscMPIInt mynneigh, *myneigh;
109: PetscFunctionBegin;
110: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
111: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
112: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
113: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
114: PetscCall(DMGetNeighbors(dmcell, &nneighbors, &neighbourranks));
115: PetscCall(DMSwarmDataExTopologyInitialize(de));
116: for (r = 0; r < nneighbors; r++) {
117: _rank = neighbourranks[r];
118: if ((_rank != rank) && (_rank > 0)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, _rank));
119: }
120: PetscCall(DMSwarmDataExTopologyFinalize(de));
121: PetscCall(DMSwarmDataExTopologyGetNeighbours(de, &mynneigh, &myneigh));
122: PetscCall(DMSwarmDataExInitializeSendCount(de));
123: for (p = 0; p < npoints; p++) {
124: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
125: for (r = 0; r < mynneigh; r++) {
126: _rank = myneigh[r];
127: PetscCall(DMSwarmDataExAddToSendCount(de, _rank, 1));
128: }
129: }
130: }
131: PetscCall(DMSwarmDataExFinalizeSendCount(de));
132: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
133: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
134: for (p = 0; p < npoints; p++) {
135: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
136: for (r = 0; r < mynneigh; r++) {
137: _rank = myneigh[r];
138: /* copy point into buffer */
139: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
140: /* insert point buffer into DMSwarmDataExchanger */
141: PetscCall(DMSwarmDataExPackData(de, _rank, 1, point_buffer));
142: }
143: }
144: }
145: PetscCall(DMSwarmDataExPackFinalize(de));
146: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
147: if (remove_sent_points) {
148: DMSwarmDataField PField;
150: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
151: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
152: /* remove points which left processor */
153: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
154: for (p = 0; p < npoints; p++) {
155: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
156: /* kill point */
157: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
158: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL)); /* you need to update npoints as the list size decreases! */
159: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point increase realloc performed */
160: p--; /* check replacement point */
161: }
162: }
163: }
164: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, npoints_prior_migration, NULL, NULL));
165: PetscCall(DMSwarmDataExBegin(de));
166: PetscCall(DMSwarmDataExEnd(de));
167: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
168: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
169: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
170: for (p = 0; p < n_points_recv; p++) {
171: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
173: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
174: }
175: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
176: PetscCall(DMSwarmDataExDestroy(de));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm, PetscBool remove_sent_points)
181: {
182: DM_Swarm *swarm = (DM_Swarm *)dm->data;
183: PetscInt p, npoints, npointsg = 0, npoints2, npoints2g, *rankval, *p_cellid, npoints_prior_migration;
184: PetscSF sfcell = NULL;
185: const PetscSFNode *LA_sfcell;
186: DM dmcell;
187: Vec pos;
188: PetscBool error_check = swarm->migrate_error_on_missing_point;
189: const char *coordname;
190: PetscMPIInt size, rank;
192: PetscFunctionBegin;
193: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
194: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
196: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
197: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
199: #if 1
200: {
201: PetscInt npoints_curr, range = 0;
202: PetscSFNode *sf_cells;
204: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
205: PetscCall(PetscMalloc1(npoints_curr, &sf_cells));
207: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
208: for (p = 0; p < npoints_curr; p++) {
209: sf_cells[p].rank = 0;
210: sf_cells[p].index = p_cellid[p];
211: if (p_cellid[p] > range) range = p_cellid[p];
212: }
213: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
215: /* PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell)); */
216: PetscCall(PetscSFCreate(PETSC_COMM_SELF, &sfcell));
217: PetscCall(PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER));
218: }
219: #endif
221: PetscCall(DMSwarmGetCoordinateField(dm, &coordname));
222: PetscCall(DMSwarmCreateLocalVectorFromField(dm, coordname, &pos));
223: PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
224: PetscCall(DMSwarmDestroyLocalVectorFromField(dm, coordname, &pos));
226: if (error_check) PetscCall(DMSwarmGetSize(dm, &npointsg));
227: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
228: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
229: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
230: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
232: for (p = 0; p < npoints; p++) {
233: p_cellid[p] = LA_sfcell[p].index;
234: rankval[p] = rank;
235: }
236: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
237: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
238: PetscCall(PetscSFDestroy(&sfcell));
240: if (size > 1) {
241: PetscCall(DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration));
242: } else {
243: DMSwarmDataField PField;
244: PetscInt npoints_curr;
246: /* remove points which the domain */
247: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
248: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
250: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
251: for (p = 0; p < npoints_curr; p++) {
252: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
253: /* kill point */
254: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
255: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL)); /* you need to update npoints as the list size decreases! */
256: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
257: p--; /* check replacement point */
258: }
259: }
260: PetscCall(DMSwarmGetLocalSize(dm, &npoints_prior_migration));
261: }
263: /* locate points newly received */
264: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
266: #if 0
267: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
268: PetscScalar *LA_coor;
269: PetscInt bs;
270: DMSwarmDataField PField;
272: PetscCall(DMSwarmGetField(dm,coordname,&bs,NULL,(void**)&LA_coor));
273: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos));
274: PetscCall(DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell));
276: PetscCall(VecDestroy(&pos));
277: PetscCall(DMSwarmRestoreField(dm,coordname,&bs,NULL,(void**)&LA_coor));
279: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
280: PetscCall(DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));
281: for (p=0; p<npoints2; p++) {
282: rankval[p] = LA_sfcell[p].index;
283: }
284: PetscCall(PetscSFDestroy(&sfcell));
285: PetscCall(DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));
287: /* remove points which left processor */
288: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField));
289: PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval));
291: PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL));
292: for (p=0; p<npoints2; p++) {
293: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
294: /* kill point */
295: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db,p));
296: PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL)); /* you need to update npoints as the list size decreases! */
297: PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval)); /* update date point increase realloc performed */
298: p--; /* check replacement point */
299: }
300: }
301: }
302: #endif
304: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
305: PetscScalar *LA_coor;
306: PetscInt npoints_from_neighbours, bs;
307: DMSwarmDataField PField;
309: npoints_from_neighbours = npoints2 - npoints_prior_migration;
311: PetscCall(DMSwarmGetField(dm, coordname, &bs, NULL, (void **)&LA_coor));
312: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, bs, bs * npoints_from_neighbours, (const PetscScalar *)&LA_coor[bs * npoints_prior_migration], &pos));
314: PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
316: PetscCall(VecDestroy(&pos));
317: PetscCall(DMSwarmRestoreField(dm, coordname, &bs, NULL, (void **)&LA_coor));
319: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
320: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
321: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
322: for (p = 0; p < npoints_from_neighbours; p++) {
323: rankval[npoints_prior_migration + p] = rank;
324: p_cellid[npoints_prior_migration + p] = LA_sfcell[p].index;
325: }
327: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
328: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
329: PetscCall(PetscSFDestroy(&sfcell));
331: /* remove points which left processor */
332: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
333: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
335: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
336: for (p = npoints_prior_migration; p < npoints2; p++) {
337: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
338: /* kill point */
339: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
340: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL)); /* you need to update npoints as the list size decreases! */
341: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
342: p--; /* check replacement point */
343: }
344: }
345: }
347: /* check for error on removed points */
348: if (error_check) {
349: PetscCall(DMSwarmGetSize(dm, &npoints2g));
350: 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);
351: }
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm, PetscBool remove_sent_points)
356: {
357: PetscFunctionBegin;
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*
362: Redundant as this assumes points can only be sent to a single rank
363: */
364: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize)
365: {
366: DM_Swarm *swarm = (DM_Swarm *)dm->data;
367: DMSwarmDataEx de;
368: PetscInt p, npoints, *rankval, n_points_recv;
369: PetscMPIInt rank, nrank, negrank;
370: void *point_buffer, *recv_points;
371: size_t sizeof_dmswarm_point;
373: PetscFunctionBegin;
374: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
375: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
376: *globalsize = npoints;
377: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
378: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
379: PetscCall(DMSwarmDataExTopologyInitialize(de));
380: for (p = 0; p < npoints; p++) {
381: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
382: if (negrank < 0) {
383: nrank = -negrank - 1;
384: PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
385: }
386: }
387: PetscCall(DMSwarmDataExTopologyFinalize(de));
388: PetscCall(DMSwarmDataExInitializeSendCount(de));
389: for (p = 0; p < npoints; p++) {
390: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
391: if (negrank < 0) {
392: nrank = -negrank - 1;
393: PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
394: }
395: }
396: PetscCall(DMSwarmDataExFinalizeSendCount(de));
397: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
398: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
399: for (p = 0; p < npoints; p++) {
400: PetscCall(PetscMPIIntCast(rankval[p], &negrank));
401: if (negrank < 0) {
402: nrank = -negrank - 1;
403: rankval[p] = nrank;
404: /* copy point into buffer */
405: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
406: /* insert point buffer into DMSwarmDataExchanger */
407: PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
408: rankval[p] = negrank;
409: }
410: }
411: PetscCall(DMSwarmDataExPackFinalize(de));
412: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
413: PetscCall(DMSwarmDataExBegin(de));
414: PetscCall(DMSwarmDataExEnd(de));
415: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
416: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
417: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
418: for (p = 0; p < n_points_recv; p++) {
419: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
421: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
422: }
423: PetscCall(DMSwarmDataExView(de));
424: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
425: PetscCall(DMSwarmDataExDestroy(de));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: typedef struct {
430: PetscMPIInt owner_rank;
431: PetscReal min[3], max[3];
432: } CollectBBox;
434: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm, PetscInt *globalsize)
435: {
436: DM_Swarm *swarm = (DM_Swarm *)dm->data;
437: DMSwarmDataEx de;
438: PetscInt p, pk, npoints, *rankval, n_points_recv, n_bbox_recv, dim, neighbour_cells;
439: PetscMPIInt rank, nrank;
440: void *point_buffer, *recv_points;
441: size_t sizeof_dmswarm_point, sizeof_bbox_ctx;
442: PetscBool isdmda;
443: CollectBBox *bbox, *recv_bbox;
444: const PetscMPIInt *dmneighborranks;
445: DM dmcell;
447: PetscFunctionBegin;
448: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
450: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
451: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
452: isdmda = PETSC_FALSE;
453: PetscCall(PetscObjectTypeCompare((PetscObject)dmcell, DMDA, &isdmda));
454: PetscCheck(isdmda, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only DMDA support for CollectBoundingBox");
456: PetscCall(DMGetDimension(dm, &dim));
457: sizeof_bbox_ctx = sizeof(CollectBBox);
458: PetscCall(PetscMalloc1(1, &bbox));
459: bbox->owner_rank = rank;
461: /* compute the bounding box based on the overlapping / stenctil size */
462: {
463: Vec lcoor;
465: PetscCall(DMGetCoordinatesLocal(dmcell, &lcoor));
466: if (dim >= 1) {
467: PetscCall(VecStrideMin(lcoor, 0, NULL, &bbox->min[0]));
468: PetscCall(VecStrideMax(lcoor, 0, NULL, &bbox->max[0]));
469: }
470: if (dim >= 2) {
471: PetscCall(VecStrideMin(lcoor, 1, NULL, &bbox->min[1]));
472: PetscCall(VecStrideMax(lcoor, 1, NULL, &bbox->max[1]));
473: }
474: if (dim == 3) {
475: PetscCall(VecStrideMin(lcoor, 2, NULL, &bbox->min[2]));
476: PetscCall(VecStrideMax(lcoor, 2, NULL, &bbox->max[2]));
477: }
478: }
479: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
480: *globalsize = npoints;
481: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
482: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
483: /* use DMDA neighbours */
484: PetscCall(DMDAGetNeighbors(dmcell, &dmneighborranks));
485: if (dim == 1) {
486: neighbour_cells = 3;
487: } else if (dim == 2) {
488: neighbour_cells = 9;
489: } else {
490: neighbour_cells = 27;
491: }
492: PetscCall(DMSwarmDataExTopologyInitialize(de));
493: for (p = 0; p < neighbour_cells; p++) {
494: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, dmneighborranks[p]));
495: }
496: PetscCall(DMSwarmDataExTopologyFinalize(de));
497: PetscCall(DMSwarmDataExInitializeSendCount(de));
498: for (p = 0; p < neighbour_cells; p++) {
499: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExAddToSendCount(de, dmneighborranks[p], 1));
500: }
501: PetscCall(DMSwarmDataExFinalizeSendCount(de));
502: /* send bounding boxes */
503: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_bbox_ctx));
504: for (p = 0; p < neighbour_cells; p++) {
505: nrank = dmneighborranks[p];
506: if ((nrank >= 0) && (nrank != rank)) {
507: /* insert bbox buffer into DMSwarmDataExchanger */
508: PetscCall(DMSwarmDataExPackData(de, nrank, 1, bbox));
509: }
510: }
511: PetscCall(DMSwarmDataExPackFinalize(de));
512: /* recv bounding boxes */
513: PetscCall(DMSwarmDataExBegin(de));
514: PetscCall(DMSwarmDataExEnd(de));
515: PetscCall(DMSwarmDataExGetRecvData(de, &n_bbox_recv, (void **)&recv_bbox));
516: /* Wrong, should not be using PETSC_COMM_WORLD */
517: for (p = 0; p < n_bbox_recv; p++) {
518: 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],
519: (double)recv_bbox[p].max[1]));
520: }
521: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
522: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
523: PetscCall(DMSwarmDataExInitializeSendCount(de));
524: for (pk = 0; pk < n_bbox_recv; pk++) {
525: PetscReal *array_x, *array_y;
527: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
528: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
529: for (p = 0; p < npoints; p++) {
530: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
531: 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));
532: }
533: }
534: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
535: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
536: }
537: PetscCall(DMSwarmDataExFinalizeSendCount(de));
538: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
539: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
540: for (pk = 0; pk < n_bbox_recv; pk++) {
541: PetscReal *array_x, *array_y;
543: PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
544: PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
545: for (p = 0; p < npoints; p++) {
546: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
547: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
548: /* copy point into buffer */
549: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
550: /* insert point buffer into DMSwarmDataExchanger */
551: PetscCall(DMSwarmDataExPackData(de, recv_bbox[pk].owner_rank, 1, point_buffer));
552: }
553: }
554: }
555: PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
556: PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
557: }
558: PetscCall(DMSwarmDataExPackFinalize(de));
559: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
560: PetscCall(DMSwarmDataExBegin(de));
561: PetscCall(DMSwarmDataExEnd(de));
562: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
563: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
564: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
565: for (p = 0; p < n_points_recv; p++) {
566: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
568: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
569: }
570: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
571: PetscCall(PetscFree(bbox));
572: PetscCall(DMSwarmDataExView(de));
573: PetscCall(DMSwarmDataExDestroy(de));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /* General collection when no order, or neighbour information is provided */
578: /*
579: User provides context and collect() method
580: Broadcast user context
582: for each context / rank {
583: collect(swarm,context,n,list)
584: }
585: */
586: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm, PetscErrorCode (*collect)(DM, void *, PetscInt *, PetscInt **), size_t ctx_size, void *ctx, PetscInt *globalsize)
587: {
588: DM_Swarm *swarm = (DM_Swarm *)dm->data;
589: DMSwarmDataEx de;
590: PetscInt p, npoints, n_points_recv;
591: PetscMPIInt size, rank, len;
592: void *point_buffer, *recv_points;
593: void *ctxlist;
594: PetscInt *n2collect, **collectlist;
595: size_t sizeof_dmswarm_point;
597: PetscFunctionBegin;
598: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
599: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
600: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
601: *globalsize = npoints;
602: /* Broadcast user context */
603: PetscCall(PetscMPIIntCast(ctx_size, &len));
604: PetscCall(PetscMalloc(ctx_size * size, &ctxlist));
605: PetscCallMPI(MPI_Allgather(ctx, len, MPI_CHAR, ctxlist, len, MPI_CHAR, PetscObjectComm((PetscObject)dm)));
606: PetscCall(PetscMalloc1(size, &n2collect));
607: PetscCall(PetscMalloc1(size, &collectlist));
608: for (PetscMPIInt r = 0; r < size; r++) {
609: PetscInt _n2collect;
610: PetscInt *_collectlist;
611: void *_ctx_r;
613: _n2collect = 0;
614: _collectlist = NULL;
615: if (r != rank) { /* don't collect data from yourself */
616: _ctx_r = (void *)((char *)ctxlist + r * ctx_size);
617: PetscCall(collect(dm, _ctx_r, &_n2collect, &_collectlist));
618: }
619: n2collect[r] = _n2collect;
620: collectlist[r] = _collectlist;
621: }
622: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
623: /* Define topology */
624: PetscCall(DMSwarmDataExTopologyInitialize(de));
625: for (PetscMPIInt r = 0; r < size; r++) {
626: if (n2collect[r] > 0) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, r));
627: }
628: PetscCall(DMSwarmDataExTopologyFinalize(de));
629: /* Define send counts */
630: PetscCall(DMSwarmDataExInitializeSendCount(de));
631: for (PetscMPIInt r = 0; r < size; r++) {
632: if (n2collect[r] > 0) PetscCall(DMSwarmDataExAddToSendCount(de, r, n2collect[r]));
633: }
634: PetscCall(DMSwarmDataExFinalizeSendCount(de));
635: /* Pack data */
636: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
637: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
638: for (PetscMPIInt r = 0; r < size; r++) {
639: for (p = 0; p < n2collect[r]; p++) {
640: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, collectlist[r][p], point_buffer));
641: /* insert point buffer into the data exchanger */
642: PetscCall(DMSwarmDataExPackData(de, r, 1, point_buffer));
643: }
644: }
645: PetscCall(DMSwarmDataExPackFinalize(de));
646: /* Scatter */
647: PetscCall(DMSwarmDataExBegin(de));
648: PetscCall(DMSwarmDataExEnd(de));
649: /* Collect data in DMSwarm container */
650: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, &recv_points));
651: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
652: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
653: for (p = 0; p < n_points_recv; p++) {
654: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
656: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
657: }
658: /* Release memory */
659: for (PetscMPIInt r = 0; r < size; r++) {
660: if (collectlist[r]) PetscCall(PetscFree(collectlist[r]));
661: }
662: PetscCall(PetscFree(collectlist));
663: PetscCall(PetscFree(n2collect));
664: PetscCall(PetscFree(ctxlist));
665: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
666: PetscCall(DMSwarmDataExView(de));
667: PetscCall(DMSwarmDataExDestroy(de));
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: /*@
672: DMSwarmGetMigrateType - Get the style of point migration
674: Logically Collective
676: Input Parameter:
677: . dm - the `DMSWARM`
679: Output Parameter:
680: . mtype - The migration type, see `DMSwarmMigrateType`
682: Level: intermediate
684: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmMigrate()`
685: @*/
686: PetscErrorCode DMSwarmGetMigrateType(DM dm, DMSwarmMigrateType *mtype)
687: {
688: PetscFunctionBegin;
690: PetscAssertPointer(mtype, 2);
691: *mtype = ((DM_Swarm *)dm->data)->migrate_type;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@
696: DMSwarmSetMigrateType - Set the style of point migration
698: Logically Collective
700: Input Parameters:
701: + dm - the `DMSWARM`
702: - mtype - The migration type, see `DMSwarmMigrateType`
704: Level: intermediate
706: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmGetMigrateType()`, `DMSwarmMigrate()`
707: @*/
708: PetscErrorCode DMSwarmSetMigrateType(DM dm, DMSwarmMigrateType mtype)
709: {
710: PetscFunctionBegin;
713: ((DM_Swarm *)dm->data)->migrate_type = mtype;
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }