1 2 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 3 #include "data_bucket.h" 4 #include "data_ex.h" 5 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "DMSwarmMigrate_Push_Basic" 9 PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points) 10 { 11 DM_Swarm *swarm = (DM_Swarm*)dm->data; 12 PetscErrorCode ierr; 13 DataEx de; 14 PetscInt p,npoints,*rankval,n_points_recv; 15 PetscMPIInt rank,nrank; 16 void *point_buffer,*recv_points; 17 size_t sizeof_dmswarm_point; 18 19 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 20 21 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 22 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 23 24 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 25 26 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 27 for (p=0; p<npoints; p++) { 28 nrank = rankval[p]; 29 if (nrank != rank) { 30 ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 31 } 32 } 33 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 34 35 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 36 for (p=0; p<npoints; p++) { 37 nrank = rankval[p]; 38 if (nrank != rank) { 39 ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 40 } 41 } 42 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 43 44 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 45 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 46 for (p=0; p<npoints; p++) { 47 nrank = rankval[p]; 48 if (nrank != rank) { 49 /* copy point into buffer */ 50 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 51 /* insert point buffer into DataExchanger */ 52 ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 53 } 54 } 55 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 56 57 58 if (remove_sent_points) { 59 /* remove points which left processor */ 60 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 61 for (p=0; p<npoints; p++) { 62 nrank = rankval[p]; 63 if (nrank != rank) { 64 /* kill point */ 65 ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 66 DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 67 p--; /* check replacement point */ 68 } 69 } 70 } 71 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 72 73 ierr = DataExBegin(de);CHKERRQ(ierr); 74 ierr = DataExEnd(de);CHKERRQ(ierr); 75 76 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 77 78 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 79 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 80 for (p=0; p<n_points_recv; p++) { 81 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 82 83 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 84 } 85 ierr = DataExView(de);CHKERRQ(ierr); 86 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 87 ierr = DataExDestroy(de);CHKERRQ(ierr); 88 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "DMSwarmMigrate_GlobalToLocal_Basic" 94 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize) 95 { 96 DM_Swarm *swarm = (DM_Swarm*)dm->data; 97 PetscErrorCode ierr; 98 DataEx de; 99 PetscInt p,npoints,*rankval,n_points_recv; 100 PetscMPIInt rank,nrank,negrank; 101 void *point_buffer,*recv_points; 102 size_t sizeof_dmswarm_point; 103 104 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 105 106 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 107 *globalsize = npoints; 108 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 109 110 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 111 112 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 113 for (p=0; p<npoints; p++) { 114 negrank = rankval[p]; 115 if (negrank < 0) { 116 nrank = -negrank - 1; 117 ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 118 } 119 } 120 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 121 122 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 123 for (p=0; p<npoints; p++) { 124 negrank = rankval[p]; 125 if (negrank < 0) { 126 nrank = -negrank - 1; 127 ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 128 } 129 } 130 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 131 132 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 133 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 134 for (p=0; p<npoints; p++) { 135 negrank = rankval[p]; 136 if (negrank < 0) { 137 nrank = -negrank - 1; 138 rankval[p] = nrank; 139 /* copy point into buffer */ 140 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 141 /* insert point buffer into DataExchanger */ 142 ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 143 rankval[p] = negrank; 144 } 145 } 146 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 147 148 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 149 150 ierr = DataExBegin(de);CHKERRQ(ierr); 151 ierr = DataExEnd(de);CHKERRQ(ierr); 152 153 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 154 155 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 156 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 157 for (p=0; p<n_points_recv; p++) { 158 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 159 160 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 161 } 162 ierr = DataExView(de);CHKERRQ(ierr); 163 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 164 ierr = DataExDestroy(de);CHKERRQ(ierr); 165 166 PetscFunctionReturn(0); 167 } 168 169 typedef struct { 170 PetscMPIInt owner_rank; 171 PetscReal min[3],max[3]; 172 } CollectBBox; 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "DMSwarmMigrate_GlobalToLocal_BoundingBox" 176 PETSC_EXTERN PetscErrorCode DMSwarmMigrate_GlobalToLocal_BoundingBox(DM dm,PetscInt *globalsize) 177 { 178 DM_Swarm *swarm = (DM_Swarm*)dm->data; 179 PetscErrorCode ierr; 180 DataEx de; 181 PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells; 182 PetscMPIInt rank,nrank; 183 void *point_buffer,*recv_points; 184 size_t sizeof_dmswarm_point,sizeof_bbox_ctx; 185 PetscBool isdmda; 186 CollectBBox *bbox,*recv_bbox; 187 const PetscMPIInt *dmneighborranks; 188 DM dmcell; 189 190 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 191 192 //ierr = DMSwarmGetDMCell(swarm,&dmcell);CHKERRQ(ierr); 193 dmcell = swarm->dmcell; 194 195 dim = 2; 196 if (dim == 2) { 197 neighbour_cells = 9; 198 } 199 200 isdmda = PETSC_FALSE; 201 PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda); 202 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox"); 203 204 sizeof_bbox_ctx = sizeof(CollectBBox); 205 PetscMalloc1(1,&bbox); 206 bbox->owner_rank = rank; 207 208 /* compute the bounding box based on the overlapping / stenctil size */ 209 //ierr = DMDAGetLocalBoundingBox(dmcell,bbox->min,bbox->max);CHKERRQ(ierr); 210 { 211 Vec lcoor; 212 213 ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr); 214 215 ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr); 216 ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr); 217 ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr); 218 ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr); 219 } 220 221 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 222 *globalsize = npoints; 223 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 224 225 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 226 227 /* use DMDA neighbours */ 228 ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr); 229 230 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 231 for (p=0; p<neighbour_cells; p++) { 232 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 233 ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr); 234 } 235 } 236 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 237 238 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 239 for (p=0; p<neighbour_cells; p++) { 240 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 241 ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr); 242 } 243 } 244 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 245 246 247 /* send bounding boxes */ 248 ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr); 249 for (p=0; p<neighbour_cells; p++) { 250 nrank = dmneighborranks[p]; 251 if ( (nrank >= 0) && (nrank != rank) ) { 252 /* insert bbox buffer into DataExchanger */ 253 ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr); 254 } 255 } 256 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 257 258 /* recv bounding boxes */ 259 ierr = DataExBegin(de);CHKERRQ(ierr); 260 ierr = DataExEnd(de);CHKERRQ(ierr); 261 262 ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr); 263 264 265 for (p=0; p<n_bbox_recv; p++) { 266 printf("[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank, 267 recv_bbox[p].min[0],recv_bbox[p].max[0],recv_bbox[p].min[1],recv_bbox[p].max[1]); 268 } 269 270 /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */ 271 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 272 for (pk=0; pk<n_bbox_recv; pk++) { 273 PetscReal *array_x,*array_y; 274 275 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 276 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 277 278 for (p=0; p<npoints; p++) { 279 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 280 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 281 282 ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr); 283 284 } 285 } 286 } 287 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 288 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 289 } 290 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 291 292 293 294 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 295 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 296 for (pk=0; pk<n_bbox_recv; pk++) { 297 PetscReal *array_x,*array_y; 298 299 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 300 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 301 302 for (p=0; p<npoints; p++) { 303 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 304 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 305 // copy point into buffer // 306 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 307 // insert point buffer into DataExchanger // 308 ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr); 309 } 310 } 311 } 312 313 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 314 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 315 } 316 317 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 318 319 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 320 321 ierr = DataExBegin(de);CHKERRQ(ierr); 322 ierr = DataExEnd(de);CHKERRQ(ierr); 323 324 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 325 326 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 327 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 328 for (p=0; p<n_points_recv; p++) { 329 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 330 331 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 332 } 333 334 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 335 PetscFree(bbox); 336 ierr = DataExView(de);CHKERRQ(ierr); 337 ierr = DataExDestroy(de);CHKERRQ(ierr); 338 339 PetscFunctionReturn(0); 340 } 341 342 343 /* General collection when no order, or neighbour information is provided */ 344 /* 345 User provides context and collect() method 346 Broadcast user context 347 348 for each context / rank { 349 collect(swarm,context,n,list) 350 } 351 352 */ 353 #undef __FUNCT__ 354 #define __FUNCT__ "DMSwarmCollect_General" 355 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize) 356 { 357 DM_Swarm *swarm = (DM_Swarm*)dm->data; 358 PetscErrorCode ierr; 359 DataEx de; 360 PetscInt p,r,npoints,n_points_recv; 361 PetscMPIInt commsize,rank; 362 void *point_buffer,*recv_points; 363 void *ctxlist; 364 PetscInt *n2collect,**collectlist; 365 size_t sizeof_dmswarm_point; 366 367 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr); 368 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 369 370 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 371 *globalsize = npoints; 372 373 /* Broadcast user context */ 374 PetscMalloc(ctx_size*commsize,&ctxlist); 375 ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 376 377 PetscMalloc1(commsize,&n2collect); 378 PetscMalloc1(commsize,&collectlist); 379 380 for (r=0; r<commsize; r++) { 381 PetscInt _n2collect; 382 PetscInt *_collectlist; 383 void *_ctx_r; 384 385 _n2collect = 0; 386 _collectlist = NULL; 387 388 if (r != rank) { /* don't collect data from yourself */ 389 _ctx_r = (void*)( (char*)ctxlist + r * ctx_size ); 390 ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr); 391 } 392 393 n2collect[r] = _n2collect; 394 collectlist[r] = _collectlist; 395 } 396 397 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 398 399 /* Define topology */ 400 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 401 for (r=0; r<commsize; r++) { 402 if (n2collect[r] > 0) { 403 ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr); 404 } 405 } 406 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 407 408 /* Define send counts */ 409 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 410 for (r=0; r<commsize; r++) { 411 if (n2collect[r] > 0) { 412 ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr); 413 } 414 } 415 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 416 417 /* Pack data */ 418 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 419 420 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 421 for (r=0; r<commsize; r++) { 422 423 for (p=0; p<n2collect[r]; p++) { 424 ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr); 425 /* insert point buffer into the data exchanger */ 426 ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr); 427 } 428 } 429 430 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 431 432 /* Scatter */ 433 ierr = DataExBegin(de);CHKERRQ(ierr); 434 ierr = DataExEnd(de);CHKERRQ(ierr); 435 436 /* Collect data in DMSwarm container */ 437 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 438 439 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 440 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 441 for (p=0; p<n_points_recv; p++) { 442 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 443 444 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 445 } 446 447 /* Release memory */ 448 for (r=0; r<commsize; r++) { 449 if (collectlist[r]) PetscFree(collectlist[r]); 450 } 451 PetscFree(collectlist); 452 PetscFree(n2collect); 453 PetscFree(ctxlist); 454 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 455 ierr = DataExView(de);CHKERRQ(ierr); 456 ierr = DataExDestroy(de);CHKERRQ(ierr); 457 458 PetscFunctionReturn(0); 459 } 460 461