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:     PetscCall(PetscMPIIntCast(rankval[p], &nrank));
     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:     PetscCall(PetscMPIIntCast(rankval[p], &nrank));
     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:     PetscCall(PetscMPIIntCast(rankval[p], &nrank));
     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:       PetscCall(PetscMPIIntCast(rankval[p], &nrank));
     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, &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:   DMSwarmCellDM      celldm;
    101:   DMSwarmDataEx      de;
    102:   PetscInt           r, p, npoints, *p_cellid, n_points_recv;
    103:   PetscMPIInt        rank, _rank;
    104:   const PetscMPIInt *neighbourranks;
    105:   void              *point_buffer, *recv_points;
    106:   size_t             sizeof_dmswarm_point;
    107:   PetscInt           nneighbors;
    108:   PetscMPIInt        mynneigh, *myneigh;
    109:   const char        *cellid;

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

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

    177:     PetscCall(DMSwarmDataBucketInsertPackedArray(swarm->db, npoints + p, data_p));
    178:   }
    179:   PetscCall(DMSwarmDataBucketDestroyPackedArray(swarm->db, &point_buffer));
    180:   PetscCall(DMSwarmDataExDestroy(de));
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm, PetscBool remove_sent_points)
    185: {
    186:   DM_Swarm          *swarm = (DM_Swarm *)dm->data;
    187:   DMSwarmCellDM      celldm;
    188:   PetscInt           p, npoints, npointsg = 0, npoints2, npoints2g, *rankval, *p_cellid, npoints_prior_migration, Nfc;
    189:   PetscSF            sfcell = NULL;
    190:   const PetscSFNode *LA_sfcell;
    191:   DM                 dmcell;
    192:   Vec                pos;
    193:   PetscBool          error_check = swarm->migrate_error_on_missing_point;
    194:   const char       **coordFields;
    195:   PetscMPIInt        size, rank;
    196:   const char        *cellid;

    198:   PetscFunctionBegin;
    199:   PetscCall(DMSwarmGetCellDM(dm, &dmcell));
    200:   PetscCheck(dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid if cell DM provided");
    201:   PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
    202:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));

    204:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    205:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));

    207: #if 1
    208:   {
    209:     PetscInt     npoints_curr, range = 0;
    210:     PetscSFNode *sf_cells;

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

    215:     PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&p_cellid));
    216:     for (p = 0; p < npoints_curr; p++) {
    217:       sf_cells[p].rank  = 0;
    218:       sf_cells[p].index = p_cellid[p];
    219:       if (p_cellid[p] > range) range = p_cellid[p];
    220:     }
    221:     PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));

    223:     /* PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell)); */
    224:     PetscCall(PetscSFCreate(PETSC_COMM_SELF, &sfcell));
    225:     PetscCall(PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER));
    226:   }
    227: #endif

    229:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
    230:   PetscCall(DMSwarmCreateLocalVectorFromFields(dm, Nfc, coordFields, &pos));
    231:   PetscCall(DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell));
    232:   PetscCall(DMSwarmDestroyLocalVectorFromFields(dm, Nfc, coordFields, &pos));

    234:   if (error_check) PetscCall(DMSwarmGetSize(dm, &npointsg));
    235:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
    236:   PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&p_cellid));
    237:   PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    238:   PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));

    240:   for (p = 0; p < npoints; p++) {
    241:     p_cellid[p] = LA_sfcell[p].index;
    242:     rankval[p]  = rank;
    243:   }
    244:   PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
    245:   PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));
    246:   PetscCall(PetscSFDestroy(&sfcell));

    248:   if (size > 1) {
    249:     PetscCall(DMSwarmMigrate_DMNeighborScatter(dm, dmcell, remove_sent_points, &npoints_prior_migration));
    250:   } else {
    251:     DMSwarmDataField PField;
    252:     PetscInt         npoints_curr;

    254:     /* remove points which the domain */
    255:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, cellid, &PField));
    256:     PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));

    258:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL));
    259:     if (remove_sent_points) {
    260:       for (p = 0; p < npoints_curr; p++) {
    261:         if (p_cellid[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
    262:           /* kill point */
    263:           PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, p));
    264:           PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints_curr, NULL, NULL)); /* you need to update npoints as the list size decreases! */
    265:           PetscCall(DMSwarmDataFieldGetEntries(PField, (void **)&p_cellid));          /* update date point in case realloc performed */
    266:           p--;                                                                        /* check replacement point */
    267:         }
    268:       }
    269:     }
    270:     PetscCall(DMSwarmGetLocalSize(dm, &npoints_prior_migration));
    271:   }

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

    276:   { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
    277:     Vec              npos;
    278:     IS               nis;
    279:     PetscInt         npoints_from_neighbours, bs;
    280:     DMSwarmDataField PField;

    282:     npoints_from_neighbours = npoints2 - npoints_prior_migration;

    284:     PetscCall(DMSwarmCreateLocalVectorFromFields(dm, Nfc, coordFields, &pos));
    285:     PetscCall(VecGetBlockSize(pos, &bs));
    286:     PetscCall(ISCreateStride(PETSC_COMM_SELF, npoints_from_neighbours * bs, npoints_prior_migration * bs, 1, &nis));
    287: /*
    288:   Device VecGetSubVector to zero sized subvector triggers
    289:   debug error with mismatching memory types due to the device
    290:   pointer being host memtype without anything to point to in
    291:   Vec"Type"GetArrays(...), and we still need to pass something to
    292:   DMLocatePoints to avoid deadlock
    293: */
    294: #if defined(PETSC_HAVE_DEVICE)
    295:     if (npoints_from_neighbours > 0) {
    296:       PetscCall(VecGetSubVector(pos, nis, &npos));
    297:       PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
    298:       PetscCall(VecRestoreSubVector(pos, nis, &npos));
    299:     } else {
    300:       PetscCall(VecCreate(PETSC_COMM_SELF, &npos));
    301:       PetscCall(VecSetSizes(npos, 0, PETSC_DETERMINE));
    302:       PetscCall(VecSetBlockSize(npos, bs));
    303:       PetscCall(VecSetType(npos, dm->vectype));
    304:       PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
    305:       PetscCall(VecDestroy(&npos));
    306:     }
    307: #else
    308:     PetscCall(VecGetSubVector(pos, nis, &npos));
    309:     PetscCall(DMLocatePoints(dmcell, npos, DM_POINTLOCATION_NONE, &sfcell));
    310:     PetscCall(VecRestoreSubVector(pos, nis, &npos));
    311: #endif
    312:     PetscCall(ISDestroy(&nis));
    313:     PetscCall(DMSwarmDestroyLocalVectorFromFields(dm, Nfc, coordFields, &pos));

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

    323:     PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&p_cellid));
    324:     PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));

    326:     PetscCall(PetscSFDestroy(&sfcell));

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

    332:     PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints2, NULL, NULL));
    333:     if (remove_sent_points) {
    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:     }
    344:   }

    346:   /* check for error on removed points */
    347:   if (error_check) {
    348:     PetscCall(DMSwarmGetSize(dm, &npoints2g));
    349:     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);
    350:   }
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    671:   Logically Collective

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

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

    679:   Level: intermediate

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

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

    695:   Logically Collective

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

    701:   Level: intermediate

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