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__ "DMSwarmCollect_DMDABoundingBox" 176 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(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 = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 193 if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 194 195 isdmda = PETSC_FALSE; 196 PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda); 197 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox"); 198 199 ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 200 if (dim == 1) { 201 neighbour_cells = 3; 202 } else if (dim == 2) { 203 neighbour_cells = 9; 204 } else { 205 neighbour_cells = 27; 206 } 207 208 sizeof_bbox_ctx = sizeof(CollectBBox); 209 PetscMalloc1(1,&bbox); 210 bbox->owner_rank = rank; 211 212 /* compute the bounding box based on the overlapping / stenctil size */ 213 //ierr = DMDAGetLocalBoundingBox(dmcell,bbox->min,bbox->max);CHKERRQ(ierr); 214 { 215 Vec lcoor; 216 217 ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr); 218 219 if (dim >= 1) { 220 ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr); 221 ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr); 222 } 223 224 if (dim >= 2) { 225 ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr); 226 ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr); 227 } 228 229 if (dim == 3) { 230 ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr); 231 ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr); 232 } 233 } 234 235 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 236 *globalsize = npoints; 237 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 238 239 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 240 241 /* use DMDA neighbours */ 242 ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr); 243 244 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 245 for (p=0; p<neighbour_cells; p++) { 246 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 247 ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr); 248 } 249 } 250 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 251 252 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 253 for (p=0; p<neighbour_cells; p++) { 254 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 255 ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr); 256 } 257 } 258 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 259 260 261 /* send bounding boxes */ 262 ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr); 263 for (p=0; p<neighbour_cells; p++) { 264 nrank = dmneighborranks[p]; 265 if ( (nrank >= 0) && (nrank != rank) ) { 266 /* insert bbox buffer into DataExchanger */ 267 ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr); 268 } 269 } 270 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 271 272 /* recv bounding boxes */ 273 ierr = DataExBegin(de);CHKERRQ(ierr); 274 ierr = DataExEnd(de);CHKERRQ(ierr); 275 276 ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr); 277 278 279 for (p=0; p<n_bbox_recv; p++) { 280 printf("[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank, 281 recv_bbox[p].min[0],recv_bbox[p].max[0],recv_bbox[p].min[1],recv_bbox[p].max[1]); 282 } 283 284 /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */ 285 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 286 for (pk=0; pk<n_bbox_recv; pk++) { 287 PetscReal *array_x,*array_y; 288 289 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 290 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 291 292 for (p=0; p<npoints; p++) { 293 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 294 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 295 296 ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr); 297 298 } 299 } 300 } 301 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 302 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 303 } 304 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 305 306 307 308 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 309 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 310 for (pk=0; pk<n_bbox_recv; pk++) { 311 PetscReal *array_x,*array_y; 312 313 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 314 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 315 316 for (p=0; p<npoints; p++) { 317 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 318 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 319 // copy point into buffer // 320 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 321 // insert point buffer into DataExchanger // 322 ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr); 323 } 324 } 325 } 326 327 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 328 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 329 } 330 331 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 332 333 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 334 335 ierr = DataExBegin(de);CHKERRQ(ierr); 336 ierr = DataExEnd(de);CHKERRQ(ierr); 337 338 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 339 340 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 341 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 342 for (p=0; p<n_points_recv; p++) { 343 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 344 345 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 346 } 347 348 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 349 PetscFree(bbox); 350 ierr = DataExView(de);CHKERRQ(ierr); 351 ierr = DataExDestroy(de);CHKERRQ(ierr); 352 353 PetscFunctionReturn(0); 354 } 355 356 357 /* General collection when no order, or neighbour information is provided */ 358 /* 359 User provides context and collect() method 360 Broadcast user context 361 362 for each context / rank { 363 collect(swarm,context,n,list) 364 } 365 366 */ 367 #undef __FUNCT__ 368 #define __FUNCT__ "DMSwarmCollect_General" 369 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize) 370 { 371 DM_Swarm *swarm = (DM_Swarm*)dm->data; 372 PetscErrorCode ierr; 373 DataEx de; 374 PetscInt p,r,npoints,n_points_recv; 375 PetscMPIInt commsize,rank; 376 void *point_buffer,*recv_points; 377 void *ctxlist; 378 PetscInt *n2collect,**collectlist; 379 size_t sizeof_dmswarm_point; 380 381 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr); 382 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 383 384 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 385 *globalsize = npoints; 386 387 /* Broadcast user context */ 388 PetscMalloc(ctx_size*commsize,&ctxlist); 389 ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 390 391 PetscMalloc1(commsize,&n2collect); 392 PetscMalloc1(commsize,&collectlist); 393 394 for (r=0; r<commsize; r++) { 395 PetscInt _n2collect; 396 PetscInt *_collectlist; 397 void *_ctx_r; 398 399 _n2collect = 0; 400 _collectlist = NULL; 401 402 if (r != rank) { /* don't collect data from yourself */ 403 _ctx_r = (void*)( (char*)ctxlist + r * ctx_size ); 404 ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr); 405 } 406 407 n2collect[r] = _n2collect; 408 collectlist[r] = _collectlist; 409 } 410 411 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 412 413 /* Define topology */ 414 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 415 for (r=0; r<commsize; r++) { 416 if (n2collect[r] > 0) { 417 ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr); 418 } 419 } 420 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 421 422 /* Define send counts */ 423 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 424 for (r=0; r<commsize; r++) { 425 if (n2collect[r] > 0) { 426 ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr); 427 } 428 } 429 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 430 431 /* Pack data */ 432 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 433 434 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 435 for (r=0; r<commsize; r++) { 436 437 for (p=0; p<n2collect[r]; p++) { 438 ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr); 439 /* insert point buffer into the data exchanger */ 440 ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr); 441 } 442 } 443 444 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 445 446 /* Scatter */ 447 ierr = DataExBegin(de);CHKERRQ(ierr); 448 ierr = DataExEnd(de);CHKERRQ(ierr); 449 450 /* Collect data in DMSwarm container */ 451 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 452 453 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 454 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 455 for (p=0; p<n_points_recv; p++) { 456 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 457 458 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 459 } 460 461 /* Release memory */ 462 for (r=0; r<commsize; r++) { 463 if (collectlist[r]) PetscFree(collectlist[r]); 464 } 465 PetscFree(collectlist); 466 PetscFree(n2collect); 467 PetscFree(ctxlist); 468 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 469 ierr = DataExView(de);CHKERRQ(ierr); 470 ierr = DataExDestroy(de);CHKERRQ(ierr); 471 472 PetscFunctionReturn(0); 473 } 474 475