1 #include <petscsf.h> 2 #include <petscdmswarm.h> 3 #include <petscdmda.h> 4 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 5 #include "../src/dm/impls/swarm/data_bucket.h" 6 #include "../src/dm/impls/swarm/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);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(ierr); 202 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr); 203 204 #if 1 205 { 206 PetscInt *p_cellid; 207 PetscInt npoints_curr,range = 0; 208 PetscSFNode *sf_cells; 209 210 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);CHKERRQ(ierr); 211 ierr = PetscMalloc1(npoints_curr, &sf_cells);CHKERRQ(ierr); 212 213 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 214 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);CHKERRQ(ierr); 215 for (p=0; p<npoints_curr; p++) { 216 217 sf_cells[p].rank = 0; 218 sf_cells[p].index = p_cellid[p]; 219 if (p_cellid[p] > range) { 220 range = p_cellid[p]; 221 } 222 } 223 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);CHKERRQ(ierr); 224 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 225 226 /*ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell);CHKERRQ(ierr);*/ 227 ierr = PetscSFCreate(PETSC_COMM_SELF,&sfcell);CHKERRQ(ierr); 228 ierr = PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);CHKERRQ(ierr); 229 } 230 #endif 231 232 ierr = DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);CHKERRQ(ierr); 233 ierr = DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);CHKERRQ(ierr); 234 ierr = DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);CHKERRQ(ierr); 235 236 if (error_check) { 237 ierr = DMSwarmGetSize(dm,&npointsg);CHKERRQ(ierr); 238 } 239 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 240 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 241 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 242 for (p=0; p<npoints; p++) { 243 rankval[p] = LA_sfcell[p].index; 244 } 245 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 246 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 247 248 if (size > 1) { 249 ierr = DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);CHKERRQ(ierr); 250 } else { 251 DMSwarmDataField PField; 252 PetscInt npoints_curr; 253 254 /* remove points which the domain */ 255 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 256 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 257 258 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);CHKERRQ(ierr); 259 for (p=0; p<npoints_curr; p++) { 260 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 261 /* kill point */ 262 ierr = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 263 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 264 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */ 265 p--; /* check replacement point */ 266 } 267 } 268 ierr = DMSwarmGetSize(dm,&npoints_prior_migration);CHKERRQ(ierr); 269 270 } 271 272 /* locate points newly received */ 273 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); 274 275 #if 0 276 { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */ 277 PetscScalar *LA_coor; 278 PetscInt bs; 279 DMSwarmDataField PField; 280 281 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 282 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);CHKERRQ(ierr); 283 ierr = DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr); 284 285 ierr = VecDestroy(&pos);CHKERRQ(ierr); 286 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 287 288 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 289 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 290 for (p=0; p<npoints2; p++) { 291 rankval[p] = LA_sfcell[p].index; 292 } 293 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 294 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 295 296 /* remove points which left processor */ 297 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 298 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 299 300 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); 301 for (p=0; p<npoints2; p++) { 302 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 303 /* kill point */ 304 ierr = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 305 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 306 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */ 307 p--; /* check replacement point */ 308 } 309 } 310 } 311 #endif 312 313 { /* this performs two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */ 314 PetscScalar *LA_coor; 315 PetscInt npoints_from_neighbours,bs; 316 DMSwarmDataField PField; 317 318 npoints_from_neighbours = npoints2 - npoints_prior_migration; 319 320 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 321 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);CHKERRQ(ierr); 322 323 ierr = DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr); 324 325 ierr = VecDestroy(&pos);CHKERRQ(ierr); 326 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr); 327 328 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 329 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 330 for (p=0; p<npoints_from_neighbours; p++) { 331 rankval[npoints_prior_migration + p] = LA_sfcell[p].index; 332 } 333 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 334 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 335 336 /* remove points which left processor */ 337 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr); 338 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); 339 340 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); 341 for (p=npoints_prior_migration; p<npoints2; p++) { 342 if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) { 343 /* kill point */ 344 ierr = DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 345 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 346 ierr = DMSwarmDataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */ 347 p--; /* check replacement point */ 348 } 349 } 350 } 351 352 { 353 PetscInt *p_cellid; 354 355 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); 356 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 357 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);CHKERRQ(ierr); 358 for (p=0; p<npoints2; p++) { 359 p_cellid[p] = rankval[p]; 360 } 361 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);CHKERRQ(ierr); 362 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 363 } 364 365 /* check for error on removed points */ 366 if (error_check) { 367 ierr = DMSwarmGetSize(dm,&npoints2g);CHKERRQ(ierr); 368 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); 369 } 370 PetscFunctionReturn(0); 371 } 372 373 PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points) 374 { 375 PetscFunctionBegin; 376 PetscFunctionReturn(0); 377 } 378 379 /* 380 Redundant as this assumes points can only be sent to a single rank 381 */ 382 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize) 383 { 384 DM_Swarm *swarm = (DM_Swarm*)dm->data; 385 PetscErrorCode ierr; 386 DMSwarmDataEx de; 387 PetscInt p,npoints,*rankval,n_points_recv; 388 PetscMPIInt rank,nrank,negrank; 389 void *point_buffer,*recv_points; 390 size_t sizeof_dmswarm_point; 391 392 PetscFunctionBegin; 393 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr); 394 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 395 *globalsize = npoints; 396 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 397 ierr = DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr); 398 ierr = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 399 for (p=0; p<npoints; p++) { 400 negrank = rankval[p]; 401 if (negrank < 0) { 402 nrank = -negrank - 1; 403 ierr = DMSwarmDataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 404 } 405 } 406 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 407 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 408 for (p=0; p<npoints; p++) { 409 negrank = rankval[p]; 410 if (negrank < 0) { 411 nrank = -negrank - 1; 412 ierr = DMSwarmDataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 413 } 414 } 415 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 416 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 417 ierr = DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 418 for (p=0; p<npoints; p++) { 419 negrank = rankval[p]; 420 if (negrank < 0) { 421 nrank = -negrank - 1; 422 rankval[p] = nrank; 423 /* copy point into buffer */ 424 ierr = DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 425 /* insert point buffer into DMSwarmDataExchanger */ 426 ierr = DMSwarmDataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 427 rankval[p] = negrank; 428 } 429 } 430 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 431 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 432 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 433 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 434 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 435 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 436 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr); 437 for (p=0; p<n_points_recv; p++) { 438 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point); 439 440 ierr = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 441 } 442 ierr = DMSwarmDataExView(de);CHKERRQ(ierr); 443 ierr = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 444 ierr = DMSwarmDataExDestroy(de);CHKERRQ(ierr); 445 PetscFunctionReturn(0); 446 } 447 448 typedef struct { 449 PetscMPIInt owner_rank; 450 PetscReal min[3],max[3]; 451 } CollectBBox; 452 453 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize) 454 { 455 DM_Swarm * swarm = (DM_Swarm*)dm->data; 456 PetscErrorCode ierr; 457 DMSwarmDataEx de; 458 PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells; 459 PetscMPIInt rank,nrank; 460 void *point_buffer,*recv_points; 461 size_t sizeof_dmswarm_point,sizeof_bbox_ctx; 462 PetscBool isdmda; 463 CollectBBox *bbox,*recv_bbox; 464 const PetscMPIInt *dmneighborranks; 465 DM dmcell; 466 467 PetscFunctionBegin; 468 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr); 469 470 ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 471 if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 472 isdmda = PETSC_FALSE; 473 PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda); 474 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox"); 475 476 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 477 sizeof_bbox_ctx = sizeof(CollectBBox); 478 PetscMalloc1(1,&bbox); 479 bbox->owner_rank = rank; 480 481 /* compute the bounding box based on the overlapping / stenctil size */ 482 { 483 Vec lcoor; 484 485 ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr); 486 if (dim >= 1) { 487 ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr); 488 ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr); 489 } 490 if (dim >= 2) { 491 ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr); 492 ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr); 493 } 494 if (dim == 3) { 495 ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr); 496 ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr); 497 } 498 } 499 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 500 *globalsize = npoints; 501 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 502 ierr = DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr); 503 /* use DMDA neighbours */ 504 ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr); 505 if (dim == 1) { 506 neighbour_cells = 3; 507 } else if (dim == 2) { 508 neighbour_cells = 9; 509 } else { 510 neighbour_cells = 27; 511 } 512 ierr = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 513 for (p=0; p<neighbour_cells; p++) { 514 if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) { 515 ierr = DMSwarmDataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr); 516 } 517 } 518 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 519 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 520 for (p=0; p<neighbour_cells; p++) { 521 if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) { 522 ierr = DMSwarmDataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr); 523 } 524 } 525 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 526 /* send bounding boxes */ 527 ierr = DMSwarmDataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr); 528 for (p=0; p<neighbour_cells; p++) { 529 nrank = dmneighborranks[p]; 530 if ((nrank >= 0) && (nrank != rank)) { 531 /* insert bbox buffer into DMSwarmDataExchanger */ 532 ierr = DMSwarmDataExPackData(de,nrank,1,bbox);CHKERRQ(ierr); 533 } 534 } 535 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 536 /* recv bounding boxes */ 537 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 538 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 539 ierr = DMSwarmDataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr); 540 /* Wrong, should not be using PETSC_COMM_WORLD */ 541 for (p=0; p<n_bbox_recv; p++) { 542 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, 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]);CHKERRQ(ierr); 544 } 545 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);CHKERRQ(ierr); 546 /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */ 547 ierr = DMSwarmDataExInitializeSendCount(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 = DMSwarmDataExAddToSendCount(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 = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 564 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 565 ierr = DMSwarmDataExPackInitialize(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 = DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 576 /* insert point buffer into DMSwarmDataExchanger */ 577 ierr = DMSwarmDataExPackData(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 = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 585 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 586 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 587 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 588 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 589 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 590 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_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 = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 595 } 596 ierr = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 597 PetscFree(bbox); 598 ierr = DMSwarmDataExView(de);CHKERRQ(ierr); 599 ierr = DMSwarmDataExDestroy(de);CHKERRQ(ierr); 600 PetscFunctionReturn(0); 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 DMSwarmDataEx de; 617 PetscInt p,r,npoints,n_points_recv; 618 PetscMPIInt size,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),&size);CHKERRMPI(ierr); 626 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRMPI(ierr); 627 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 628 *globalsize = npoints; 629 /* Broadcast user context */ 630 PetscMalloc(ctx_size*size,&ctxlist); 631 ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr); 632 ierr = PetscMalloc1(size,&n2collect);CHKERRQ(ierr); 633 ierr = PetscMalloc1(size,&collectlist);CHKERRQ(ierr); 634 for (r=0; r<size; 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 = DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr); 649 /* Define topology */ 650 ierr = DMSwarmDataExTopologyInitialize(de);CHKERRQ(ierr); 651 for (r=0; r<size; r++) { 652 if (n2collect[r] > 0) { 653 ierr = DMSwarmDataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr); 654 } 655 } 656 ierr = DMSwarmDataExTopologyFinalize(de);CHKERRQ(ierr); 657 /* Define send counts */ 658 ierr = DMSwarmDataExInitializeSendCount(de);CHKERRQ(ierr); 659 for (r=0; r<size; r++) { 660 if (n2collect[r] > 0) { 661 ierr = DMSwarmDataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr); 662 } 663 } 664 ierr = DMSwarmDataExFinalizeSendCount(de);CHKERRQ(ierr); 665 /* Pack data */ 666 ierr = DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 667 ierr = DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 668 for (r=0; r<size; r++) { 669 for (p=0; p<n2collect[r]; p++) { 670 ierr = DMSwarmDataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr); 671 /* insert point buffer into the data exchanger */ 672 ierr = DMSwarmDataExPackData(de,r,1,point_buffer);CHKERRQ(ierr); 673 } 674 } 675 ierr = DMSwarmDataExPackFinalize(de);CHKERRQ(ierr); 676 /* Scatter */ 677 ierr = DMSwarmDataExBegin(de);CHKERRQ(ierr); 678 ierr = DMSwarmDataExEnd(de);CHKERRQ(ierr); 679 /* Collect data in DMSwarm container */ 680 ierr = DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 681 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 682 ierr = DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_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 = DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 687 } 688 /* Release memory */ 689 for (r=0; r<size; 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 = DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 696 ierr = DMSwarmDataExView(de);CHKERRQ(ierr); 697 ierr = DMSwarmDataExDestroy(de);CHKERRQ(ierr); 698 PetscFunctionReturn(0); 699 } 700