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: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
23: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
24: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
25: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0, &de);
26: DMSwarmDataExTopologyInitialize(de);
27: for (p = 0; p < npoints; ++p) {
28: nrank = rankval[p];
29: if (nrank != rank) {
30: DMSwarmDataExTopologyAddNeighbour(de,nrank);
31: }
32: }
33: DMSwarmDataExTopologyFinalize(de);
34: DMSwarmDataExInitializeSendCount(de);
35: for (p=0; p<npoints; p++) {
36: nrank = rankval[p];
37: if (nrank != rank) {
38: DMSwarmDataExAddToSendCount(de,nrank,1);
39: }
40: }
41: DMSwarmDataExFinalizeSendCount(de);
42: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
43: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
44: for (p=0; p<npoints; p++) {
45: nrank = rankval[p];
46: if (nrank != rank) {
47: /* copy point into buffer */
48: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
49: /* insert point buffer into DMSwarmDataExchanger */
50: DMSwarmDataExPackData(de,nrank,1,point_buffer);
51: }
52: }
53: DMSwarmDataExPackFinalize(de);
54: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
56: if (remove_sent_points) {
57: DMSwarmDataField gfield;
59: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);
60: DMSwarmDataFieldGetAccess(gfield);
61: DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);
63: /* remove points which left processor */
64: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
65: for (p=0; p<npoints; p++) {
66: nrank = rankval[p];
67: if (nrank != rank) {
68: /* kill point */
69: DMSwarmDataFieldRestoreAccess(gfield);
71: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
73: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
74: DMSwarmDataFieldGetAccess(gfield);
75: DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);
76: p--; /* check replacement point */
77: }
78: }
79: DMSwarmDataFieldRestoreEntries(gfield,(void**)&rankval);
80: DMSwarmDataFieldRestoreAccess(gfield);
81: }
82: DMSwarmDataExBegin(de);
83: DMSwarmDataExEnd(de);
84: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
85: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
86: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
87: for (p=0; p<n_points_recv; p++) {
88: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
90: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
91: }
92: if (debug) DMSwarmDataExView(de);
93: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
94: DMSwarmDataExDestroy(de);
95: return 0;
96: }
98: PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm,DM dmcell,PetscBool remove_sent_points,PetscInt *npoints_prior_migration)
99: {
100: DM_Swarm *swarm = (DM_Swarm*)dm->data;
101: DMSwarmDataEx de;
102: PetscInt r,p,npoints,*rankval,n_points_recv;
103: PetscMPIInt rank,_rank;
104: const PetscMPIInt *neighbourranks;
105: void *point_buffer,*recv_points;
106: size_t sizeof_dmswarm_point;
107: PetscInt nneighbors;
108: PetscMPIInt mynneigh,*myneigh;
110: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
111: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
112: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
113: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
114: DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);
115: DMSwarmDataExTopologyInitialize(de);
116: for (r=0; r<nneighbors; r++) {
117: _rank = neighbourranks[r];
118: if ((_rank != rank) && (_rank > 0)) {
119: DMSwarmDataExTopologyAddNeighbour(de,_rank);
120: }
121: }
122: DMSwarmDataExTopologyFinalize(de);
123: DMSwarmDataExTopologyGetNeighbours(de,&mynneigh,&myneigh);
124: DMSwarmDataExInitializeSendCount(de);
125: for (p=0; p<npoints; p++) {
126: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
127: for (r=0; r<mynneigh; r++) {
128: _rank = myneigh[r];
129: DMSwarmDataExAddToSendCount(de,_rank,1);
130: }
131: }
132: }
133: DMSwarmDataExFinalizeSendCount(de);
134: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
135: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
136: for (p=0; p<npoints; p++) {
137: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
138: for (r=0; r<mynneigh; r++) {
139: _rank = myneigh[r];
140: /* copy point into buffer */
141: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
142: /* insert point buffer into DMSwarmDataExchanger */
143: DMSwarmDataExPackData(de,_rank,1,point_buffer);
144: }
145: }
146: }
147: DMSwarmDataExPackFinalize(de);
148: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
149: if (remove_sent_points) {
150: DMSwarmDataField PField;
152: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
153: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
154: /* remove points which left processor */
155: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
156: for (p=0; p<npoints; p++) {
157: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
158: /* kill point */
159: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
160: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
161: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
162: p--; /* check replacement point */
163: }
164: }
165: }
166: DMSwarmDataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);
167: DMSwarmDataExBegin(de);
168: DMSwarmDataExEnd(de);
169: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
170: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
171: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
172: for (p=0; p<n_points_recv; p++) {
173: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
175: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
176: }
177: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
178: DMSwarmDataExDestroy(de);
179: return 0;
180: }
182: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points)
183: {
184: DM_Swarm *swarm = (DM_Swarm*)dm->data;
185: PetscInt p,npoints,npointsg=0,npoints2,npoints2g,*rankval,npoints_prior_migration;
186: PetscSF sfcell = NULL;
187: const PetscSFNode *LA_sfcell;
188: DM dmcell;
189: Vec pos;
190: PetscBool error_check = swarm->migrate_error_on_missing_point;
191: PetscMPIInt size,rank;
193: DMSwarmGetCellDM(dm,&dmcell);
196: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
197: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
199: #if 1
200: {
201: PetscInt *p_cellid;
202: PetscInt npoints_curr,range = 0;
203: PetscSFNode *sf_cells;
205: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
206: PetscMalloc1(npoints_curr, &sf_cells);
208: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
209: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
210: for (p=0; p<npoints_curr; p++) {
212: sf_cells[p].rank = 0;
213: sf_cells[p].index = p_cellid[p];
214: if (p_cellid[p] > range) {
215: range = p_cellid[p];
216: }
217: }
218: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
219: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
221: /* PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell); */
222: PetscSFCreate(PETSC_COMM_SELF,&sfcell);
223: PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
224: }
225: #endif
227: DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
228: DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
229: DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
231: if (error_check) {
232: DMSwarmGetSize(dm,&npointsg);
233: }
234: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
235: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
236: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
237: for (p=0; p<npoints; p++) {
238: rankval[p] = LA_sfcell[p].index;
239: }
240: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
241: PetscSFDestroy(&sfcell);
243: if (size > 1) {
244: DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);
245: } else {
246: DMSwarmDataField PField;
247: PetscInt npoints_curr;
249: /* remove points which the domain */
250: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
251: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
253: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
254: for (p=0; p<npoints_curr; p++) {
255: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
256: /* kill point */
257: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
258: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL); /* you need to update npoints as the list size decreases! */
259: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
260: p--; /* check replacement point */
261: }
262: }
263: DMSwarmGetSize(dm,&npoints_prior_migration);
265: }
267: /* locate points newly received */
268: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
270: #if 0
271: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
272: PetscScalar *LA_coor;
273: PetscInt bs;
274: DMSwarmDataField PField;
276: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
277: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);
278: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
280: VecDestroy(&pos);
281: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
283: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
284: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
285: for (p=0; p<npoints2; p++) {
286: rankval[p] = LA_sfcell[p].index;
287: }
288: PetscSFDestroy(&sfcell);
289: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
291: /* remove points which left processor */
292: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
293: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
295: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
296: for (p=0; p<npoints2; p++) {
297: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
298: /* kill point */
299: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
300: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
301: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
302: p--; /* check replacement point */
303: }
304: }
305: }
306: #endif
308: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
309: PetscScalar *LA_coor;
310: PetscInt npoints_from_neighbours,bs;
311: DMSwarmDataField PField;
313: npoints_from_neighbours = npoints2 - npoints_prior_migration;
315: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
316: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);
318: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
320: VecDestroy(&pos);
321: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
323: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
324: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
325: for (p=0; p<npoints_from_neighbours; p++) {
326: rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
327: }
328: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
329: PetscSFDestroy(&sfcell);
331: /* remove points which left processor */
332: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
333: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
335: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
336: for (p=npoints_prior_migration; p<npoints2; p++) {
337: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
338: /* kill point */
339: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
340: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
341: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
342: p--; /* check replacement point */
343: }
344: }
345: }
347: {
348: PetscInt *p_cellid;
350: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
351: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
352: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
353: for (p=0; p<npoints2; p++) {
354: p_cellid[p] = rankval[p];
355: }
356: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
357: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
358: }
360: /* check for error on removed points */
361: if (error_check) {
362: DMSwarmGetSize(dm,&npoints2g);
364: }
365: return 0;
366: }
368: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
369: {
370: return 0;
371: }
373: /*
374: Redundant as this assumes points can only be sent to a single rank
375: */
376: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
377: {
378: DM_Swarm *swarm = (DM_Swarm*)dm->data;
379: DMSwarmDataEx de;
380: PetscInt p,npoints,*rankval,n_points_recv;
381: PetscMPIInt rank,nrank,negrank;
382: void *point_buffer,*recv_points;
383: size_t sizeof_dmswarm_point;
385: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
386: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
387: *globalsize = npoints;
388: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
389: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
390: DMSwarmDataExTopologyInitialize(de);
391: for (p=0; p<npoints; p++) {
392: negrank = rankval[p];
393: if (negrank < 0) {
394: nrank = -negrank - 1;
395: DMSwarmDataExTopologyAddNeighbour(de,nrank);
396: }
397: }
398: DMSwarmDataExTopologyFinalize(de);
399: DMSwarmDataExInitializeSendCount(de);
400: for (p=0; p<npoints; p++) {
401: negrank = rankval[p];
402: if (negrank < 0) {
403: nrank = -negrank - 1;
404: DMSwarmDataExAddToSendCount(de,nrank,1);
405: }
406: }
407: DMSwarmDataExFinalizeSendCount(de);
408: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
409: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
410: for (p=0; p<npoints; p++) {
411: negrank = rankval[p];
412: if (negrank < 0) {
413: nrank = -negrank - 1;
414: rankval[p] = nrank;
415: /* copy point into buffer */
416: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
417: /* insert point buffer into DMSwarmDataExchanger */
418: DMSwarmDataExPackData(de,nrank,1,point_buffer);
419: rankval[p] = negrank;
420: }
421: }
422: DMSwarmDataExPackFinalize(de);
423: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
424: DMSwarmDataExBegin(de);
425: DMSwarmDataExEnd(de);
426: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
427: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
428: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
429: for (p=0; p<n_points_recv; p++) {
430: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
432: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
433: }
434: DMSwarmDataExView(de);
435: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
436: DMSwarmDataExDestroy(de);
437: return 0;
438: }
440: typedef struct {
441: PetscMPIInt owner_rank;
442: PetscReal min[3],max[3];
443: } CollectBBox;
445: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
446: {
447: DM_Swarm * swarm = (DM_Swarm*)dm->data;
448: DMSwarmDataEx de;
449: PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
450: PetscMPIInt rank,nrank;
451: void *point_buffer,*recv_points;
452: size_t sizeof_dmswarm_point,sizeof_bbox_ctx;
453: PetscBool isdmda;
454: CollectBBox *bbox,*recv_bbox;
455: const PetscMPIInt *dmneighborranks;
456: DM dmcell;
458: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
460: DMSwarmGetCellDM(dm,&dmcell);
462: isdmda = PETSC_FALSE;
463: PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
466: DMGetDimension(dm,&dim);
467: sizeof_bbox_ctx = sizeof(CollectBBox);
468: PetscMalloc1(1,&bbox);
469: bbox->owner_rank = rank;
471: /* compute the bounding box based on the overlapping / stenctil size */
472: {
473: Vec lcoor;
475: DMGetCoordinatesLocal(dmcell,&lcoor);
476: if (dim >= 1) {
477: VecStrideMin(lcoor,0,NULL,&bbox->min[0]);
478: VecStrideMax(lcoor,0,NULL,&bbox->max[0]);
479: }
480: if (dim >= 2) {
481: VecStrideMin(lcoor,1,NULL,&bbox->min[1]);
482: VecStrideMax(lcoor,1,NULL,&bbox->max[1]);
483: }
484: if (dim == 3) {
485: VecStrideMin(lcoor,2,NULL,&bbox->min[2]);
486: VecStrideMax(lcoor,2,NULL,&bbox->max[2]);
487: }
488: }
489: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
490: *globalsize = npoints;
491: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
492: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
493: /* use DMDA neighbours */
494: DMDAGetNeighbors(dmcell,&dmneighborranks);
495: if (dim == 1) {
496: neighbour_cells = 3;
497: } else if (dim == 2) {
498: neighbour_cells = 9;
499: } else {
500: neighbour_cells = 27;
501: }
502: DMSwarmDataExTopologyInitialize(de);
503: for (p=0; p<neighbour_cells; p++) {
504: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
505: DMSwarmDataExTopologyAddNeighbour(de,dmneighborranks[p]);
506: }
507: }
508: DMSwarmDataExTopologyFinalize(de);
509: DMSwarmDataExInitializeSendCount(de);
510: for (p=0; p<neighbour_cells; p++) {
511: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
512: DMSwarmDataExAddToSendCount(de,dmneighborranks[p],1);
513: }
514: }
515: DMSwarmDataExFinalizeSendCount(de);
516: /* send bounding boxes */
517: DMSwarmDataExPackInitialize(de,sizeof_bbox_ctx);
518: for (p=0; p<neighbour_cells; p++) {
519: nrank = dmneighborranks[p];
520: if ((nrank >= 0) && (nrank != rank)) {
521: /* insert bbox buffer into DMSwarmDataExchanger */
522: DMSwarmDataExPackData(de,nrank,1,bbox);
523: }
524: }
525: DMSwarmDataExPackFinalize(de);
526: /* recv bounding boxes */
527: DMSwarmDataExBegin(de);
528: DMSwarmDataExEnd(de);
529: DMSwarmDataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);
530: /* Wrong, should not be using PETSC_COMM_WORLD */
531: for (p=0; p<n_bbox_recv; p++) {
532: 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,
533: (double)recv_bbox[p].min[0],(double)recv_bbox[p].max[0],(double)recv_bbox[p].min[1],(double)recv_bbox[p].max[1]));
534: }
535: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
536: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
537: DMSwarmDataExInitializeSendCount(de);
538: for (pk=0; pk<n_bbox_recv; pk++) {
539: PetscReal *array_x,*array_y;
541: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
542: 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: DMSwarmDataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);
547: }
548: }
549: }
550: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
551: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
552: }
553: DMSwarmDataExFinalizeSendCount(de);
554: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
555: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
556: for (pk=0; pk<n_bbox_recv; pk++) {
557: PetscReal *array_x,*array_y;
559: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
560: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
561: for (p=0; p<npoints; p++) {
562: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
563: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
564: /* copy point into buffer */
565: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
566: /* insert point buffer into DMSwarmDataExchanger */
567: DMSwarmDataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);
568: }
569: }
570: }
571: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
572: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
573: }
574: DMSwarmDataExPackFinalize(de);
575: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
576: DMSwarmDataExBegin(de);
577: DMSwarmDataExEnd(de);
578: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
579: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
580: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
581: for (p=0; p<n_points_recv; p++) {
582: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
584: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
585: }
586: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
587: PetscFree(bbox);
588: DMSwarmDataExView(de);
589: DMSwarmDataExDestroy(de);
590: return 0;
591: }
593: /* General collection when no order, or neighbour information is provided */
594: /*
595: User provides context and collect() method
596: Broadcast user context
598: for each context / rank {
599: collect(swarm,context,n,list)
600: }
601: */
602: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
603: {
604: DM_Swarm *swarm = (DM_Swarm*)dm->data;
605: DMSwarmDataEx de;
606: PetscInt p,r,npoints,n_points_recv;
607: PetscMPIInt size,rank;
608: void *point_buffer,*recv_points;
609: void *ctxlist;
610: PetscInt *n2collect,**collectlist;
611: size_t sizeof_dmswarm_point;
613: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
614: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
615: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
616: *globalsize = npoints;
617: /* Broadcast user context */
618: PetscMalloc(ctx_size*size,&ctxlist);
619: MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));
620: PetscMalloc1(size,&n2collect);
621: PetscMalloc1(size,&collectlist);
622: for (r=0; r<size; r++) {
623: PetscInt _n2collect;
624: PetscInt *_collectlist;
625: void *_ctx_r;
627: _n2collect = 0;
628: _collectlist = NULL;
629: if (r != rank) { /* don't collect data from yourself */
630: _ctx_r = (void*)( (char*)ctxlist + r * ctx_size);
631: collect(dm,_ctx_r,&_n2collect,&_collectlist);
632: }
633: n2collect[r] = _n2collect;
634: collectlist[r] = _collectlist;
635: }
636: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
637: /* Define topology */
638: DMSwarmDataExTopologyInitialize(de);
639: for (r=0; r<size; r++) {
640: if (n2collect[r] > 0) {
641: DMSwarmDataExTopologyAddNeighbour(de,(PetscMPIInt)r);
642: }
643: }
644: DMSwarmDataExTopologyFinalize(de);
645: /* Define send counts */
646: DMSwarmDataExInitializeSendCount(de);
647: for (r=0; r<size; r++) {
648: if (n2collect[r] > 0) {
649: DMSwarmDataExAddToSendCount(de,r,n2collect[r]);
650: }
651: }
652: DMSwarmDataExFinalizeSendCount(de);
653: /* Pack data */
654: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
655: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
656: for (r=0; r<size; r++) {
657: for (p=0; p<n2collect[r]; p++) {
658: DMSwarmDataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);
659: /* insert point buffer into the data exchanger */
660: DMSwarmDataExPackData(de,r,1,point_buffer);
661: }
662: }
663: DMSwarmDataExPackFinalize(de);
664: /* Scatter */
665: DMSwarmDataExBegin(de);
666: DMSwarmDataExEnd(de);
667: /* Collect data in DMSwarm container */
668: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
669: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
670: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
671: for (p=0; p<n_points_recv; p++) {
672: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
674: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
675: }
676: /* Release memory */
677: for (r=0; r<size; r++) {
678: if (collectlist[r]) PetscFree(collectlist[r]);
679: }
680: PetscFree(collectlist);
681: PetscFree(n2collect);
682: PetscFree(ctxlist);
683: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
684: DMSwarmDataExView(de);
685: DMSwarmDataExDestroy(de);
686: return 0;
687: }