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: nrank = (PetscMPIInt)rankval[p];
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: nrank = (PetscMPIInt)rankval[p];
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: nrank = (PetscMPIInt)rankval[p];
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: nrank = (PetscMPIInt)rankval[p];
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, (void **)&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, (void **)&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: PetscMPIInt size, rank;
191: PetscFunctionBegin;
192: PetscCall(DMSwarmGetCellDM(dm, &dmcell));
193: PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
195: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
196: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
198: #if 1
199: {
200: PetscInt npoints_curr, range = 0;
201: PetscSFNode *sf_cells;
203: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
204: PetscCall(PetscMalloc1(npoints_curr, &sf_cells));
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));
214: /* PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell)); */
215: PetscCall(PetscSFCreate(PETSC_COMM_SELF, &sfcell));
216: PetscCall(PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER));
217: }
218: #endif
220: PetscCall(DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos));
221: PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
222: PetscCall(DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos));
224: if (error_check) PetscCall(DMSwarmGetSize(dm, &npointsg));
225: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
226: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
227: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
228: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
230: for (p = 0; p < npoints; p++) {
231: p_cellid[p] = LA_sfcell[p].index;
232: rankval[p] = rank;
233: }
234: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
235: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
236: PetscCall(PetscSFDestroy(&sfcell));
238: if (size > 1) {
239: PetscCall(DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration));
240: } else {
241: DMSwarmDataField PField;
242: PetscInt npoints_curr;
244: /* remove points which the domain */
245: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
246: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
248: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
249: for (p = 0; p < npoints_curr; p++) {
250: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
251: /* kill point */
252: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
253: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL)); /* you need to update npoints as the list size decreases! */
254: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
255: p--; /* check replacement point */
256: }
257: }
258: PetscCall(DMSwarmGetLocalSize(dm, &npoints_prior_migration));
259: }
261: /* locate points newly received */
262: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
264: #if 0
265: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
266: PetscScalar *LA_coor;
267: PetscInt bs;
268: DMSwarmDataField PField;
270: PetscCall(DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor));
271: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos));
272: PetscCall(DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell));
274: PetscCall(VecDestroy(&pos));
275: PetscCall(DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor));
277: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
278: PetscCall(DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));
279: for (p=0; p<npoints2; p++) {
280: rankval[p] = LA_sfcell[p].index;
281: }
282: PetscCall(PetscSFDestroy(&sfcell));
283: PetscCall(DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));
285: /* remove points which left processor */
286: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField));
287: PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval));
289: PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL));
290: for (p=0; p<npoints2; p++) {
291: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
292: /* kill point */
293: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db,p));
294: PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL)); /* you need to update npoints as the list size decreases! */
295: PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval)); /* update date point increase realloc performed */
296: p--; /* check replacement point */
297: }
298: }
299: }
300: #endif
302: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
303: PetscScalar *LA_coor;
304: PetscInt npoints_from_neighbours, bs;
305: DMSwarmDataField PField;
307: npoints_from_neighbours = npoints2 - npoints_prior_migration;
309: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor));
310: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, bs, bs * npoints_from_neighbours, (const PetscScalar *)&LA_coor[bs * npoints_prior_migration], &pos));
312: PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
314: PetscCall(VecDestroy(&pos));
315: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor));
317: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
318: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
319: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
320: for (p = 0; p < npoints_from_neighbours; p++) {
321: rankval[npoints_prior_migration + p] = rank;
322: p_cellid[npoints_prior_migration + p] = LA_sfcell[p].index;
323: }
325: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
326: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
327: PetscCall(PetscSFDestroy(&sfcell));
329: /* remove points which left processor */
330: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
331: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
333: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
334: for (p = npoints_prior_migration; p < npoints2; p++) {
335: if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
336: /* kill point */
337: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
338: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL)); /* you need to update npoints as the list size decreases! */
339: PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid)); /* update date point in case realloc performed */
340: p--; /* check replacement point */
341: }
342: }
343: }
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 = (PetscMPIInt)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 = (PetscMPIInt)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 = (PetscMPIInt)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, npoints, n_points_recv;
589: PetscMPIInt size, rank, len;
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(PetscMPIIntCast(ctx_size, &len));
602: PetscCall(PetscMalloc(ctx_size * size, &ctxlist));
603: PetscCallMPI(MPI_Allgather(ctx, len, MPI_CHAR, ctxlist, len, MPI_CHAR, PetscObjectComm((PetscObject)dm)));
604: PetscCall(PetscMalloc1(size, &n2collect));
605: PetscCall(PetscMalloc1(size, &collectlist));
606: for (PetscMPIInt r = 0; r < size; r++) {
607: PetscInt _n2collect;
608: PetscInt *_collectlist;
609: void *_ctx_r;
611: _n2collect = 0;
612: _collectlist = NULL;
613: if (r != rank) { /* don't collect data from yourself */
614: _ctx_r = (void *)((char *)ctxlist + r * ctx_size);
615: PetscCall(collect(dm, _ctx_r, &_n2collect, &_collectlist));
616: }
617: n2collect[r] = _n2collect;
618: collectlist[r] = _collectlist;
619: }
620: PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
621: /* Define topology */
622: PetscCall(DMSwarmDataExTopologyInitialize(de));
623: for (PetscMPIInt r = 0; r < size; r++) {
624: if (n2collect[r] > 0) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, (PetscMPIInt)r));
625: }
626: PetscCall(DMSwarmDataExTopologyFinalize(de));
627: /* Define send counts */
628: PetscCall(DMSwarmDataExInitializeSendCount(de));
629: for (PetscMPIInt r = 0; r < size; r++) {
630: if (n2collect[r] > 0) PetscCall(DMSwarmDataExAddToSendCount(de, r, n2collect[r]));
631: }
632: PetscCall(DMSwarmDataExFinalizeSendCount(de));
633: /* Pack data */
634: PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
635: PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
636: for (PetscMPIInt r = 0; r < size; r++) {
637: for (p = 0; p < n2collect[r]; p++) {
638: PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, collectlist[r][p], point_buffer));
639: /* insert point buffer into the data exchanger */
640: PetscCall(DMSwarmDataExPackData(de, r, 1, point_buffer));
641: }
642: }
643: PetscCall(DMSwarmDataExPackFinalize(de));
644: /* Scatter */
645: PetscCall(DMSwarmDataExBegin(de));
646: PetscCall(DMSwarmDataExEnd(de));
647: /* Collect data in DMSwarm container */
648: PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
649: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
650: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
651: for (p = 0; p < n_points_recv; p++) {
652: void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);
654: PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
655: }
656: /* Release memory */
657: for (PetscMPIInt r = 0; r < size; r++) {
658: if (collectlist[r]) PetscCall(PetscFree(collectlist[r]));
659: }
660: PetscCall(PetscFree(collectlist));
661: PetscCall(PetscFree(n2collect));
662: PetscCall(PetscFree(ctxlist));
663: PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
664: PetscCall(DMSwarmDataExView(de));
665: PetscCall(DMSwarmDataExDestroy(de));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: DMSwarmGetMigrateType - Get the style of point migration
672: Logically Collective
674: Input Parameter:
675: . dm - the `DMSWARM`
677: Output Parameter:
678: . mtype - The migration type, see `DMSwarmMigrateType`
680: Level: intermediate
682: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmMigrate()`
683: @*/
684: PetscErrorCode DMSwarmGetMigrateType(DM dm, DMSwarmMigrateType *mtype)
685: {
686: PetscFunctionBegin;
688: PetscAssertPointer(mtype, 2);
689: *mtype = ((DM_Swarm *)dm->data)->migrate_type;
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: /*@
694: DMSwarmSetMigrateType - Set the style of point migration
696: Logically Collective
698: Input Parameters:
699: + dm - the `DMSWARM`
700: - mtype - The migration type, see `DMSwarmMigrateType`
702: Level: intermediate
704: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmGetMigrateType()`, `DMSwarmMigrate()`
705: @*/
706: PetscErrorCode DMSwarmSetMigrateType(DM dm, DMSwarmMigrateType mtype)
707: {
708: PetscFunctionBegin;
711: ((DM_Swarm *)dm->data)->migrate_type = mtype;
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }