1 #include <petscsf.h> 2 #include <petscdmswarm.h> 3 #include <petscdmda.h> 4 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 5 #include "data_bucket.h" 6 #include "data_ex.h" 7 8 /* 9 User loads desired location (MPI rank) into field DMSwarm_rank 10 */ 11 PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points) 12 { 13 DM_Swarm *swarm = (DM_Swarm*)dm->data; 14 PetscErrorCode ierr; 15 DMSwarmDataEx de; 16 PetscInt p,npoints,*rankval,n_points_recv; 17 PetscMPIInt rank,nrank; 18 void *point_buffer,*recv_points; 19 size_t sizeof_dmswarm_point; 20 21 PetscFunctionBegin; 22 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 23 24 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 25 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 26 ierr = DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0, &de);CHKERRQ(ierr); 27 ierr = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 28 for (p = 0; p < npoints; ++p) { 29 nrank = rankval[p]; 30 if (nrank != rank) { 31 ierr = DMSwarmDataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 32 } 33 } 34 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 35 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 36 for (p=0; p<npoints; p++) { 37 nrank = rankval[p]; 38 if (nrank != rank) { 39 ierr = DMSwarmDataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 40 } 41 } 42 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 43 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 44 ierr = DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 45 for (p=0; p<npoints; p++) { 46 nrank = rankval[p]; 47 if (nrank != rank) { 48 /* copy point into buffer */ 49 ierr = DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 50 /* insert point buffer into DMSwarmDataExchanger */ 51 ierr = DMSwarmDataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 52 } 53 } 54 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 55 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 56 57 if (remove_sent_points) { 58 DMSwarmDataField gfield; 59 60 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);CHKERRQ(ierr); 61 ierr = DMSwarmDataFieldGetAccess(gfield);CHKERRQ(ierr); 62 ierr = DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);CHKERRQ(ierr); 63 64 /* remove points which left processor */ 65 ierr = DMSwarmDataBucketGetSizes(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 = DMSwarmDataFieldRestoreAccess(gfield);CHKERRQ(ierr); 71 72 ierr = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 73 74 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 75 ierr = DMSwarmDataFieldGetAccess(gfield);CHKERRQ(ierr); 76 ierr = DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);CHKERRQ(ierr); 77 p--; /* check replacement point */ 78 } 79 } 80 ierr = DMSwarmDataFieldRestoreEntries(gfield,(void**)&rankval);CHKERRQ(ierr); 81 ierr = DMSwarmDataFieldRestoreAccess(gfield);CHKERRQ(ierr); 82 } 83 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 84 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 85 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 86 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 87 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr); 88 for (p=0; p<n_points_recv; p++) { 89 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 90 91 ierr = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 92 } 93 ierr = DMSwarmDataExView(de);CHKERRQ(ierr); 94 ierr = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 95 ierr = DMSwarmDataExDestroy(de);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 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 DMSwarmDataEx 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 PetscFunctionBegin; 113 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 114 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 115 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 116 ierr = DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr); 117 ierr = DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);CHKERRQ(ierr); 118 ierr = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 119 for (r=0; r<nneighbors; r++) { 120 _rank = neighbourranks[r]; 121 if ((_rank != rank) && (_rank > 0)) { 122 ierr = DMSwarmDataExTopologyAddNeighbour(de,_rank);CHKERRQ(ierr); 123 } 124 } 125 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 126 ierr = DMSwarmDataExTopologyGetNeighbours(de,&mynneigh,&myneigh);CHKERRQ(ierr); 127 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 128 for (p=0; p<npoints; p++) { 129 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 130 for (r=0; r<mynneigh; r++) { 131 _rank = myneigh[r]; 132 ierr = DMSwarmDataExAddToSendCount(de,_rank,1);CHKERRQ(ierr); 133 } 134 } 135 } 136 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 137 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 138 ierr = DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 139 for (p=0; p<npoints; p++) { 140 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 141 for (r=0; r<mynneigh; r++) { 142 _rank = myneigh[r]; 143 /* copy point into buffer */ 144 ierr = DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 145 /* insert point buffer into DMSwarmDataExchanger */ 146 ierr = DMSwarmDataExPackData(de,_rank,1,point_buffer);CHKERRQ(ierr); 147 } 148 } 149 } 150 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 151 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 152 if (remove_sent_points) { 153 DMSwarmDataField PField; 154 155 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 156 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 157 /* remove points which left processor */ 158 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 159 for (p=0; p<npoints; p++) { 160 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 161 /* kill point */ 162 ierr = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 163 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 164 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */ 165 p--; /* check replacement point */ 166 } 167 } 168 } 169 ierr = DMSwarmDataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);CHKERRQ(ierr); 170 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 171 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 172 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 173 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 174 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr); 175 for (p=0; p<n_points_recv; p++) { 176 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 177 178 ierr = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 179 } 180 ierr = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 181 ierr = DMSwarmDataExDestroy(de);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points) 186 { 187 DM_Swarm *swarm = (DM_Swarm*)dm->data; 188 PetscErrorCode ierr; 189 PetscInt p,npoints,npointsg=0,npoints2,npoints2g,*rankval,npoints_prior_migration; 190 PetscSF sfcell = NULL; 191 const PetscSFNode *LA_sfcell; 192 DM dmcell; 193 Vec pos; 194 PetscBool error_check = swarm->migrate_error_on_missing_point; 195 PetscMPIInt size,rank; 196 197 PetscFunctionBegin; 198 ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 199 if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 200 201 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr); 202 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 203 204 #if 1 205 { 206 PetscInt *p_cellid; 207 PetscInt npoints_curr,range = 0; 208 PetscSFNode *sf_cells; 209 210 211 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);CHKERRQ(ierr); 212 ierr = PetscMalloc1(npoints_curr, &sf_cells);CHKERRQ(ierr); 213 214 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 215 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);CHKERRQ(ierr); 216 for (p=0; p<npoints_curr; p++) { 217 218 sf_cells[p].rank = 0; 219 sf_cells[p].index = p_cellid[p]; 220 if (p_cellid[p] > range) { 221 range = p_cellid[p]; 222 } 223 } 224 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);CHKERRQ(ierr); 225 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 226 227 /*ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell);CHKERRQ(ierr);*/ 228 ierr = PetscSFCreate(PETSC_COMM_SELF,&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 = DMSwarmDataBucketGetSizes(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 (size > 1) { 250 ierr = DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);CHKERRQ(ierr); 251 } else { 252 DMSwarmDataField PField; 253 PetscInt npoints_curr; 254 255 /* remove points which the domain */ 256 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 257 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 258 259 ierr = DMSwarmDataBucketGetSizes(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 = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 264 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 265 ierr = DMSwarmDataFieldGetEntries(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 = DMSwarmDataBucketGetSizes(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 DMSwarmDataField 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 = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 299 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 300 301 ierr = DMSwarmDataBucketGetSizes(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 = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 306 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 307 ierr = DMSwarmDataFieldGetEntries(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 DMSwarmDataField 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 = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 339 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 340 341 ierr = DMSwarmDataBucketGetSizes(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 = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 346 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 347 ierr = DMSwarmDataFieldGetEntries(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 = DMSwarmDataBucketGetSizes(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 DMSwarmDataEx 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 = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 396 *globalsize = npoints; 397 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 398 ierr = DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr); 399 ierr = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 400 for (p=0; p<npoints; p++) { 401 negrank = rankval[p]; 402 if (negrank < 0) { 403 nrank = -negrank - 1; 404 ierr = DMSwarmDataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 405 } 406 } 407 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 408 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 409 for (p=0; p<npoints; p++) { 410 negrank = rankval[p]; 411 if (negrank < 0) { 412 nrank = -negrank - 1; 413 ierr = DMSwarmDataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 414 } 415 } 416 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 417 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 418 ierr = DMSwarmDataExPackInitialize(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 = DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 426 /* insert point buffer into DMSwarmDataExchanger */ 427 ierr = DMSwarmDataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 428 rankval[p] = negrank; 429 } 430 } 431 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 432 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 433 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 434 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 435 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 436 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 437 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_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 = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 442 } 443 ierr = DMSwarmDataExView(de);CHKERRQ(ierr); 444 ierr = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 445 ierr = DMSwarmDataExDestroy(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 DMSwarmDataEx 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 = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 501 *globalsize = npoints; 502 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 503 ierr = DMSwarmDataExCreate(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 = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 514 for (p=0; p<neighbour_cells; p++) { 515 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 516 ierr = DMSwarmDataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr); 517 } 518 } 519 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 520 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 521 for (p=0; p<neighbour_cells; p++) { 522 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 523 ierr = DMSwarmDataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr); 524 } 525 } 526 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 527 /* send bounding boxes */ 528 ierr = DMSwarmDataExPackInitialize(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 DMSwarmDataExchanger */ 533 ierr = DMSwarmDataExPackData(de,nrank,1,bbox);CHKERRQ(ierr); 534 } 535 } 536 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 537 /* recv bounding boxes */ 538 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 539 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 540 ierr = DMSwarmDataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr); 541 /* Wrong, should not be using PETSC_COMM_WORLD */ 542 for (p=0; p<n_bbox_recv; p++) { 543 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[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]);CHKERRQ(ierr); 545 } 546 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);CHKERRQ(ierr); 547 /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */ 548 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 549 for (pk=0; pk<n_bbox_recv; pk++) { 550 PetscReal *array_x,*array_y; 551 552 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 553 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 554 for (p=0; p<npoints; p++) { 555 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 556 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 557 ierr = DMSwarmDataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr); 558 } 559 } 560 } 561 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 562 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 563 } 564 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 565 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 566 ierr = DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 567 for (pk=0; pk<n_bbox_recv; pk++) { 568 PetscReal *array_x,*array_y; 569 570 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 571 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 572 for (p=0; p<npoints; p++) { 573 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 574 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 575 /* copy point into buffer */ 576 ierr = DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 577 /* insert point buffer into DMSwarmDataExchanger */ 578 ierr = DMSwarmDataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr); 579 } 580 } 581 } 582 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 583 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 584 } 585 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 586 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 587 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 588 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 589 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 590 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 591 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr); 592 for (p=0; p<n_points_recv; p++) { 593 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 594 595 ierr = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 596 } 597 ierr = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 598 PetscFree(bbox); 599 ierr = DMSwarmDataExView(de);CHKERRQ(ierr); 600 ierr = DMSwarmDataExDestroy(de);CHKERRQ(ierr); 601 PetscFunctionReturn(0); 602 } 603 604 605 /* General collection when no order, or neighbour information is provided */ 606 /* 607 User provides context and collect() method 608 Broadcast user context 609 610 for each context / rank { 611 collect(swarm,context,n,list) 612 } 613 */ 614 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize) 615 { 616 DM_Swarm *swarm = (DM_Swarm*)dm->data; 617 PetscErrorCode ierr; 618 DMSwarmDataEx de; 619 PetscInt p,r,npoints,n_points_recv; 620 PetscMPIInt size,rank; 621 void *point_buffer,*recv_points; 622 void *ctxlist; 623 PetscInt *n2collect,**collectlist; 624 size_t sizeof_dmswarm_point; 625 626 PetscFunctionBegin; 627 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr); 628 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 629 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 630 *globalsize = npoints; 631 /* Broadcast user context */ 632 PetscMalloc(ctx_size*size,&ctxlist); 633 ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 634 ierr = PetscMalloc1(size,&n2collect);CHKERRQ(ierr); 635 ierr = PetscMalloc1(size,&collectlist);CHKERRQ(ierr); 636 for (r=0; r<size; r++) { 637 PetscInt _n2collect; 638 PetscInt *_collectlist; 639 void *_ctx_r; 640 641 _n2collect = 0; 642 _collectlist = NULL; 643 if (r != rank) { /* don't collect data from yourself */ 644 _ctx_r = (void*)( (char*)ctxlist + r * ctx_size ); 645 ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr); 646 } 647 n2collect[r] = _n2collect; 648 collectlist[r] = _collectlist; 649 } 650 ierr = DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr); 651 /* Define topology */ 652 ierr = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 653 for (r=0; r<size; r++) { 654 if (n2collect[r] > 0) { 655 ierr = DMSwarmDataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr); 656 } 657 } 658 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 659 /* Define send counts */ 660 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 661 for (r=0; r<size; r++) { 662 if (n2collect[r] > 0) { 663 ierr = DMSwarmDataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr); 664 } 665 } 666 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 667 /* Pack data */ 668 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 669 ierr = DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 670 for (r=0; r<size; r++) { 671 for (p=0; p<n2collect[r]; p++) { 672 ierr = DMSwarmDataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr); 673 /* insert point buffer into the data exchanger */ 674 ierr = DMSwarmDataExPackData(de,r,1,point_buffer);CHKERRQ(ierr); 675 } 676 } 677 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 678 /* Scatter */ 679 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 680 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 681 /* Collect data in DMSwarm container */ 682 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 683 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 684 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr); 685 for (p=0; p<n_points_recv; p++) { 686 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 687 688 ierr = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 689 } 690 /* Release memory */ 691 for (r=0; r<size; r++) { 692 if (collectlist[r]) PetscFree(collectlist[r]); 693 } 694 ierr = PetscFree(collectlist);CHKERRQ(ierr); 695 ierr = PetscFree(n2collect);CHKERRQ(ierr); 696 ierr = PetscFree(ctxlist);CHKERRQ(ierr); 697 ierr = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 698 ierr = DMSwarmDataExView(de);CHKERRQ(ierr); 699 ierr = DMSwarmDataExDestroy(de);CHKERRQ(ierr); 700 PetscFunctionReturn(0); 701 } 702 703