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,len,*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 if (commsize == 1) PetscFunctionReturn(0); 215 216 ierr = DMSwarmCreateGlobalVectorFromField(dm,DMSwarmPICField_coor,&pos);CHKERRQ(ierr); 217 ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr); 218 ierr = DMSwarmDestroyGlobalVectorFromField(dm,DMSwarmPICField_coor,&pos);CHKERRQ(ierr); 219 220 if (error_check) { 221 ierr = DMSwarmGetSize(dm,&npointsg);CHKERRQ(ierr); 222 } 223 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 224 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 225 ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr); 226 for (p=0; p<npoints; p++) { 227 rankval[p] = LA_iscell[p]; 228 } 229 ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr); 230 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 231 ierr = ISDestroy(&iscell);CHKERRQ(ierr); 232 233 ierr = DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);CHKERRQ(ierr); 234 235 /* locate points newly recevied */ 236 ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); 237 #if 0 238 len = npoints2 - npoints_prior_migration; 239 if (len > 0) { 240 PetscScalar *LA_coor; 241 PetscInt bs; 242 243 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 244 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*len,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);CHKERRQ(ierr); 245 246 ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr); 247 248 ierr = VecDestroy(&pos);CHKERRQ(ierr); 249 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 250 251 ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr); 252 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 253 for (p=0; p<len; p++) { 254 rankval[npoints_prior_migration+p] = LA_iscell[p]; 255 } 256 ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr); 257 ierr = ISDestroy(&iscell);CHKERRQ(ierr); 258 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 259 260 /* remove points which left processor */ 261 ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 262 ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 263 264 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 265 for (p=npoints_prior_migration; p<npoints2; p++) { 266 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 267 /* kill point */ 268 ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 269 ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 270 ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */ 271 p--; /* check replacement point */ 272 } 273 } 274 275 } 276 #endif 277 278 { 279 PetscScalar *LA_coor; 280 PetscInt bs; 281 DataField PField; 282 283 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 284 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);CHKERRQ(ierr); 285 ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr); 286 287 ierr = VecDestroy(&pos);CHKERRQ(ierr); 288 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 289 290 ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr); 291 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 292 for (p=0; p<npoints2; p++) { 293 rankval[p] = LA_iscell[p]; 294 } 295 ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr); 296 ierr = ISDestroy(&iscell);CHKERRQ(ierr); 297 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 298 299 /* remove points which left processor */ 300 ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 301 ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 302 303 ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); 304 for (p=0; p<npoints2; p++) { 305 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 306 /* kill point */ 307 ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 308 ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 309 ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */ 310 p--; /* check replacement point */ 311 } 312 } 313 314 } 315 316 /* check for error on removed points */ 317 if (error_check) { 318 ierr = DMSwarmGetSize(dm,&npoints2g);CHKERRQ(ierr); 319 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); 320 } 321 PetscFunctionReturn(0); 322 } 323 324 #undef __FUNCT__ 325 #define __FUNCT__ "DMSwarmMigrate_CellDMExact" 326 PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points) 327 { 328 // DM_Swarm *swarm = (DM_Swarm*)dm->data; 329 // PetscErrorCode ierr; 330 331 PetscFunctionReturn(0); 332 } 333 334 /* 335 Redundant as this assumes points can only be sent to a single rank 336 */ 337 #undef __FUNCT__ 338 #define __FUNCT__ "DMSwarmMigrate_GlobalToLocal_Basic" 339 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize) 340 { 341 DM_Swarm *swarm = (DM_Swarm*)dm->data; 342 PetscErrorCode ierr; 343 DataEx de; 344 PetscInt p,npoints,*rankval,n_points_recv; 345 PetscMPIInt rank,nrank,negrank; 346 void *point_buffer,*recv_points; 347 size_t sizeof_dmswarm_point; 348 349 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 350 351 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 352 *globalsize = npoints; 353 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 354 355 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 356 357 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 358 for (p=0; p<npoints; p++) { 359 negrank = rankval[p]; 360 if (negrank < 0) { 361 nrank = -negrank - 1; 362 ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 363 } 364 } 365 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 366 367 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 368 for (p=0; p<npoints; p++) { 369 negrank = rankval[p]; 370 if (negrank < 0) { 371 nrank = -negrank - 1; 372 ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 373 } 374 } 375 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 376 377 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 378 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 379 for (p=0; p<npoints; p++) { 380 negrank = rankval[p]; 381 if (negrank < 0) { 382 nrank = -negrank - 1; 383 rankval[p] = nrank; 384 /* copy point into buffer */ 385 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 386 /* insert point buffer into DataExchanger */ 387 ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 388 rankval[p] = negrank; 389 } 390 } 391 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 392 393 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 394 395 ierr = DataExBegin(de);CHKERRQ(ierr); 396 ierr = DataExEnd(de);CHKERRQ(ierr); 397 398 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 399 400 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 401 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 402 for (p=0; p<n_points_recv; p++) { 403 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 404 405 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 406 } 407 ierr = DataExView(de);CHKERRQ(ierr); 408 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 409 ierr = DataExDestroy(de);CHKERRQ(ierr); 410 411 PetscFunctionReturn(0); 412 } 413 414 typedef struct { 415 PetscMPIInt owner_rank; 416 PetscReal min[3],max[3]; 417 } CollectBBox; 418 419 #undef __FUNCT__ 420 #define __FUNCT__ "DMSwarmCollect_DMDABoundingBox" 421 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize) 422 { 423 DM_Swarm *swarm = (DM_Swarm*)dm->data; 424 PetscErrorCode ierr; 425 DataEx de; 426 PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells; 427 PetscMPIInt rank,nrank; 428 void *point_buffer,*recv_points; 429 size_t sizeof_dmswarm_point,sizeof_bbox_ctx; 430 PetscBool isdmda; 431 CollectBBox *bbox,*recv_bbox; 432 const PetscMPIInt *dmneighborranks; 433 DM dmcell; 434 435 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 436 437 ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 438 if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 439 440 isdmda = PETSC_FALSE; 441 PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda); 442 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox"); 443 444 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 445 446 sizeof_bbox_ctx = sizeof(CollectBBox); 447 PetscMalloc1(1,&bbox); 448 bbox->owner_rank = rank; 449 450 /* compute the bounding box based on the overlapping / stenctil size */ 451 //ierr = DMDAGetLocalBoundingBox(dmcell,bbox->min,bbox->max);CHKERRQ(ierr); 452 { 453 Vec lcoor; 454 455 ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr); 456 457 if (dim >= 1) { 458 ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr); 459 ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr); 460 } 461 462 if (dim >= 2) { 463 ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr); 464 ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr); 465 } 466 467 if (dim == 3) { 468 ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr); 469 ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr); 470 } 471 } 472 473 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 474 *globalsize = npoints; 475 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 476 477 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 478 479 /* use DMDA neighbours */ 480 ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr); 481 if (dim == 1) { 482 neighbour_cells = 3; 483 } else if (dim == 2) { 484 neighbour_cells = 9; 485 } else { 486 neighbour_cells = 27; 487 } 488 489 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 490 for (p=0; p<neighbour_cells; p++) { 491 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 492 ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr); 493 } 494 } 495 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 496 497 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 498 for (p=0; p<neighbour_cells; p++) { 499 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 500 ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr); 501 } 502 } 503 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 504 505 506 /* send bounding boxes */ 507 ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr); 508 for (p=0; p<neighbour_cells; p++) { 509 nrank = dmneighborranks[p]; 510 if ( (nrank >= 0) && (nrank != rank) ) { 511 /* insert bbox buffer into DataExchanger */ 512 ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr); 513 } 514 } 515 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 516 517 /* recv bounding boxes */ 518 ierr = DataExBegin(de);CHKERRQ(ierr); 519 ierr = DataExEnd(de);CHKERRQ(ierr); 520 521 ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr); 522 523 524 for (p=0; p<n_bbox_recv; p++) { 525 printf("[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank, 526 recv_bbox[p].min[0],recv_bbox[p].max[0],recv_bbox[p].min[1],recv_bbox[p].max[1]); 527 } 528 529 /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */ 530 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 531 for (pk=0; pk<n_bbox_recv; pk++) { 532 PetscReal *array_x,*array_y; 533 534 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 535 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 536 537 for (p=0; p<npoints; p++) { 538 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 539 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 540 541 ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr); 542 543 } 544 } 545 } 546 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 547 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 548 } 549 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 550 551 552 553 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 554 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 555 for (pk=0; pk<n_bbox_recv; pk++) { 556 PetscReal *array_x,*array_y; 557 558 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 559 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 560 561 for (p=0; p<npoints; p++) { 562 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 563 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 564 // copy point into buffer // 565 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 566 // insert point buffer into DataExchanger // 567 ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr); 568 } 569 } 570 } 571 572 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 573 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 574 } 575 576 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 577 578 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 579 580 ierr = DataExBegin(de);CHKERRQ(ierr); 581 ierr = DataExEnd(de);CHKERRQ(ierr); 582 583 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 584 585 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 586 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 587 for (p=0; p<n_points_recv; p++) { 588 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 589 590 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 591 } 592 593 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 594 PetscFree(bbox); 595 ierr = DataExView(de);CHKERRQ(ierr); 596 ierr = DataExDestroy(de);CHKERRQ(ierr); 597 598 PetscFunctionReturn(0); 599 } 600 601 602 /* General collection when no order, or neighbour information is provided */ 603 /* 604 User provides context and collect() method 605 Broadcast user context 606 607 for each context / rank { 608 collect(swarm,context,n,list) 609 } 610 611 */ 612 #undef __FUNCT__ 613 #define __FUNCT__ "DMSwarmCollect_General" 614 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize) 615 { 616 DM_Swarm *swarm = (DM_Swarm*)dm->data; 617 PetscErrorCode ierr; 618 DataEx de; 619 PetscInt p,r,npoints,n_points_recv; 620 PetscMPIInt commsize,rank; 621 void *point_buffer,*recv_points; 622 void *ctxlist; 623 PetscInt *n2collect,**collectlist; 624 size_t sizeof_dmswarm_point; 625 626 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr); 627 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 628 629 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 630 *globalsize = npoints; 631 632 /* Broadcast user context */ 633 PetscMalloc(ctx_size*commsize,&ctxlist); 634 ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 635 636 PetscMalloc1(commsize,&n2collect); 637 PetscMalloc1(commsize,&collectlist); 638 639 for (r=0; r<commsize; r++) { 640 PetscInt _n2collect; 641 PetscInt *_collectlist; 642 void *_ctx_r; 643 644 _n2collect = 0; 645 _collectlist = NULL; 646 647 if (r != rank) { /* don't collect data from yourself */ 648 _ctx_r = (void*)( (char*)ctxlist + r * ctx_size ); 649 ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr); 650 } 651 652 n2collect[r] = _n2collect; 653 collectlist[r] = _collectlist; 654 } 655 656 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 657 658 /* Define topology */ 659 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 660 for (r=0; r<commsize; r++) { 661 if (n2collect[r] > 0) { 662 ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr); 663 } 664 } 665 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 666 667 /* Define send counts */ 668 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 669 for (r=0; r<commsize; r++) { 670 if (n2collect[r] > 0) { 671 ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr); 672 } 673 } 674 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 675 676 /* Pack data */ 677 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 678 679 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 680 for (r=0; r<commsize; r++) { 681 682 for (p=0; p<n2collect[r]; p++) { 683 ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr); 684 /* insert point buffer into the data exchanger */ 685 ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr); 686 } 687 } 688 689 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 690 691 /* Scatter */ 692 ierr = DataExBegin(de);CHKERRQ(ierr); 693 ierr = DataExEnd(de);CHKERRQ(ierr); 694 695 /* Collect data in DMSwarm container */ 696 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 697 698 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 699 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 700 for (p=0; p<n_points_recv; p++) { 701 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 702 703 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 704 } 705 706 /* Release memory */ 707 for (r=0; r<commsize; r++) { 708 if (collectlist[r]) PetscFree(collectlist[r]); 709 } 710 PetscFree(collectlist); 711 PetscFree(n2collect); 712 PetscFree(ctxlist); 713 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 714 ierr = DataExView(de);CHKERRQ(ierr); 715 ierr = DataExDestroy(de);CHKERRQ(ierr); 716 717 PetscFunctionReturn(0); 718 } 719 720