Actual source code: plexcheckinterface.c

      1: #include <petsc/private/dmpleximpl.h>

      3: /* TODO PetscArrayExchangeBegin/End */
      4: /* TODO blocksize */
      5: /* TODO move to API ? */
      6: static PetscErrorCode ExchangeArrayByRank_Private(PetscObject obj, MPI_Datatype dt, PetscInt nsranks, const PetscMPIInt sranks[], PetscInt ssize[], const void *sarr[], PetscInt nrranks, const PetscMPIInt rranks[], PetscInt *rsize_out[], void **rarr_out[])
      7: {
      8:   PetscInt r;
      9:   PetscInt *rsize;
     10:   void **rarr;
     11:   MPI_Request *sreq, *rreq;
     12:   PetscMPIInt tag, unitsize;
     13:   MPI_Comm comm;

     17:   MPI_Type_size(dt, &unitsize);
     18:   PetscObjectGetComm(obj, &comm);
     19:   PetscMalloc2(nrranks, &rsize, nrranks, &rarr);
     20:   PetscMalloc2(nrranks, &rreq, nsranks, &sreq);
     21:   /* exchange array size */
     22:   PetscObjectGetNewTag(obj,&tag);
     23:   for (r=0; r<nrranks; r++) {
     24:     MPI_Irecv(&rsize[r], 1, MPIU_INT, rranks[r], tag, comm, &rreq[r]);
     25:   }
     26:   for (r=0; r<nsranks; r++) {
     27:     MPI_Isend(&ssize[r], 1, MPIU_INT, sranks[r], tag, comm, &sreq[r]);
     28:   }
     29:   MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE);
     30:   MPI_Waitall(nsranks, sreq, MPI_STATUSES_IGNORE);
     31:   /* exchange array */
     32:   PetscObjectGetNewTag(obj,&tag);
     33:   for (r=0; r<nrranks; r++) {
     34:     PetscMalloc(rsize[r]*unitsize, &rarr[r]);
     35:     MPI_Irecv(rarr[r], rsize[r], dt, rranks[r], tag, comm, &rreq[r]);
     36:   }
     37:   for (r=0; r<nsranks; r++) {
     38:     MPI_Isend(sarr[r], ssize[r], dt, sranks[r], tag, comm, &sreq[r]);
     39:   }
     40:   MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE);
     41:   MPI_Waitall(nsranks, sreq, MPI_STATUSES_IGNORE);
     42:   PetscFree2(rreq, sreq);
     43:   *rsize_out = rsize;
     44:   *rarr_out = rarr;
     45:   return(0);
     46: }

     48: /* TODO VecExchangeBegin/End */
     49: /* TODO move to API ? */
     50: static PetscErrorCode ExchangeVecByRank_Private(PetscObject obj, PetscInt nsranks, const PetscMPIInt sranks[], Vec svecs[], PetscInt nrranks, const PetscMPIInt rranks[], Vec *rvecs[])
     51: {
     52:   PetscInt r;
     53:   PetscInt *ssize, *rsize;
     54:   PetscScalar **rarr;
     55:   const PetscScalar **sarr;
     56:   Vec *rvecs_;
     57:   MPI_Request *sreq, *rreq;

     61:   PetscMalloc4(nsranks, &ssize, nsranks, &sarr, nrranks, &rreq, nsranks, &sreq);
     62:   for (r=0; r<nsranks; r++) {
     63:     VecGetLocalSize(svecs[r], &ssize[r]);
     64:     VecGetArrayRead(svecs[r], &sarr[r]);
     65:   }
     66:   ExchangeArrayByRank_Private(obj, MPIU_SCALAR, nsranks, sranks, ssize, (const void**)sarr, nrranks, rranks, &rsize, (void***)&rarr);
     67:   PetscMalloc1(nrranks, &rvecs_);
     68:   for (r=0; r<nrranks; r++) {
     69:     /* set array in two steps to mimic PETSC_OWN_POINTER */
     70:     VecCreateSeqWithArray(PETSC_COMM_SELF, 1, rsize[r], NULL, &rvecs_[r]);
     71:     VecReplaceArray(rvecs_[r], rarr[r]);
     72:   }
     73:   for (r=0; r<nsranks; r++) {
     74:     VecRestoreArrayRead(svecs[r], &sarr[r]);
     75:   }
     76:   PetscFree2(rsize, rarr);
     77:   PetscFree4(ssize, sarr, rreq, sreq);
     78:   *rvecs = rvecs_;
     79:   return(0);
     80: }

     82: static PetscErrorCode SortByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[])
     83: {
     84:   PetscInt            nleaves;
     85:   PetscInt            nranks;
     86:   const PetscMPIInt   *ranks;
     87:   const PetscInt      *roffset, *rmine, *rremote;
     88:   PetscInt            n, o, r;
     89:   PetscErrorCode      ierr;

     92:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
     93:   nleaves = roffset[nranks];
     94:   PetscMalloc2(nleaves, rmine1, nleaves, rremote1);
     95:   for (r=0; r<nranks; r++) {
     96:     /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote
     97:        - to unify order with the other side */
     98:     o = roffset[r];
     99:     n = roffset[r+1] - o;
    100:     PetscArraycpy(&(*rmine1)[o], &rmine[o], n);
    101:     PetscArraycpy(&(*rremote1)[o], &rremote[o], n);
    102:     PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);
    103:   }
    104:   return(0);
    105: }

    107: static PetscErrorCode GetRecursiveConeCoordinatesPerRank_Private(DM dm, PetscSF sf, PetscInt rmine[], Vec *coordinatesPerRank[])
    108: {
    109:   IS                  pointsPerRank, conesPerRank;
    110:   PetscInt            nranks;
    111:   const PetscMPIInt   *ranks;
    112:   const PetscInt      *roffset;
    113:   PetscInt            n, o, r;
    114:   PetscErrorCode      ierr;

    117:   DMGetCoordinatesLocalSetUp(dm);
    118:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);
    119:   PetscMalloc1(nranks, coordinatesPerRank);
    120:   for (r=0; r<nranks; r++) {
    121:     o = roffset[r];
    122:     n = roffset[r+1] - o;
    123:     ISCreateGeneral(PETSC_COMM_SELF, n, &rmine[o], PETSC_USE_POINTER, &pointsPerRank);
    124:     DMPlexGetConeRecursiveVertices(dm, pointsPerRank, &conesPerRank);
    125:     DMGetCoordinatesLocalTuple(dm, conesPerRank, NULL, &(*coordinatesPerRank)[r]);
    126:     ISDestroy(&pointsPerRank);
    127:     ISDestroy(&conesPerRank);
    128:   }
    129:   return(0);
    130: }

    132: static PetscErrorCode PetscSFComputeMultiRootOriginalNumberingByRank_Private(PetscSF sf, PetscSF imsf, PetscInt *irmine1[])
    133: {
    134:   PetscInt            *mRootsOrigNumbering;
    135:   PetscInt            nileaves, niranks;
    136:   const PetscInt      *iroffset, *irmine, *degree;
    137:   PetscInt            i, n, o, r;
    138:   PetscErrorCode      ierr;

    141:   PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL);
    142:   PetscSFGetRootRanks(imsf, &niranks, NULL, &iroffset, &irmine, NULL);
    143:   if (PetscUnlikely(nileaves != iroffset[niranks])) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nileaves != iroffset[niranks])");
    144:   PetscSFComputeDegreeBegin(sf, &degree);
    145:   PetscSFComputeDegreeEnd(sf, &degree);
    146:   PetscSFComputeMultiRootOriginalNumbering(sf, degree, NULL, &mRootsOrigNumbering);
    147:   PetscMalloc1(nileaves, irmine1);
    148:   for (r=0; r<niranks; r++) {
    149:     o = iroffset[r];
    150:     n = iroffset[r+1] - o;
    151:     for (i=0; i<n; i++) (*irmine1)[o+i] = mRootsOrigNumbering[irmine[o+i]];
    152:   }
    153:   PetscFree(mRootsOrigNumbering);
    154:   return(0);
    155: }

    157: /*@
    158:   DMPlexCheckInterfaceCones - Check that points on inter-partition interfaces have conforming order of cone points.

    160:   Input Parameters:
    161: . dm - The DMPlex object

    163:   Notes:
    164:   For example, if there is an edge (rank,index)=(0,2) connecting points cone(0,2)=[(0,0),(0,1)] in this order, and the point SF containts connections 0 <- (1,0), 1 <- (1,1) and 2 <- (1,2),
    165:   then this check would pass if the edge (1,2) has cone(1,2)=[(1,0),(1,1)]. By contrast, if cone(1,2)=[(1,1),(1,0)], then this check would fail.

    167:   This is mainly intended for debugging/testing purposes. Does not check cone orientation, for this purpose use DMPlexCheckFaces().

    169:   For the complete list of DMPlexCheck* functions, see DMSetFromOptions().

    171:   Developer Note:
    172:   Interface cones are expanded into vertices and then their coordinates are compared.

    174:   Level: developer

    176: .seealso: DMPlexGetCone(), DMPlexGetConeSize(), DMGetPointSF(), DMGetCoordinates(), DMSetFromOptions()
    177: @*/
    178: PetscErrorCode DMPlexCheckInterfaceCones(DM dm)
    179: {
    180:   PetscSF             sf;
    181:   PetscInt            nleaves, nranks, nroots;
    182:   const PetscInt      *mine, *roffset, *rmine, *rremote;
    183:   const PetscSFNode   *remote;
    184:   const PetscMPIInt   *ranks;
    185:   PetscSF             msf, imsf;
    186:   PetscInt            nileaves, niranks;
    187:   const PetscMPIInt   *iranks;
    188:   const PetscInt      *iroffset, *irmine, *irremote;
    189:   PetscInt            *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
    190:   PetscInt            *mine_orig_numbering;
    191:   Vec                 *sntCoordinatesPerRank;
    192:   Vec                 *refCoordinatesPerRank;
    193:   Vec                 *recCoordinatesPerRank=NULL;
    194:   PetscInt            r;
    195:   PetscMPIInt         commsize, myrank;
    196:   PetscBool           same;
    197:   PetscBool           verbose=PETSC_FALSE;
    198:   MPI_Comm            comm;
    199:   PetscErrorCode      ierr;

    203:   PetscObjectGetComm((PetscObject)dm, &comm);
    204:   MPI_Comm_rank(comm, &myrank);
    205:   MPI_Comm_size(comm, &commsize);
    206:   if (commsize < 2) return(0);
    207:   DMGetPointSF(dm, &sf);
    208:   if (!sf) return(0);
    209:   PetscSFGetGraph(sf, &nroots, &nleaves, &mine, &remote);
    210:   if (nroots < 0) return(0);
    211:   if (!dm->coordinates && !dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM coordinates must be set");
    212:   PetscSFSetUp(sf);
    213:   PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);

    215:   /* Expand sent cones per rank */
    216:   SortByRemote_Private(sf, &rmine1, &rremote1);
    217:   GetRecursiveConeCoordinatesPerRank_Private(dm, sf, rmine1, &sntCoordinatesPerRank);

    219:   /* Create inverse SF */
    220:   PetscSFGetMultiSF(sf,&msf);
    221:   PetscSFCreateInverseSF(msf,&imsf);
    222:   PetscSFSetUp(imsf);
    223:   PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL);
    224:   PetscSFGetRootRanks(imsf, &niranks, &iranks, &iroffset, &irmine, &irremote);

    226:   /* Compute original numbering of multi-roots (referenced points) */
    227:   PetscSFComputeMultiRootOriginalNumberingByRank_Private(sf, imsf, &mine_orig_numbering);

    229:   /* Expand coordinates of the referred cones per rank */
    230:   GetRecursiveConeCoordinatesPerRank_Private(dm, imsf, mine_orig_numbering, &refCoordinatesPerRank);

    232:   /* Send the coordinates */
    233:   ExchangeVecByRank_Private((PetscObject)sf, nranks, ranks, sntCoordinatesPerRank, niranks, iranks, &recCoordinatesPerRank);

    235:   /* verbose output */
    236:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_check_cones_conform_on_interfaces_verbose", &verbose, NULL);
    237:   if (verbose) {
    238:     PetscViewer sv, v = PETSC_VIEWER_STDOUT_WORLD;
    239:     PetscViewerASCIIPrintf(v, "============\nDMPlexCheckInterfaceCones output\n============\n");
    240:     PetscViewerASCIIPushSynchronized(v);
    241:     PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", myrank);
    242:     for (r=0; r<nranks; r++) {
    243:       PetscViewerASCIISynchronizedPrintf(v, "  r=%D ranks[r]=%d sntCoordinatesPerRank[r]:\n", r, ranks[r]);
    244:       PetscViewerASCIIPushTab(v);
    245:       PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);
    246:       VecView(sntCoordinatesPerRank[r], sv);
    247:       PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);
    248:       PetscViewerASCIIPopTab(v);
    249:     }
    250:     PetscViewerASCIISynchronizedPrintf(v, "  ----------\n");
    251:     for (r=0; r<niranks; r++) {
    252:       PetscViewerASCIISynchronizedPrintf(v, "  r=%D iranks[r]=%d refCoordinatesPerRank[r]:\n", r, iranks[r]);
    253:       PetscViewerASCIIPushTab(v);
    254:       PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);
    255:       VecView(refCoordinatesPerRank[r], sv);
    256:       PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);
    257:       PetscViewerASCIIPopTab(v);
    258:     }
    259:     PetscViewerASCIISynchronizedPrintf(v, "  ----------\n");
    260:     for (r=0; r<niranks; r++) {
    261:       PetscViewerASCIISynchronizedPrintf(v, "  r=%D iranks[r]=%d recCoordinatesPerRank[r]:\n", r, iranks[r]);
    262:       PetscViewerASCIIPushTab(v);
    263:       PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);
    264:       VecView(recCoordinatesPerRank[r], sv);
    265:       PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);
    266:       PetscViewerASCIIPopTab(v);
    267:     }
    268:   }

    270:   /* Compare recCoordinatesPerRank with refCoordinatesPerRank */
    271:   for (r=0; r<niranks; r++) {
    272:     VecEqual(refCoordinatesPerRank[r], recCoordinatesPerRank[r], &same);
    273:     if (!same) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interface cones do not conform for remote rank %d", iranks[r]);
    274:   }

    276:   /* destroy sent stuff */
    277:   for (r=0; r<nranks; r++) {
    278:     VecDestroy(&sntCoordinatesPerRank[r]);
    279:   }
    280:   PetscFree(sntCoordinatesPerRank);
    281:   PetscFree2(rmine1, rremote1);
    282:   PetscSFDestroy(&imsf);

    284:   /* destroy referenced stuff */
    285:   for (r=0; r<niranks; r++) {
    286:     VecDestroy(&refCoordinatesPerRank[r]);
    287:   }
    288:   PetscFree(refCoordinatesPerRank);
    289:   PetscFree(mine_orig_numbering);

    291:   /* destroy received stuff */
    292:   for (r=0; r<niranks; r++) {
    293:     VecDestroy(&recCoordinatesPerRank[r]);
    294:   }
    295:   PetscFree(recCoordinatesPerRank);
    296:   return(0);
    297: }