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