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