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