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