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