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