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