xref: /petsc/src/dm/impls/swarm/swarm_migrate.c (revision 2bf68e3e0f2a61f71e7c65bee250bfa1c8ce0cdb) 
1 
2 #include <petscsf.h>
3 #include <petscdmswarm.h>
4 #include <petsc/private/dmswarmimpl.h>    /*I   "petscdmswarm.h"   I*/
5 #include "data_bucket.h"
6 #include "data_ex.h"
7 
8 
9 /*
10  User loads desired location (MPI rank) into field DMSwarm_rank
11 */
12 #undef __FUNCT__
13 #define __FUNCT__ "DMSwarmMigrate_Push_Basic"
14 PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points)
15 {
16   DM_Swarm *swarm = (DM_Swarm*)dm->data;
17   PetscErrorCode ierr;
18   DataEx de;
19   PetscInt p,npoints,*rankval,n_points_recv;
20   PetscMPIInt rank,nrank;
21   void *point_buffer,*recv_points;
22   size_t sizeof_dmswarm_point;
23 
24   PetscFunctionBegin;
25   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
26 
27   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
28   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
29   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0, &de);CHKERRQ(ierr);
30   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
31   for (p = 0; p < npoints; ++p) {
32     nrank = rankval[p];
33     if (nrank != rank) {
34       ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr);
35     }
36   }
37   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
38   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
39   for (p=0; p<npoints; p++) {
40     nrank = rankval[p];
41     if (nrank != rank) {
42       ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr);
43     }
44   }
45   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
46   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
47   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
48   for (p=0; p<npoints; p++) {
49     nrank = rankval[p];
50     if (nrank != rank) {
51       /* copy point into buffer */
52       ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
53       /* insert point buffer into DataExchanger */
54       ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr);
55     }
56   }
57   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
58   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
59 
60   if (remove_sent_points) {
61     DataField gfield;
62 
63     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);CHKERRQ(ierr);
64     ierr = DataFieldGetAccess(gfield);CHKERRQ(ierr);
65     ierr = DataFieldGetEntries(gfield,(void**)&rankval);CHKERRQ(ierr);
66 
67     /* remove points which left processor */
68     ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
69     for (p=0; p<npoints; p++) {
70       nrank = rankval[p];
71       if (nrank != rank) {
72         /* kill point */
73         ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
74 
75         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
76 
77         ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
78         ierr = DataFieldGetAccess(gfield);CHKERRQ(ierr);
79         ierr = DataFieldGetEntries(gfield,(void**)&rankval);CHKERRQ(ierr);
80         p--; /* check replacement point */
81       }
82     }
83     ierr = DataFieldRestoreEntries(gfield,(void**)&rankval);CHKERRQ(ierr);
84     ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
85   }
86   ierr = DataExBegin(de);CHKERRQ(ierr);
87   ierr = DataExEnd(de);CHKERRQ(ierr);
88   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
89   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
90   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
91   for (p=0; p<n_points_recv; p++) {
92     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
93 
94     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
95   }
96   ierr = DataExView(de);CHKERRQ(ierr);
97   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
98   ierr = DataExDestroy(de);CHKERRQ(ierr);
99   PetscFunctionReturn(0);
100 }
101 
102 #undef __FUNCT__
103 #define __FUNCT__ "DMSwarmMigrate_DMNeighborScatter"
104 PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm,DM dmcell,PetscBool remove_sent_points,PetscInt *npoints_prior_migration)
105 {
106   DM_Swarm *swarm = (DM_Swarm*)dm->data;
107   PetscErrorCode ierr;
108   DataEx de;
109   PetscInt r,p,npoints,*rankval,n_points_recv;
110   PetscMPIInt rank,_rank;
111   const PetscMPIInt *neighbourranks;
112   void *point_buffer,*recv_points;
113   size_t sizeof_dmswarm_point;
114   PetscInt nneighbors;
115   PetscMPIInt mynneigh,*myneigh;
116 
117   PetscFunctionBegin;
118   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
119   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
120   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
121   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
122   ierr = DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);CHKERRQ(ierr);
123   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
124   for (r=0; r<nneighbors; r++) {
125     _rank = neighbourranks[r];
126     if ((_rank != rank) && (_rank > 0)) {
127       ierr = DataExTopologyAddNeighbour(de,_rank);CHKERRQ(ierr);
128     }
129   }
130   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
131   ierr = DataExTopologyGetNeighbours(de,&mynneigh,&myneigh);CHKERRQ(ierr);
132   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
133   for (p=0; p<npoints; p++) {
134     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
135       for (r=0; r<mynneigh; r++) {
136         _rank = myneigh[r];
137         ierr = DataExAddToSendCount(de,_rank,1);CHKERRQ(ierr);
138       }
139     }
140   }
141   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
142   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
143   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
144   for (p=0; p<npoints; p++) {
145     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
146       for (r=0; r<mynneigh; r++) {
147         _rank = myneigh[r];
148         /* copy point into buffer */
149         ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
150         /* insert point buffer into DataExchanger */
151         ierr = DataExPackData(de,_rank,1,point_buffer);CHKERRQ(ierr);
152       }
153     }
154   }
155   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
156   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
157   if (remove_sent_points) {
158     DataField PField;
159 
160     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr);
161     ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr);
162     /* remove points which left processor */
163     ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
164     for (p=0; p<npoints; p++) {
165       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
166         /* kill point */
167         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
168         ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
169         ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */
170         p--; /* check replacement point */
171       }
172     }
173   }
174   ierr = DataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);CHKERRQ(ierr);
175   ierr = DataExBegin(de);CHKERRQ(ierr);
176   ierr = DataExEnd(de);CHKERRQ(ierr);
177   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
178   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
179   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
180   for (p=0; p<n_points_recv; p++) {
181     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
182 
183     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
184   }
185   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
186   ierr = DataExDestroy(de);CHKERRQ(ierr);
187   PetscFunctionReturn(0);
188 }
189 
190 #undef __FUNCT__
191 #define __FUNCT__ "DMSwarmMigrate_CellDMScatter"
192 PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points)
193 {
194   DM_Swarm *swarm = (DM_Swarm*)dm->data;
195   PetscErrorCode ierr;
196   PetscInt p,npoints,npointsg=0,npoints2,npoints2g,*rankval,npoints_prior_migration;
197   PetscSF sfcell = NULL;
198   const PetscSFNode *LA_sfcell;
199   DM dmcell;
200   Vec pos;
201   PetscBool error_check = swarm->migrate_error_on_missing_point;
202   PetscMPIInt commsize,rank;
203 
204   PetscFunctionBegin;
205   ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr);
206   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
207 
208   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr);
209   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
210 
211   ierr = DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);CHKERRQ(ierr);
212   ierr = DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);CHKERRQ(ierr);
213   ierr = DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);CHKERRQ(ierr);
214 
215   if (error_check) {
216     ierr = DMSwarmGetSize(dm,&npointsg);CHKERRQ(ierr);
217   }
218   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
219   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
220   ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
221   for (p=0; p<npoints; p++) {
222     rankval[p] = LA_sfcell[p].index;
223   }
224   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
225   ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
226 
227   if (commsize > 1) {
228     ierr = DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);CHKERRQ(ierr);
229   } else {
230     ierr = DMSwarmGetSize(dm,&npoints_prior_migration);CHKERRQ(ierr);
231   }
232 
233   /* locate points newly recevied */
234   ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr);
235 #if 0
236   len = npoints2 - npoints_prior_migration;
237   if (len > 0) {
238     PetscScalar *LA_coor;
239     PetscInt bs;
240 
241     ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
242     ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*len,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);CHKERRQ(ierr);
243 
244     ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr);
245 
246     ierr = VecDestroy(&pos);CHKERRQ(ierr);
247     ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
248 
249     ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr);
250     ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
251     for (p=0; p<len; p++) {
252       rankval[npoints_prior_migration+p] = LA_iscell[p];
253     }
254     ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr);
255     ierr = ISDestroy(&iscell);CHKERRQ(ierr);
256     ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
257 
258     /* remove points which left processor */
259     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr);
260     ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr);
261 
262     ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
263     for (p=npoints_prior_migration; p<npoints2; p++) {
264       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
265         /* kill point */
266         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
267         ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
268         ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */
269         p--; /* check replacement point */
270       }
271     }
272 
273   }
274 #endif
275 
276   {
277     PetscScalar *LA_coor;
278     PetscInt bs;
279     DataField PField;
280 
281     ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
282     ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);CHKERRQ(ierr);
283     ierr = DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr);
284 
285     ierr = VecDestroy(&pos);CHKERRQ(ierr);
286     ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
287 
288     ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
289     ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
290     for (p=0; p<npoints2; p++) {
291       rankval[p] = LA_sfcell[p].index;
292     }
293     ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
294     ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
295 
296     /* remove points which left processor */
297     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr);
298     ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr);
299 
300     ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr);
301     for (p=0; p<npoints2; p++) {
302       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
303         /* kill point */
304         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
305         ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
306         ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */
307         p--; /* check replacement point */
308       }
309     }
310   }
311   /* check for error on removed points */
312   if (error_check) {
313     ierr = DMSwarmGetSize(dm,&npoints2g);CHKERRQ(ierr);
314     if (npointsg != npoints2g) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points from the DMSwarm must remain constant during migration (initial %D - final %D)",npointsg,npoints2g);
315   }
316   PetscFunctionReturn(0);
317 }
318 
319 #undef __FUNCT__
320 #define __FUNCT__ "DMSwarmMigrate_CellDMExact"
321 PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
322 {
323   PetscFunctionBegin;
324   PetscFunctionReturn(0);
325 }
326 
327 /*
328  Redundant as this assumes points can only be sent to a single rank
329 */
330 #undef __FUNCT__
331 #define __FUNCT__ "DMSwarmMigrate_GlobalToLocal_Basic"
332 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
333 {
334   DM_Swarm *swarm = (DM_Swarm*)dm->data;
335   PetscErrorCode ierr;
336   DataEx de;
337   PetscInt p,npoints,*rankval,n_points_recv;
338   PetscMPIInt rank,nrank,negrank;
339   void *point_buffer,*recv_points;
340   size_t sizeof_dmswarm_point;
341 
342   PetscFunctionBegin;
343   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
344   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
345   *globalsize = npoints;
346   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
347   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
348   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
349   for (p=0; p<npoints; p++) {
350     negrank = rankval[p];
351     if (negrank < 0) {
352       nrank = -negrank - 1;
353       ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr);
354     }
355   }
356   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
357   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
358   for (p=0; p<npoints; p++) {
359     negrank = rankval[p];
360     if (negrank < 0) {
361       nrank = -negrank - 1;
362       ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr);
363     }
364   }
365   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
366   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
367   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
368   for (p=0; p<npoints; p++) {
369     negrank = rankval[p];
370     if (negrank < 0) {
371       nrank = -negrank - 1;
372       rankval[p] = nrank;
373       /* copy point into buffer */
374       ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
375       /* insert point buffer into DataExchanger */
376       ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr);
377       rankval[p] = negrank;
378     }
379   }
380   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
381   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
382   ierr = DataExBegin(de);CHKERRQ(ierr);
383   ierr = DataExEnd(de);CHKERRQ(ierr);
384   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
385   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
386   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
387   for (p=0; p<n_points_recv; p++) {
388     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
389 
390     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
391   }
392   ierr = DataExView(de);CHKERRQ(ierr);
393   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
394   ierr = DataExDestroy(de);CHKERRQ(ierr);
395   PetscFunctionReturn(0);
396 }
397 
398 typedef struct {
399   PetscMPIInt owner_rank;
400   PetscReal min[3],max[3];
401 } CollectBBox;
402 
403 #undef __FUNCT__
404 #define __FUNCT__ "DMSwarmCollect_DMDABoundingBox"
405 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
406 {
407   DM_Swarm *swarm = (DM_Swarm*)dm->data;
408   PetscErrorCode ierr;
409   DataEx de;
410   PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
411   PetscMPIInt rank,nrank;
412   void *point_buffer,*recv_points;
413   size_t sizeof_dmswarm_point,sizeof_bbox_ctx;
414   PetscBool isdmda;
415   CollectBBox *bbox,*recv_bbox;
416   const PetscMPIInt *dmneighborranks;
417   DM dmcell;
418 
419   PetscFunctionBegin;
420   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
421 
422   ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr);
423   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
424   isdmda = PETSC_FALSE;
425   PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
426   if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox");
427 
428   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
429   sizeof_bbox_ctx = sizeof(CollectBBox);
430   PetscMalloc1(1,&bbox);
431   bbox->owner_rank = rank;
432 
433   /* compute the bounding box based on the overlapping / stenctil size */
434   {
435     Vec lcoor;
436 
437     ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr);
438     if (dim >= 1) {
439       ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr);
440       ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr);
441     }
442     if (dim >= 2) {
443       ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr);
444       ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr);
445     }
446     if (dim == 3) {
447       ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr);
448       ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr);
449     }
450   }
451   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
452   *globalsize = npoints;
453   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
454   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
455   /* use DMDA neighbours */
456   ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr);
457   if (dim == 1) {
458     neighbour_cells = 3;
459   } else if (dim == 2) {
460     neighbour_cells = 9;
461   } else {
462     neighbour_cells = 27;
463   }
464   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
465   for (p=0; p<neighbour_cells; p++) {
466     if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
467       ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr);
468     }
469   }
470   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
471   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
472   for (p=0; p<neighbour_cells; p++) {
473     if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
474       ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr);
475     }
476   }
477   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
478   /* send bounding boxes */
479   ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr);
480   for (p=0; p<neighbour_cells; p++) {
481     nrank = dmneighborranks[p];
482     if ( (nrank >= 0) && (nrank != rank) ) {
483       /* insert bbox buffer into DataExchanger */
484       ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr);
485     }
486   }
487   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
488   /* recv bounding boxes */
489   ierr = DataExBegin(de);CHKERRQ(ierr);
490   ierr = DataExEnd(de);CHKERRQ(ierr);
491   ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr);
492   for (p=0; p<n_bbox_recv; p++) {
493     PetscPrintf(PETSC_COMM_SELF,"[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank,
494            (double)recv_bbox[p].min[0],(double)recv_bbox[p].max[0],(double)recv_bbox[p].min[1],(double)recv_bbox[p].max[1]);
495   }
496   /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
497   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
498   for (pk=0; pk<n_bbox_recv; pk++) {
499     PetscReal *array_x,*array_y;
500 
501     ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
502     ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
503     for (p=0; p<npoints; p++) {
504       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
505         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
506           ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr);
507         }
508       }
509     }
510     ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
511     ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
512   }
513   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
514   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
515   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
516   for (pk=0; pk<n_bbox_recv; pk++) {
517     PetscReal *array_x,*array_y;
518 
519     ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
520     ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
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]) ) {
524           /* copy point into buffer */
525           ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
526           /* insert point buffer into DataExchanger */
527           ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr);
528         }
529       }
530     }
531     ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
532     ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
533   }
534   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
535   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
536   ierr = DataExBegin(de);CHKERRQ(ierr);
537   ierr = DataExEnd(de);CHKERRQ(ierr);
538   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
539   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
540   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
541   for (p=0; p<n_points_recv; p++) {
542     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
543 
544     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
545   }
546   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
547   PetscFree(bbox);
548   ierr = DataExView(de);CHKERRQ(ierr);
549   ierr = DataExDestroy(de);CHKERRQ(ierr);
550   PetscFunctionReturn(0);
551 }
552 
553 
554 /* General collection when no order, or neighbour information is provided */
555 /*
556  User provides context and collect() method
557  Broadcast user context
558 
559  for each context / rank {
560    collect(swarm,context,n,list)
561  }
562 */
563 #undef __FUNCT__
564 #define __FUNCT__ "DMSwarmCollect_General"
565 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
566 {
567   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
568   PetscErrorCode ierr;
569   DataEx         de;
570   PetscInt       p,r,npoints,n_points_recv;
571   PetscMPIInt    commsize,rank;
572   void           *point_buffer,*recv_points;
573   void           *ctxlist;
574   PetscInt       *n2collect,**collectlist;
575   size_t         sizeof_dmswarm_point;
576 
577   PetscFunctionBegin;
578   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr);
579   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
580   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
581   *globalsize = npoints;
582   /* Broadcast user context */
583   PetscMalloc(ctx_size*commsize,&ctxlist);
584   ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
585   ierr = PetscMalloc1(commsize,&n2collect);CHKERRQ(ierr);
586   ierr = PetscMalloc1(commsize,&collectlist);CHKERRQ(ierr);
587   for (r=0; r<commsize; r++) {
588     PetscInt _n2collect;
589     PetscInt *_collectlist;
590     void     *_ctx_r;
591 
592     _n2collect   = 0;
593     _collectlist = NULL;
594     if (r != rank) { /* don't collect data from yourself */
595       _ctx_r = (void*)( (char*)ctxlist + r * ctx_size );
596       ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr);
597     }
598     n2collect[r]   = _n2collect;
599     collectlist[r] = _collectlist;
600   }
601   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
602   /* Define topology */
603   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
604   for (r=0; r<commsize; r++) {
605     if (n2collect[r] > 0) {
606       ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr);
607     }
608   }
609   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
610   /* Define send counts */
611   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
612   for (r=0; r<commsize; r++) {
613     if (n2collect[r] > 0) {
614       ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr);
615     }
616   }
617   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
618   /* Pack data */
619   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
620   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
621   for (r=0; r<commsize; r++) {
622     for (p=0; p<n2collect[r]; p++) {
623       ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr);
624       /* insert point buffer into the data exchanger */
625       ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr);
626     }
627   }
628   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
629   /* Scatter */
630   ierr = DataExBegin(de);CHKERRQ(ierr);
631   ierr = DataExEnd(de);CHKERRQ(ierr);
632   /* Collect data in DMSwarm container */
633   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
634   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
635   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
636   for (p=0; p<n_points_recv; p++) {
637     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
638 
639     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
640   }
641   /* Release memory */
642   for (r=0; r<commsize; r++) {
643     if (collectlist[r]) PetscFree(collectlist[r]);
644   }
645   ierr = PetscFree(collectlist);CHKERRQ(ierr);
646   ierr = PetscFree(n2collect);CHKERRQ(ierr);
647   ierr = PetscFree(ctxlist);CHKERRQ(ierr);
648   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
649   ierr = DataExView(de);CHKERRQ(ierr);
650   ierr = DataExDestroy(de);CHKERRQ(ierr);
651   PetscFunctionReturn(0);
652 }
653 
654