Actual source code: swarm_migrate.c

      1: #include <petscsf.h>
      2: #include <petscdmswarm.h>
      3: #include <petscdmda.h>
      4: #include <petsc/private/dmswarmimpl.h>
      5: #include "../src/dm/impls/swarm/data_bucket.h"
      6: #include "../src/dm/impls/swarm/data_ex.h"

      8: /*
      9:  User loads desired location (MPI rank) into field DMSwarm_rank

     11:  It should be storing the rank information as MPIInt not Int
     12: */
     13: PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm, PetscBool remove_sent_points)
     14: {
     15:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
     16:   DMSwarmDataEx 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:   PetscBool     debug = PETSC_FALSE;

     23:   PetscFunctionBegin;
     24:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

     26:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
     27:   PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
     28:   PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
     29:   PetscCall(DMSwarmDataExTopologyInitialize(de));
     30:   for (p = 0; p < npoints; ++p) {
     31:     nrank = (PetscMPIInt)rankval[p];
     32:     if (nrank != rank) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
     33:   }
     34:   PetscCall(DMSwarmDataExTopologyFinalize(de));
     35:   PetscCall(DMSwarmDataExInitializeSendCount(de));
     36:   for (p = 0; p < npoints; p++) {
     37:     nrank = (PetscMPIInt)rankval[p];
     38:     if (nrank != rank) PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
     39:   }
     40:   PetscCall(DMSwarmDataExFinalizeSendCount(de));
     41:   PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
     42:   PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
     43:   for (p = 0; p < npoints; p++) {
     44:     nrank = (PetscMPIInt)rankval[p];
     45:     if (nrank != rank) {
     46:       /* copy point into buffer */
     47:       PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
     48:       /* insert point buffer into DMSwarmDataExchanger */
     49:       PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
     50:     }
     51:   }
     52:   PetscCall(DMSwarmDataExPackFinalize(de));
     53:   PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));

     55:   if (remove_sent_points) {
     56:     DMSwarmDataField gfield;

     58:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmField_rank, &gfield));
     59:     PetscCall(DMSwarmDataFieldGetAccess(gfield));
     60:     PetscCall(DMSwarmDataFieldGetEntries(gfield, (void **)&rankval));

     62:     /* remove points which left processor */
     63:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
     64:     for (p = 0; p < npoints; p++) {
     65:       nrank = (PetscMPIInt)rankval[p];
     66:       if (nrank != rank) {
     67:         /* kill point */
     68:         PetscCall(DMSwarmDataFieldRestoreAccess(gfield));

     70:         PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));

     72:         PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL)); /* you need to update npoints as the list size decreases! */
     73:         PetscCall(DMSwarmDataFieldGetAccess(gfield));
     74:         PetscCall(DMSwarmDataFieldGetEntries(gfield, (void **)&rankval));
     75:         p--; /* check replacement point */
     76:       }
     77:     }
     78:     PetscCall(DMSwarmDataFieldRestoreEntries(gfield, (void **)&rankval));
     79:     PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
     80:   }
     81:   PetscCall(DMSwarmDataExBegin(de));
     82:   PetscCall(DMSwarmDataExEnd(de));
     83:   PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
     84:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
     85:   PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
     86:   for (p = 0; p < n_points_recv; p++) {
     87:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

     89:     PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
     90:   }
     91:   if (debug) PetscCall(DMSwarmDataExView(de));
     92:   PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
     93:   PetscCall(DMSwarmDataExDestroy(de));
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: static PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm, DM dmcell, PetscBool remove_sent_points, PetscInt *npoints_prior_migration)
     98: {
     99:   DM_Swarm          *swarm = (DM_Swarm *)dm->data;
    100:   DMSwarmDataEx      de;
    101:   PetscInt           r, p, npoints, *p_cellid, n_points_recv;
    102:   PetscMPIInt        rank, _rank;
    103:   const PetscMPIInt *neighbourranks;
    104:   void              *point_buffer, *recv_points;
    105:   size_t             sizeof_dmswarm_point;
    106:   PetscInt           nneighbors;
    107:   PetscMPIInt        mynneigh, *myneigh;

    109:   PetscFunctionBegin;
    110:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    111:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    112:   PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
    113:   PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
    114:   PetscCall(DMGetNeighbors(dmcell, &nneighbors, &neighbourranks));
    115:   PetscCall(DMSwarmDataExTopologyInitialize(de));
    116:   for (r = 0; r < nneighbors; r++) {
    117:     _rank = neighbourranks[r];
    118:     if ((_rank != rank) && (_rank > 0)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, _rank));
    119:   }
    120:   PetscCall(DMSwarmDataExTopologyFinalize(de));
    121:   PetscCall(DMSwarmDataExTopologyGetNeighbours(de, &mynneigh, &myneigh));
    122:   PetscCall(DMSwarmDataExInitializeSendCount(de));
    123:   for (p = 0; p < npoints; p++) {
    124:     if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
    125:       for (r = 0; r < mynneigh; r++) {
    126:         _rank = myneigh[r];
    127:         PetscCall(DMSwarmDataExAddToSendCount(de, _rank, 1));
    128:       }
    129:     }
    130:   }
    131:   PetscCall(DMSwarmDataExFinalizeSendCount(de));
    132:   PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
    133:   PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
    134:   for (p = 0; p < npoints; p++) {
    135:     if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
    136:       for (r = 0; r < mynneigh; r++) {
    137:         _rank = myneigh[r];
    138:         /* copy point into buffer */
    139:         PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
    140:         /* insert point buffer into DMSwarmDataExchanger */
    141:         PetscCall(DMSwarmDataExPackData(de, _rank, 1, point_buffer));
    142:       }
    143:     }
    144:   }
    145:   PetscCall(DMSwarmDataExPackFinalize(de));
    146:   PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
    147:   if (remove_sent_points) {
    148:     DMSwarmDataField PField;

    150:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
    151:     PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));
    152:     /* remove points which left processor */
    153:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    154:     for (p = 0; p < npoints; p++) {
    155:       if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
    156:         /* kill point */
    157:         PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
    158:         PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL)); /* you need to update npoints as the list size decreases! */
    159:         PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));     /* update date point increase realloc performed */
    160:         p--;                                                                   /* check replacement point */
    161:       }
    162:     }
    163:   }
    164:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, npoints_prior_migration, NULL, NULL));
    165:   PetscCall(DMSwarmDataExBegin(de));
    166:   PetscCall(DMSwarmDataExEnd(de));
    167:   PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
    168:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    169:   PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
    170:   for (p = 0; p < n_points_recv; p++) {
    171:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

    173:     PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
    174:   }
    175:   PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
    176:   PetscCall(DMSwarmDataExDestroy(de));
    177:   PetscFunctionReturn(PETSC_SUCCESS);
    178: }

    180: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm, PetscBool remove_sent_points)
    181: {
    182:   DM_Swarm          *swarm = (DM_Swarm *)dm->data;
    183:   PetscInt           p, npoints, npointsg = 0, npoints2, npoints2g, *rankval, *p_cellid, npoints_prior_migration;
    184:   PetscSF            sfcell = NULL;
    185:   const PetscSFNode *LA_sfcell;
    186:   DM                 dmcell;
    187:   Vec                pos;
    188:   PetscBool          error_check = swarm->migrate_error_on_missing_point;
    189:   PetscMPIInt        size, rank;

    191:   PetscFunctionBegin;
    192:   PetscCall(DMSwarmGetCellDM(dm, &dmcell));
    193:   PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");

    195:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    196:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

    198: #if 1
    199:   {
    200:     PetscInt     npoints_curr, range = 0;
    201:     PetscSFNode *sf_cells;

    203:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
    204:     PetscCall(PetscMalloc1(npoints_curr, &sf_cells));

    206:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
    207:     for (p = 0; p < npoints_curr; p++) {
    208:       sf_cells[p].rank  = 0;
    209:       sf_cells[p].index = p_cellid[p];
    210:       if (p_cellid[p] > range) range = p_cellid[p];
    211:     }
    212:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));

    214:     /* PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell)); */
    215:     PetscCall(PetscSFCreate(PETSC_COMM_SELF, &sfcell));
    216:     PetscCall(PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER));
    217:   }
    218: #endif

    220:   PetscCall(DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos));
    221:   PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
    222:   PetscCall(DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos));

    224:   if (error_check) PetscCall(DMSwarmGetSize(dm, &npointsg));
    225:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    226:   PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
    227:   PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    228:   PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));

    230:   for (p = 0; p < npoints; p++) {
    231:     p_cellid[p] = LA_sfcell[p].index;
    232:     rankval[p]  = rank;
    233:   }
    234:   PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    235:   PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
    236:   PetscCall(PetscSFDestroy(&sfcell));

    238:   if (size > 1) {
    239:     PetscCall(DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration));
    240:   } else {
    241:     DMSwarmDataField PField;
    242:     PetscInt         npoints_curr;

    244:     /* remove points which the domain */
    245:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
    246:     PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));

    248:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
    249:     for (p = 0; p < npoints_curr; p++) {
    250:       if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
    251:         /* kill point */
    252:         PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
    253:         PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL)); /* you need to update npoints as the list size decreases! */
    254:         PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));          /* update date point in case realloc performed */
    255:         p--;                                                                        /* check replacement point */
    256:       }
    257:     }
    258:     PetscCall(DMSwarmGetLocalSize(dm, &npoints_prior_migration));
    259:   }

    261:   /* locate points newly received */
    262:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));

    264: #if 0
    265:   { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
    266:     PetscScalar *LA_coor;
    267:     PetscInt bs;
    268:     DMSwarmDataField PField;

    270:     PetscCall(DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor));
    271:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos));
    272:     PetscCall(DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell));

    274:     PetscCall(VecDestroy(&pos));
    275:     PetscCall(DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor));

    277:     PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
    278:     PetscCall(DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));
    279:     for (p=0; p<npoints2; p++) {
    280:       rankval[p] = LA_sfcell[p].index;
    281:     }
    282:     PetscCall(PetscSFDestroy(&sfcell));
    283:     PetscCall(DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval));

    285:     /* remove points which left processor */
    286:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField));
    287:     PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval));

    289:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL));
    290:     for (p=0; p<npoints2; p++) {
    291:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
    292:         /* kill point */
    293:         PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db,p));
    294:         PetscCall(DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL)); /* you need to update npoints as the list size decreases! */
    295:         PetscCall(DMSwarmDataFieldGetEntries(PField,(void**)&rankval)); /* update date point increase realloc performed */
    296:         p--; /* check replacement point */
    297:       }
    298:     }
    299:   }
    300: #endif

    302:   { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
    303:     PetscScalar     *LA_coor;
    304:     PetscInt         npoints_from_neighbours, bs;
    305:     DMSwarmDataField PField;

    307:     npoints_from_neighbours = npoints2 - npoints_prior_migration;

    309:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor));
    310:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, bs, bs * npoints_from_neighbours, (const PetscScalar *)&LA_coor[bs * npoints_prior_migration], &pos));

    312:     PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));

    314:     PetscCall(VecDestroy(&pos));
    315:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&LA_coor));

    317:     PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
    318:     PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    319:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
    320:     for (p = 0; p < npoints_from_neighbours; p++) {
    321:       rankval[npoints_prior_migration + p]  = rank;
    322:       p_cellid[npoints_prior_migration + p] = LA_sfcell[p].index;
    323:     }

    325:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&p_cellid));
    326:     PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    327:     PetscCall(PetscSFDestroy(&sfcell));

    329:     /* remove points which left processor */
    330:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, DMSwarmPICField_cellid, &PField));
    331:     PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));

    333:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
    334:     for (p = npoints_prior_migration; p < npoints2; p++) {
    335:       if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
    336:         /* kill point */
    337:         PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
    338:         PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL)); /* you need to update npoints as the list size decreases! */
    339:         PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));      /* update date point in case realloc performed */
    340:         p--;                                                                    /* check replacement point */
    341:       }
    342:     }
    343:   }

    345:   /* check for error on removed points */
    346:   if (error_check) {
    347:     PetscCall(DMSwarmGetSize(dm, &npoints2g));
    348:     PetscCheck(npointsg == npoints2g, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points from the DMSwarm must remain constant during migration (initial %" PetscInt_FMT " - final %" PetscInt_FMT ")", npointsg, npoints2g);
    349:   }
    350:   PetscFunctionReturn(PETSC_SUCCESS);
    351: }

    353: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm, PetscBool remove_sent_points)
    354: {
    355:   PetscFunctionBegin;
    356:   PetscFunctionReturn(PETSC_SUCCESS);
    357: }

    359: /*
    360:  Redundant as this assumes points can only be sent to a single rank
    361: */
    362: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize)
    363: {
    364:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
    365:   DMSwarmDataEx de;
    366:   PetscInt      p, npoints, *rankval, n_points_recv;
    367:   PetscMPIInt   rank, nrank, negrank;
    368:   void         *point_buffer, *recv_points;
    369:   size_t        sizeof_dmswarm_point;

    371:   PetscFunctionBegin;
    372:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    373:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    374:   *globalsize = npoints;
    375:   PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    376:   PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
    377:   PetscCall(DMSwarmDataExTopologyInitialize(de));
    378:   for (p = 0; p < npoints; p++) {
    379:     negrank = (PetscMPIInt)rankval[p];
    380:     if (negrank < 0) {
    381:       nrank = -negrank - 1;
    382:       PetscCall(DMSwarmDataExTopologyAddNeighbour(de, nrank));
    383:     }
    384:   }
    385:   PetscCall(DMSwarmDataExTopologyFinalize(de));
    386:   PetscCall(DMSwarmDataExInitializeSendCount(de));
    387:   for (p = 0; p < npoints; p++) {
    388:     negrank = (PetscMPIInt)rankval[p];
    389:     if (negrank < 0) {
    390:       nrank = -negrank - 1;
    391:       PetscCall(DMSwarmDataExAddToSendCount(de, nrank, 1));
    392:     }
    393:   }
    394:   PetscCall(DMSwarmDataExFinalizeSendCount(de));
    395:   PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
    396:   PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
    397:   for (p = 0; p < npoints; p++) {
    398:     negrank = (PetscMPIInt)rankval[p];
    399:     if (negrank < 0) {
    400:       nrank      = -negrank - 1;
    401:       rankval[p] = nrank;
    402:       /* copy point into buffer */
    403:       PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
    404:       /* insert point buffer into DMSwarmDataExchanger */
    405:       PetscCall(DMSwarmDataExPackData(de, nrank, 1, point_buffer));
    406:       rankval[p] = negrank;
    407:     }
    408:   }
    409:   PetscCall(DMSwarmDataExPackFinalize(de));
    410:   PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    411:   PetscCall(DMSwarmDataExBegin(de));
    412:   PetscCall(DMSwarmDataExEnd(de));
    413:   PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
    414:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    415:   PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
    416:   for (p = 0; p < n_points_recv; p++) {
    417:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

    419:     PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
    420:   }
    421:   PetscCall(DMSwarmDataExView(de));
    422:   PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
    423:   PetscCall(DMSwarmDataExDestroy(de));
    424:   PetscFunctionReturn(PETSC_SUCCESS);
    425: }

    427: typedef struct {
    428:   PetscMPIInt owner_rank;
    429:   PetscReal   min[3], max[3];
    430: } CollectBBox;

    432: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm, PetscInt *globalsize)
    433: {
    434:   DM_Swarm          *swarm = (DM_Swarm *)dm->data;
    435:   DMSwarmDataEx      de;
    436:   PetscInt           p, pk, npoints, *rankval, n_points_recv, n_bbox_recv, dim, neighbour_cells;
    437:   PetscMPIInt        rank, nrank;
    438:   void              *point_buffer, *recv_points;
    439:   size_t             sizeof_dmswarm_point, sizeof_bbox_ctx;
    440:   PetscBool          isdmda;
    441:   CollectBBox       *bbox, *recv_bbox;
    442:   const PetscMPIInt *dmneighborranks;
    443:   DM                 dmcell;

    445:   PetscFunctionBegin;
    446:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

    448:   PetscCall(DMSwarmGetCellDM(dm, &dmcell));
    449:   PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
    450:   isdmda = PETSC_FALSE;
    451:   PetscCall(PetscObjectTypeCompare((PetscObject)dmcell, DMDA, &isdmda));
    452:   PetscCheck(isdmda, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only DMDA support for CollectBoundingBox");

    454:   PetscCall(DMGetDimension(dm, &dim));
    455:   sizeof_bbox_ctx = sizeof(CollectBBox);
    456:   PetscCall(PetscMalloc1(1, &bbox));
    457:   bbox->owner_rank = rank;

    459:   /* compute the bounding box based on the overlapping / stenctil size */
    460:   {
    461:     Vec lcoor;

    463:     PetscCall(DMGetCoordinatesLocal(dmcell, &lcoor));
    464:     if (dim >= 1) {
    465:       PetscCall(VecStrideMin(lcoor, 0, NULL, &bbox->min[0]));
    466:       PetscCall(VecStrideMax(lcoor, 0, NULL, &bbox->max[0]));
    467:     }
    468:     if (dim >= 2) {
    469:       PetscCall(VecStrideMin(lcoor, 1, NULL, &bbox->min[1]));
    470:       PetscCall(VecStrideMax(lcoor, 1, NULL, &bbox->max[1]));
    471:     }
    472:     if (dim == 3) {
    473:       PetscCall(VecStrideMin(lcoor, 2, NULL, &bbox->min[2]));
    474:       PetscCall(VecStrideMax(lcoor, 2, NULL, &bbox->max[2]));
    475:     }
    476:   }
    477:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    478:   *globalsize = npoints;
    479:   PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    480:   PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
    481:   /* use DMDA neighbours */
    482:   PetscCall(DMDAGetNeighbors(dmcell, &dmneighborranks));
    483:   if (dim == 1) {
    484:     neighbour_cells = 3;
    485:   } else if (dim == 2) {
    486:     neighbour_cells = 9;
    487:   } else {
    488:     neighbour_cells = 27;
    489:   }
    490:   PetscCall(DMSwarmDataExTopologyInitialize(de));
    491:   for (p = 0; p < neighbour_cells; p++) {
    492:     if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, dmneighborranks[p]));
    493:   }
    494:   PetscCall(DMSwarmDataExTopologyFinalize(de));
    495:   PetscCall(DMSwarmDataExInitializeSendCount(de));
    496:   for (p = 0; p < neighbour_cells; p++) {
    497:     if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) PetscCall(DMSwarmDataExAddToSendCount(de, dmneighborranks[p], 1));
    498:   }
    499:   PetscCall(DMSwarmDataExFinalizeSendCount(de));
    500:   /* send bounding boxes */
    501:   PetscCall(DMSwarmDataExPackInitialize(de, sizeof_bbox_ctx));
    502:   for (p = 0; p < neighbour_cells; p++) {
    503:     nrank = dmneighborranks[p];
    504:     if ((nrank >= 0) && (nrank != rank)) {
    505:       /* insert bbox buffer into DMSwarmDataExchanger */
    506:       PetscCall(DMSwarmDataExPackData(de, nrank, 1, bbox));
    507:     }
    508:   }
    509:   PetscCall(DMSwarmDataExPackFinalize(de));
    510:   /* recv bounding boxes */
    511:   PetscCall(DMSwarmDataExBegin(de));
    512:   PetscCall(DMSwarmDataExEnd(de));
    513:   PetscCall(DMSwarmDataExGetRecvData(de, &n_bbox_recv, (void **)&recv_bbox));
    514:   /*  Wrong, should not be using PETSC_COMM_WORLD */
    515:   for (p = 0; p < n_bbox_recv; p++) {
    516:     PetscCall(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, (double)recv_bbox[p].min[0], (double)recv_bbox[p].max[0], (double)recv_bbox[p].min[1],
    517:                                       (double)recv_bbox[p].max[1]));
    518:   }
    519:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, stdout));
    520:   /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
    521:   PetscCall(DMSwarmDataExInitializeSendCount(de));
    522:   for (pk = 0; pk < n_bbox_recv; pk++) {
    523:     PetscReal *array_x, *array_y;

    525:     PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
    526:     PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
    527:     for (p = 0; p < npoints; p++) {
    528:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
    529:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) PetscCall(DMSwarmDataExAddToSendCount(de, recv_bbox[pk].owner_rank, 1));
    530:       }
    531:     }
    532:     PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
    533:     PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
    534:   }
    535:   PetscCall(DMSwarmDataExFinalizeSendCount(de));
    536:   PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
    537:   PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
    538:   for (pk = 0; pk < n_bbox_recv; pk++) {
    539:     PetscReal *array_x, *array_y;

    541:     PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
    542:     PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));
    543:     for (p = 0; p < npoints; p++) {
    544:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
    545:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
    546:           /* copy point into buffer */
    547:           PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, p, point_buffer));
    548:           /* insert point buffer into DMSwarmDataExchanger */
    549:           PetscCall(DMSwarmDataExPackData(de, recv_bbox[pk].owner_rank, 1, point_buffer));
    550:         }
    551:       }
    552:     }
    553:     PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
    554:     PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));
    555:   }
    556:   PetscCall(DMSwarmDataExPackFinalize(de));
    557:   PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    558:   PetscCall(DMSwarmDataExBegin(de));
    559:   PetscCall(DMSwarmDataExEnd(de));
    560:   PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
    561:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    562:   PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
    563:   for (p = 0; p < n_points_recv; p++) {
    564:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

    566:     PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
    567:   }
    568:   PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
    569:   PetscCall(PetscFree(bbox));
    570:   PetscCall(DMSwarmDataExView(de));
    571:   PetscCall(DMSwarmDataExDestroy(de));
    572:   PetscFunctionReturn(PETSC_SUCCESS);
    573: }

    575: /* General collection when no order, or neighbour information is provided */
    576: /*
    577:  User provides context and collect() method
    578:  Broadcast user context

    580:  for each context / rank {
    581:    collect(swarm,context,n,list)
    582:  }
    583: */
    584: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm, PetscErrorCode (*collect)(DM, void *, PetscInt *, PetscInt **), size_t ctx_size, void *ctx, PetscInt *globalsize)
    585: {
    586:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
    587:   DMSwarmDataEx de;
    588:   PetscInt      p, npoints, n_points_recv;
    589:   PetscMPIInt   size, rank, len;
    590:   void         *point_buffer, *recv_points;
    591:   void         *ctxlist;
    592:   PetscInt     *n2collect, **collectlist;
    593:   size_t        sizeof_dmswarm_point;

    595:   PetscFunctionBegin;
    596:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    597:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    598:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    599:   *globalsize = npoints;
    600:   /* Broadcast user context */
    601:   PetscCall(PetscMPIIntCast(ctx_size, &len));
    602:   PetscCall(PetscMalloc(ctx_size * size, &ctxlist));
    603:   PetscCallMPI(MPI_Allgather(ctx, len, MPI_CHAR, ctxlist, len, MPI_CHAR, PetscObjectComm((PetscObject)dm)));
    604:   PetscCall(PetscMalloc1(size, &n2collect));
    605:   PetscCall(PetscMalloc1(size, &collectlist));
    606:   for (PetscMPIInt r = 0; r < size; r++) {
    607:     PetscInt  _n2collect;
    608:     PetscInt *_collectlist;
    609:     void     *_ctx_r;

    611:     _n2collect   = 0;
    612:     _collectlist = NULL;
    613:     if (r != rank) { /* don't collect data from yourself */
    614:       _ctx_r = (void *)((char *)ctxlist + r * ctx_size);
    615:       PetscCall(collect(dm, _ctx_r, &_n2collect, &_collectlist));
    616:     }
    617:     n2collect[r]   = _n2collect;
    618:     collectlist[r] = _collectlist;
    619:   }
    620:   PetscCall(DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm), 0, &de));
    621:   /* Define topology */
    622:   PetscCall(DMSwarmDataExTopologyInitialize(de));
    623:   for (PetscMPIInt r = 0; r < size; r++) {
    624:     if (n2collect[r] > 0) PetscCall(DMSwarmDataExTopologyAddNeighbour(de, (PetscMPIInt)r));
    625:   }
    626:   PetscCall(DMSwarmDataExTopologyFinalize(de));
    627:   /* Define send counts */
    628:   PetscCall(DMSwarmDataExInitializeSendCount(de));
    629:   for (PetscMPIInt r = 0; r < size; r++) {
    630:     if (n2collect[r] > 0) PetscCall(DMSwarmDataExAddToSendCount(de, r, n2collect[r]));
    631:   }
    632:   PetscCall(DMSwarmDataExFinalizeSendCount(de));
    633:   /* Pack data */
    634:   PetscCall(DMSwarmDataBucketCreatePackedArray(swarm->db, &sizeof_dmswarm_point, &point_buffer));
    635:   PetscCall(DMSwarmDataExPackInitialize(de, sizeof_dmswarm_point));
    636:   for (PetscMPIInt r = 0; r < size; r++) {
    637:     for (p = 0; p < n2collect[r]; p++) {
    638:       PetscCall(DMSwarmDataBucketFillPackedArray(swarm->db, collectlist[r][p], point_buffer));
    639:       /* insert point buffer into the data exchanger */
    640:       PetscCall(DMSwarmDataExPackData(de, r, 1, point_buffer));
    641:     }
    642:   }
    643:   PetscCall(DMSwarmDataExPackFinalize(de));
    644:   /* Scatter */
    645:   PetscCall(DMSwarmDataExBegin(de));
    646:   PetscCall(DMSwarmDataExEnd(de));
    647:   /* Collect data in DMSwarm container */
    648:   PetscCall(DMSwarmDataExGetRecvData(de, &n_points_recv, (void **)&recv_points));
    649:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    650:   PetscCall(DMSwarmDataBucketSetSizes(swarm->db, npoints + n_points_recv, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
    651:   for (p = 0; p < n_points_recv; p++) {
    652:     void *data_p = (void *)((char *)recv_points + p * sizeof_dmswarm_point);

    654:     PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
    655:   }
    656:   /* Release memory */
    657:   for (PetscMPIInt r = 0; r < size; r++) {
    658:     if (collectlist[r]) PetscCall(PetscFree(collectlist[r]));
    659:   }
    660:   PetscCall(PetscFree(collectlist));
    661:   PetscCall(PetscFree(n2collect));
    662:   PetscCall(PetscFree(ctxlist));
    663:   PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
    664:   PetscCall(DMSwarmDataExView(de));
    665:   PetscCall(DMSwarmDataExDestroy(de));
    666:   PetscFunctionReturn(PETSC_SUCCESS);
    667: }

    669: /*@
    670:   DMSwarmGetMigrateType - Get the style of point migration

    672:   Logically Collective

    674:   Input Parameter:
    675: . dm - the `DMSWARM`

    677:   Output Parameter:
    678: . mtype - The migration type, see `DMSwarmMigrateType`

    680:   Level: intermediate

    682: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmMigrate()`
    683: @*/
    684: PetscErrorCode DMSwarmGetMigrateType(DM dm, DMSwarmMigrateType *mtype)
    685: {
    686:   PetscFunctionBegin;
    688:   PetscAssertPointer(mtype, 2);
    689:   *mtype = ((DM_Swarm *)dm->data)->migrate_type;
    690:   PetscFunctionReturn(PETSC_SUCCESS);
    691: }

    693: /*@
    694:   DMSwarmSetMigrateType - Set the style of point migration

    696:   Logically Collective

    698:   Input Parameters:
    699: + dm    - the `DMSWARM`
    700: - mtype - The migration type, see `DMSwarmMigrateType`

    702:   Level: intermediate

    704: .seealso: `DM`, `DMSWARM`, `DMSwarmMigrateType`, `DMSwarmGetMigrateType()`, `DMSwarmMigrate()`
    705: @*/
    706: PetscErrorCode DMSwarmSetMigrateType(DM dm, DMSwarmMigrateType mtype)
    707: {
    708:   PetscFunctionBegin;
    711:   ((DM_Swarm *)dm->data)->migrate_type = mtype;
    712:   PetscFunctionReturn(PETSC_SUCCESS);
    713: }