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