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) DMSwarmDataExTopologyAddNeighbour(de, nrank);
 30:   }
 31:   DMSwarmDataExTopologyFinalize(de);
 32:   DMSwarmDataExInitializeSendCount(de);
 33:   for (p = 0; p < npoints; p++) {
 34:     nrank = rankval[p];
 35:     if (nrank != rank) DMSwarmDataExAddToSendCount(de, nrank, 1);
 36:   }
 37:   DMSwarmDataExFinalizeSendCount(de);
 38:   DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
 39:   DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
 40:   for (p = 0; p < npoints; p++) {
 41:     nrank = rankval[p];
 42:     if (nrank != rank) {
 43:       /* copy point into buffer */
 44:       DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
 45:       /* insert point buffer into DMSwarmDataExchanger */
 46:       DMSwarmDataExPackData(de, nrank, 1, point_buffer);
 47:     }
 48:   }
 49:   DMSwarmDataExPackFinalize(de);
 50:   DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);

 52:   if (remove_sent_points) {
 53:     DMSwarmDataField gfield;

 55:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &gfield);
 56:     DMSwarmDataFieldGetAccess(gfield);
 57:     DMSwarmDataFieldGetEntries(gfield, (void **)&rankval);

 59:     /* remove points which left processor */
 60:     DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
 61:     for (p = 0; p < npoints; p++) {
 62:       nrank = rankval[p];
 63:       if (nrank != rank) {
 64:         /* kill point */
 65:         DMSwarmDataFieldRestoreAccess(gfield);

 67:         DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);

 69:         DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL); /* you need to update npoints as the list size decreases! */
 70:         DMSwarmDataFieldGetAccess(gfield);
 71:         DMSwarmDataFieldGetEntries(gfield, (void **)&rankval);
 72:         p--; /* check replacement point */
 73:       }
 74:     }
 75:     DMSwarmDataFieldRestoreEntries(gfield, (void **)&rankval);
 76:     DMSwarmDataFieldRestoreAccess(gfield);
 77:   }
 78:   DMSwarmDataExBegin(de);
 79:   DMSwarmDataExEnd(de);
 80:   DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
 81:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
 82:   DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
 83:   for (p = 0; p < n_points_recv; p++) {
 84:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

 86:     DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
 87:   }
 88:   if (debug) DMSwarmDataExView(de);
 89:   DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
 90:   DMSwarmDataExDestroy(de);
 91:   return 0;
 92: }

 94: PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm, DM dmcell, PetscBool remove_sent_points, PetscInt *npoints_prior_migration)
 95: {
 96:   DM_Swarm          *swarm = (DM_Swarm *)dm->data;
 97:   DMSwarmDataEx      de;
 98:   PetscInt           r, p, npoints, *rankval, n_points_recv;
 99:   PetscMPIInt        rank, _rank;
100:   const PetscMPIInt *neighbourranks;
101:   void              *point_buffer, *recv_points;
102:   size_t             sizeof_dmswarm_point;
103:   PetscInt           nneighbors;
104:   PetscMPIInt        mynneigh, *myneigh;

106:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
107:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
108:   DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
109:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
110:   DMGetNeighbors(dmcell, &nneighbors, &neighbourranks);
111:   DMSwarmDataExTopologyInitialize(de);
112:   for (r = 0; r < nneighbors; r++) {
113:     _rank = neighbourranks[r];
114:     if ((_rank != rank) && (_rank > 0)) DMSwarmDataExTopologyAddNeighbour(de, _rank);
115:   }
116:   DMSwarmDataExTopologyFinalize(de);
117:   DMSwarmDataExTopologyGetNeighbours(de, &mynneigh, &myneigh);
118:   DMSwarmDataExInitializeSendCount(de);
119:   for (p = 0; p < npoints; p++) {
120:     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
121:       for (r = 0; r < mynneigh; r++) {
122:         _rank = myneigh[r];
123:         DMSwarmDataExAddToSendCount(de, _rank, 1);
124:       }
125:     }
126:   }
127:   DMSwarmDataExFinalizeSendCount(de);
128:   DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
129:   DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
130:   for (p = 0; p < npoints; p++) {
131:     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
132:       for (r = 0; r < mynneigh; r++) {
133:         _rank = myneigh[r];
134:         /* copy point into buffer */
135:         DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
136:         /* insert point buffer into DMSwarmDataExchanger */
137:         DMSwarmDataExPackData(de, _rank, 1, point_buffer);
138:       }
139:     }
140:   }
141:   DMSwarmDataExPackFinalize(de);
142:   DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
143:   if (remove_sent_points) {
144:     DMSwarmDataField PField;

146:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &PField);
147:     DMSwarmDataFieldGetEntries(PField, (void **)&rankval);
148:     /* remove points which left processor */
149:     DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
150:     for (p = 0; p < npoints; p++) {
151:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
152:         /* kill point */
153:         DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);
154:         DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL); /* you need to update npoints as the list size decreases! */
155:         DMSwarmDataFieldGetEntries(PField, (void **)&rankval);      /* update date point increase realloc performed */
156:         p--;                                                                   /* check replacement point */
157:       }
158:     }
159:   }
160:   DMSwarmDataBucketGetSizes(swarm->db, npoints_prior_migration, NULL, NULL);
161:   DMSwarmDataExBegin(de);
162:   DMSwarmDataExEnd(de);
163:   DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
164:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
165:   DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
166:   for (p = 0; p < n_points_recv; p++) {
167:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

169:     DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
170:   }
171:   DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
172:   DMSwarmDataExDestroy(de);
173:   return 0;
174: }

176: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm, PetscBool remove_sent_points)
177: {
178:   DM_Swarm          *swarm = (DM_Swarm *)dm->data;
179:   PetscInt           p, npoints, npointsg = 0, npoints2, npoints2g, *rankval, npoints_prior_migration;
180:   PetscSF            sfcell = NULL;
181:   const PetscSFNode *LA_sfcell;
182:   DM                 dmcell;
183:   Vec                pos;
184:   PetscBool          error_check = swarm->migrate_error_on_missing_point;
185:   PetscMPIInt        size, rank;

187:   DMSwarmGetCellDM(dm, &dmcell);

190:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
191:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);

193: #if 1
194:   {
195:     PetscInt    *p_cellid;
196:     PetscInt     npoints_curr, range = 0;
197:     PetscSFNode *sf_cells;

199:     DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL);
200:     PetscMalloc1(npoints_curr, &sf_cells);

202:     DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
203:     DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
204:     for (p = 0; p < npoints_curr; p++) {
205:       sf_cells[p].rank  = 0;
206:       sf_cells[p].index = p_cellid[p];
207:       if (p_cellid[p] > range) range = p_cellid[p];
208:     }
209:     DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
210:     DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);

212:     /* PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell); */
213:     PetscSFCreate(PETSC_COMM_SELF, &sfcell);
214:     PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
215:   }
216: #endif

218:   DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
219:   DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
220:   DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);

222:   if (error_check) DMSwarmGetSize(dm, &npointsg);
223:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
224:   DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
225:   PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
226:   for (p = 0; p < npoints; p++) rankval[p] = LA_sfcell[p].index;
227:   DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
228:   PetscSFDestroy(&sfcell);

230:   if (size > 1) {
231:     DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration);
232:   } else {
233:     DMSwarmDataField PField;
234:     PetscInt         npoints_curr;

236:     /* remove points which the domain */
237:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &PField);
238:     DMSwarmDataFieldGetEntries(PField, (void **)&rankval);

240:     DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL);
241:     for (p = 0; p < npoints_curr; p++) {
242:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
243:         /* kill point */
244:         DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);
245:         DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL); /* you need to update npoints as the list size decreases! */
246:         DMSwarmDataFieldGetEntries(PField, (void **)&rankval);           /* update date point increase realloc performed */
247:         p--;                                                                        /* check replacement point */
248:       }
249:     }
250:     DMSwarmGetSize(dm, &npoints_prior_migration);
251:   }

253:   /* locate points newly received */
254:   DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL);

256: #if 0
257:   { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
258:     PetscScalar *LA_coor;
259:     PetscInt bs;
260:     DMSwarmDataField PField;

262:     DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
263:     VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);
264:     DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);

266:     VecDestroy(&pos);
267:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);

269:     PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
270:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
271:     for (p=0; p<npoints2; p++) {
272:       rankval[p] = LA_sfcell[p].index;
273:     }
274:     PetscSFDestroy(&sfcell);
275:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);

277:     /* remove points which left processor */
278:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
279:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);

281:     DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
282:     for (p=0; p<npoints2; p++) {
283:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
284:         /* kill point */
285:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
286:         DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
287:         DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
288:         p--; /* check replacement point */
289:       }
290:     }
291:   }
292: #endif

294:   { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
295:     PetscScalar     *LA_coor;
296:     PetscInt         npoints_from_neighbours, bs;
297:     DMSwarmDataField PField;

299:     npoints_from_neighbours = npoints2 - npoints_prior_migration;

301:     DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor);
302:     VecCreateSeqWithArray(PETSC_COMM_SELF, bs, bs * npoints_from_neighbours, (const PetscScalar *)&LA_coor[bs * npoints_prior_migration], &pos);

304:     DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);

306:     VecDestroy(&pos);
307:     DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor);

309:     PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
310:     DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
311:     for (p = 0; p < npoints_from_neighbours; p++) rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
312:     DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
313:     PetscSFDestroy(&sfcell);

315:     /* remove points which left processor */
316:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &PField);
317:     DMSwarmDataFieldGetEntries(PField, (void **)&rankval);

319:     DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL);
320:     for (p = npoints_prior_migration; p < npoints2; p++) {
321:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
322:         /* kill point */
323:         DMSwarmDataBucketRemovePointAtIndex(swarm->db, p);
324:         DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL); /* you need to update npoints as the list size decreases! */
325:         DMSwarmDataFieldGetEntries(PField, (void **)&rankval);       /* update date point increase realloc performed */
326:         p--;                                                                    /* check replacement point */
327:       }
328:     }
329:   }

331:   {
332:     PetscInt *p_cellid;

334:     DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL);
335:     DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
336:     DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
337:     for (p = 0; p < npoints2; p++) p_cellid[p] = rankval[p];
338:     DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid);
339:     DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
340:   }

342:   /* check for error on removed points */
343:   if (error_check) {
344:     DMSwarmGetSize(dm, &npoints2g);
346:   }
347:   return 0;
348: }

350: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm, PetscBool remove_sent_points)
351: {
352:   return 0;
353: }

355: /*
356:  Redundant as this assumes points can only be sent to a single rank
357: */
358: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize)
359: {
360:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
361:   DMSwarmDataEx de;
362:   PetscInt      p, npoints, *rankval, n_points_recv;
363:   PetscMPIInt   rank, nrank, negrank;
364:   void         *point_buffer, *recv_points;
365:   size_t        sizeof_dmswarm_point;

367:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
368:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
369:   *globalsize = npoints;
370:   DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
371:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
372:   DMSwarmDataExTopologyInitialize(de);
373:   for (p = 0; p < npoints; p++) {
374:     negrank = rankval[p];
375:     if (negrank < 0) {
376:       nrank = -negrank - 1;
377:       DMSwarmDataExTopologyAddNeighbour(de, nrank);
378:     }
379:   }
380:   DMSwarmDataExTopologyFinalize(de);
381:   DMSwarmDataExInitializeSendCount(de);
382:   for (p = 0; p < npoints; p++) {
383:     negrank = rankval[p];
384:     if (negrank < 0) {
385:       nrank = -negrank - 1;
386:       DMSwarmDataExAddToSendCount(de, nrank, 1);
387:     }
388:   }
389:   DMSwarmDataExFinalizeSendCount(de);
390:   DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
391:   DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
392:   for (p = 0; p < npoints; p++) {
393:     negrank = rankval[p];
394:     if (negrank < 0) {
395:       nrank      = -negrank - 1;
396:       rankval[p] = nrank;
397:       /* copy point into buffer */
398:       DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
399:       /* insert point buffer into DMSwarmDataExchanger */
400:       DMSwarmDataExPackData(de, nrank, 1, point_buffer);
401:       rankval[p] = negrank;
402:     }
403:   }
404:   DMSwarmDataExPackFinalize(de);
405:   DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
406:   DMSwarmDataExBegin(de);
407:   DMSwarmDataExEnd(de);
408:   DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
409:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
410:   DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
411:   for (p = 0; p < n_points_recv; p++) {
412:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

414:     DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
415:   }
416:   DMSwarmDataExView(de);
417:   DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
418:   DMSwarmDataExDestroy(de);
419:   return 0;
420: }

422: typedef struct {
423:   PetscMPIInt owner_rank;
424:   PetscReal   min[3], max[3];
425: } CollectBBox;

427: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm, PetscInt *globalsize)
428: {
429:   DM_Swarm          *swarm = (DM_Swarm *)dm->data;
430:   DMSwarmDataEx      de;
431:   PetscInt           p, pk, npoints, *rankval, n_points_recv, n_bbox_recv, dim, neighbour_cells;
432:   PetscMPIInt        rank, nrank;
433:   void              *point_buffer, *recv_points;
434:   size_t             sizeof_dmswarm_point, sizeof_bbox_ctx;
435:   PetscBool          isdmda;
436:   CollectBBox       *bbox, *recv_bbox;
437:   const PetscMPIInt *dmneighborranks;
438:   DM                 dmcell;

440:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);

442:   DMSwarmGetCellDM(dm, &dmcell);
444:   isdmda = PETSC_FALSE;
445:   PetscObjectTypeCompare((PetscObject)dmcell, DMDA, &isdmda);

448:   DMGetDimension(dm, &dim);
449:   sizeof_bbox_ctx = sizeof(CollectBBox);
450:   PetscMalloc1(1, &bbox);
451:   bbox->owner_rank = rank;

453:   /* compute the bounding box based on the overlapping / stenctil size */
454:   {
455:     Vec lcoor;

457:     DMGetCoordinatesLocal(dmcell, &lcoor);
458:     if (dim >= 1) {
459:       VecStrideMin(lcoor, 0, NULL, &bbox->min[0]);
460:       VecStrideMax(lcoor, 0, NULL, &bbox->max[0]);
461:     }
462:     if (dim >= 2) {
463:       VecStrideMin(lcoor, 1, NULL, &bbox->min[1]);
464:       VecStrideMax(lcoor, 1, NULL, &bbox->max[1]);
465:     }
466:     if (dim == 3) {
467:       VecStrideMin(lcoor, 2, NULL, &bbox->min[2]);
468:       VecStrideMax(lcoor, 2, NULL, &bbox->max[2]);
469:     }
470:   }
471:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
472:   *globalsize = npoints;
473:   DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
474:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
475:   /* use DMDA neighbours */
476:   DMDAGetNeighbors(dmcell, &dmneighborranks);
477:   if (dim == 1) {
478:     neighbour_cells = 3;
479:   } else if (dim == 2) {
480:     neighbour_cells = 9;
481:   } else {
482:     neighbour_cells = 27;
483:   }
484:   DMSwarmDataExTopologyInitialize(de);
485:   for (p = 0; p < neighbour_cells; p++) {
486:     if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) DMSwarmDataExTopologyAddNeighbour(de, dmneighborranks[p]);
487:   }
488:   DMSwarmDataExTopologyFinalize(de);
489:   DMSwarmDataExInitializeSendCount(de);
490:   for (p = 0; p < neighbour_cells; p++) {
491:     if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) DMSwarmDataExAddToSendCount(de, dmneighborranks[p], 1);
492:   }
493:   DMSwarmDataExFinalizeSendCount(de);
494:   /* send bounding boxes */
495:   DMSwarmDataExPackInitialize(de, sizeof_bbox_ctx);
496:   for (p = 0; p < neighbour_cells; p++) {
497:     nrank = dmneighborranks[p];
498:     if ((nrank >= 0) && (nrank != rank)) {
499:       /* insert bbox buffer into DMSwarmDataExchanger */
500:       DMSwarmDataExPackData(de, nrank, 1, bbox);
501:     }
502:   }
503:   DMSwarmDataExPackFinalize(de);
504:   /* recv bounding boxes */
505:   DMSwarmDataExBegin(de);
506:   DMSwarmDataExEnd(de);
507:   DMSwarmDataExGetRecvData(de, &n_bbox_recv, (void **)&recv_bbox);
508:   /*  Wrong, should not be using PETSC_COMM_WORLD */
509:   for (p = 0; p < n_bbox_recv; p++) {
510:     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],
511:                                       (double)recv_bbox[p].max[1]));
512:   }
513:   PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout);
514:   /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
515:   DMSwarmDataExInitializeSendCount(de);
516:   for (pk = 0; pk < n_bbox_recv; pk++) {
517:     PetscReal *array_x, *array_y;

519:     DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x);
520:     DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y);
521:     for (p = 0; p < npoints; p++) {
522:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
523:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) DMSwarmDataExAddToSendCount(de, recv_bbox[pk].owner_rank, 1);
524:       }
525:     }
526:     DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y);
527:     DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x);
528:   }
529:   DMSwarmDataExFinalizeSendCount(de);
530:   DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
531:   DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
532:   for (pk = 0; pk < n_bbox_recv; pk++) {
533:     PetscReal *array_x, *array_y;

535:     DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x);
536:     DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y);
537:     for (p = 0; p < npoints; p++) {
538:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
539:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
540:           /* copy point into buffer */
541:           DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer);
542:           /* insert point buffer into DMSwarmDataExchanger */
543:           DMSwarmDataExPackData(de, recv_bbox[pk].owner_rank, 1, point_buffer);
544:         }
545:       }
546:     }
547:     DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y);
548:     DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x);
549:   }
550:   DMSwarmDataExPackFinalize(de);
551:   DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
552:   DMSwarmDataExBegin(de);
553:   DMSwarmDataExEnd(de);
554:   DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
555:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
556:   DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
557:   for (p = 0; p < n_points_recv; p++) {
558:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

560:     DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
561:   }
562:   DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
563:   PetscFree(bbox);
564:   DMSwarmDataExView(de);
565:   DMSwarmDataExDestroy(de);
566:   return 0;
567: }

569: /* General collection when no order, or neighbour information is provided */
570: /*
571:  User provides context and collect() method
572:  Broadcast user context

574:  for each context / rank {
575:    collect(swarm,context,n,list)
576:  }
577: */
578: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm, PetscErrorCode (*collect)(DM, void *, PetscInt *, PetscInt **), size_t ctx_size, void *ctx, PetscInt *globalsize)
579: {
580:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
581:   DMSwarmDataEx de;
582:   PetscInt      p, r, npoints, n_points_recv;
583:   PetscMPIInt   size, rank;
584:   void         *point_buffer, *recv_points;
585:   void         *ctxlist;
586:   PetscInt     *n2collect, **collectlist;
587:   size_t        sizeof_dmswarm_point;

589:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
590:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
591:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
592:   *globalsize = npoints;
593:   /* Broadcast user context */
594:   PetscMalloc(ctx_size * size, &ctxlist);
595:   MPI_Allgather(ctx, ctx_size, MPI_CHAR, ctxlist, ctx_size, MPI_CHAR, PetscObjectComm((PetscObject)dm));
596:   PetscMalloc1(size, &n2collect);
597:   PetscMalloc1(size, &collectlist);
598:   for (r = 0; r < size; r++) {
599:     PetscInt  _n2collect;
600:     PetscInt *_collectlist;
601:     void     *_ctx_r;

603:     _n2collect   = 0;
604:     _collectlist = NULL;
605:     if (r != rank) { /* don't collect data from yourself */
606:       _ctx_r = (void *)((char *)ctxlist + r * ctx_size);
607:       collect(dm, _ctx_r, &_n2collect, &_collectlist);
608:     }
609:     n2collect[r]   = _n2collect;
610:     collectlist[r] = _collectlist;
611:   }
612:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de);
613:   /* Define topology */
614:   DMSwarmDataExTopologyInitialize(de);
615:   for (r = 0; r < size; r++) {
616:     if (n2collect[r] > 0) DMSwarmDataExTopologyAddNeighbour(de, (PetscMPIInt)r);
617:   }
618:   DMSwarmDataExTopologyFinalize(de);
619:   /* Define send counts */
620:   DMSwarmDataExInitializeSendCount(de);
621:   for (r = 0; r < size; r++) {
622:     if (n2collect[r] > 0) DMSwarmDataExAddToSendCount(de, r, n2collect[r]);
623:   }
624:   DMSwarmDataExFinalizeSendCount(de);
625:   /* Pack data */
626:   DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer);
627:   DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point);
628:   for (r = 0; r < size; r++) {
629:     for (p = 0; p < n2collect[r]; p++) {
630:       DMSwarmDataBucketFillPackedArray(swarm->db, collectlist[r][p], point_buffer);
631:       /* insert point buffer into the data exchanger */
632:       DMSwarmDataExPackData(de, r, 1, point_buffer);
633:     }
634:   }
635:   DMSwarmDataExPackFinalize(de);
636:   /* Scatter */
637:   DMSwarmDataExBegin(de);
638:   DMSwarmDataExEnd(de);
639:   /* Collect data in DMSwarm container */
640:   DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points);
641:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
642:   DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
643:   for (p = 0; p < n_points_recv; p++) {
644:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

646:     DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p);
647:   }
648:   /* Release memory */
649:   for (r = 0; r < size; r++) {
650:     if (collectlist[r]) PetscFree(collectlist[r]);
651:   }
652:   PetscFree(collectlist);
653:   PetscFree(n2collect);
654:   PetscFree(ctxlist);
655:   DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer);
656:   DMSwarmDataExView(de);
657:   DMSwarmDataExDestroy(de);
658:   return 0;
659: }