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