1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 /* TODO PetscArrayExchangeBegin/End */ 4 /* TODO blocksize */ 5 /* TODO move to API ? */ 6 static PetscErrorCode ExchangeArrayByRank_Private(PetscObject obj, MPI_Datatype dt, PetscMPIInt nsranks, const PetscMPIInt sranks[], PetscInt ssize[], const void *sarr[], PetscMPIInt nrranks, const PetscMPIInt rranks[], PetscInt *rsize_out[], void **rarr_out[]) 7 { 8 PetscMPIInt r; 9 PetscInt *rsize; 10 void **rarr; 11 MPI_Request *sreq, *rreq; 12 PetscMPIInt tag, unitsize; 13 MPI_Comm comm; 14 15 PetscFunctionBegin; 16 PetscCallMPI(MPI_Type_size(dt, &unitsize)); 17 PetscCall(PetscObjectGetComm(obj, &comm)); 18 PetscCall(PetscMalloc2(nrranks, &rsize, nrranks, &rarr)); 19 PetscCall(PetscMalloc2(nrranks, &rreq, nsranks, &sreq)); 20 /* exchange array size */ 21 PetscCall(PetscObjectGetNewTag(obj, &tag)); 22 for (r = 0; r < nrranks; r++) PetscCallMPI(MPIU_Irecv(&rsize[r], 1, MPIU_INT, rranks[r], tag, comm, &rreq[r])); 23 for (r = 0; r < nsranks; r++) PetscCallMPI(MPIU_Isend(&ssize[r], 1, MPIU_INT, sranks[r], tag, comm, &sreq[r])); 24 PetscCallMPI(MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE)); 25 PetscCallMPI(MPI_Waitall(nsranks, sreq, MPI_STATUSES_IGNORE)); 26 /* exchange array */ 27 PetscCall(PetscObjectGetNewTag(obj, &tag)); 28 for (r = 0; r < nrranks; r++) { 29 PetscCall(PetscMalloc(rsize[r] * unitsize, &rarr[r])); 30 PetscCallMPI(MPIU_Irecv(rarr[r], rsize[r], dt, rranks[r], tag, comm, &rreq[r])); 31 } 32 for (r = 0; r < nsranks; r++) PetscCallMPI(MPIU_Isend(sarr[r], ssize[r], dt, sranks[r], tag, comm, &sreq[r])); 33 PetscCallMPI(MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE)); 34 PetscCallMPI(MPI_Waitall(nsranks, sreq, MPI_STATUSES_IGNORE)); 35 PetscCall(PetscFree2(rreq, sreq)); 36 *rsize_out = rsize; 37 *rarr_out = rarr; 38 PetscFunctionReturn(PETSC_SUCCESS); 39 } 40 41 /* TODO VecExchangeBegin/End */ 42 /* TODO move to API ? */ 43 static PetscErrorCode ExchangeVecByRank_Private(PetscObject obj, PetscMPIInt nsranks, const PetscMPIInt sranks[], Vec svecs[], PetscMPIInt nrranks, const PetscMPIInt rranks[], Vec *rvecs[]) 44 { 45 PetscMPIInt r; 46 PetscInt *ssize, *rsize; 47 PetscScalar **rarr; 48 const PetscScalar **sarr; 49 Vec *rvecs_; 50 MPI_Request *sreq, *rreq; 51 52 PetscFunctionBegin; 53 PetscCall(PetscMalloc4(nsranks, &ssize, nsranks, &sarr, nrranks, &rreq, nsranks, &sreq)); 54 for (r = 0; r < nsranks; r++) { 55 PetscCall(VecGetLocalSize(svecs[r], &ssize[r])); 56 PetscCall(VecGetArrayRead(svecs[r], &sarr[r])); 57 } 58 PetscCall(ExchangeArrayByRank_Private(obj, MPIU_SCALAR, nsranks, sranks, ssize, (const void **)sarr, nrranks, rranks, &rsize, (void ***)&rarr)); 59 PetscCall(PetscMalloc1(nrranks, &rvecs_)); 60 for (r = 0; r < nrranks; r++) { 61 /* set array in two steps to mimic PETSC_OWN_POINTER */ 62 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, rsize[r], NULL, &rvecs_[r])); 63 PetscCall(VecReplaceArray(rvecs_[r], rarr[r])); 64 } 65 for (r = 0; r < nsranks; r++) PetscCall(VecRestoreArrayRead(svecs[r], &sarr[r])); 66 PetscCall(PetscFree2(rsize, rarr)); 67 PetscCall(PetscFree4(ssize, sarr, rreq, sreq)); 68 *rvecs = rvecs_; 69 PetscFunctionReturn(PETSC_SUCCESS); 70 } 71 72 static PetscErrorCode SortByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[]) 73 { 74 PetscInt nleaves; 75 PetscMPIInt nranks; 76 const PetscMPIInt *ranks; 77 const PetscInt *roffset, *rmine, *rremote; 78 PetscInt n, o; 79 80 PetscFunctionBegin; 81 PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote)); 82 nleaves = roffset[nranks]; 83 PetscCall(PetscMalloc2(nleaves, rmine1, nleaves, rremote1)); 84 for (PetscMPIInt r = 0; r < nranks; r++) { 85 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote 86 - to unify order with the other side */ 87 o = roffset[r]; 88 n = roffset[r + 1] - o; 89 PetscCall(PetscArraycpy(&(*rmine1)[o], &rmine[o], n)); 90 PetscCall(PetscArraycpy(&(*rremote1)[o], &rremote[o], n)); 91 PetscCall(PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o])); 92 } 93 PetscFunctionReturn(PETSC_SUCCESS); 94 } 95 96 static PetscErrorCode GetRecursiveConeCoordinatesPerRank_Private(DM dm, PetscSF sf, PetscInt rmine[], Vec *coordinatesPerRank[]) 97 { 98 IS pointsPerRank, conesPerRank; 99 PetscMPIInt nranks; 100 const PetscMPIInt *ranks; 101 const PetscInt *roffset; 102 PetscInt n, o; 103 104 PetscFunctionBegin; 105 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 106 PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL)); 107 PetscCall(PetscMalloc1(nranks, coordinatesPerRank)); 108 for (PetscMPIInt r = 0; r < nranks; r++) { 109 o = roffset[r]; 110 n = roffset[r + 1] - o; 111 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, &rmine[o], PETSC_USE_POINTER, &pointsPerRank)); 112 PetscCall(DMPlexGetConeRecursiveVertices(dm, pointsPerRank, &conesPerRank)); 113 PetscCall(DMGetCoordinatesLocalTuple(dm, conesPerRank, NULL, &(*coordinatesPerRank)[r])); 114 PetscCall(ISDestroy(&pointsPerRank)); 115 PetscCall(ISDestroy(&conesPerRank)); 116 } 117 PetscFunctionReturn(PETSC_SUCCESS); 118 } 119 120 static PetscErrorCode PetscSFComputeMultiRootOriginalNumberingByRank_Private(PetscSF sf, PetscSF imsf, PetscInt *irmine1[]) 121 { 122 PetscInt *mRootsOrigNumbering; 123 PetscMPIInt niranks; 124 PetscInt nileaves; 125 const PetscInt *iroffset, *irmine, *degree; 126 PetscInt n, o; 127 128 PetscFunctionBegin; 129 PetscCall(PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL)); 130 PetscCall(PetscSFGetRootRanks(imsf, &niranks, NULL, &iroffset, &irmine, NULL)); 131 PetscCheck(nileaves == iroffset[niranks], PETSC_COMM_SELF, PETSC_ERR_PLIB, "nileaves != iroffset[niranks])"); 132 PetscCall(PetscSFComputeDegreeBegin(sf, °ree)); 133 PetscCall(PetscSFComputeDegreeEnd(sf, °ree)); 134 PetscCall(PetscSFComputeMultiRootOriginalNumbering(sf, degree, NULL, &mRootsOrigNumbering)); 135 PetscCall(PetscMalloc1(nileaves, irmine1)); 136 for (PetscMPIInt r = 0; r < niranks; r++) { 137 o = iroffset[r]; 138 n = iroffset[r + 1] - o; 139 for (PetscInt i = 0; i < n; i++) (*irmine1)[o + i] = mRootsOrigNumbering[irmine[o + i]]; 140 } 141 PetscCall(PetscFree(mRootsOrigNumbering)); 142 PetscFunctionReturn(PETSC_SUCCESS); 143 } 144 145 /*@ 146 DMPlexCheckInterfaceCones - Check that points on inter-partition interfaces have conforming order of cone points. 147 148 Input Parameter: 149 . dm - The `DMPLEX` object 150 151 Level: developer 152 153 Notes: 154 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 `PetscSF` contains connections 0 <- (1,0), 1 <- (1,1) and 2 <- (1,2), 155 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. 156 157 This is mainly intended for debugging/testing purposes. Does not check cone orientation, for this purpose use `DMPlexCheckFaces()`. 158 159 For the complete list of DMPlexCheck* functions, see `DMSetFromOptions()`. 160 161 Developer Notes: 162 Interface cones are expanded into vertices and then their coordinates are compared. 163 164 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexGetCone()`, `DMPlexGetConeSize()`, `DMGetPointSF()`, `DMGetCoordinates()`, `DMSetFromOptions()` 165 @*/ 166 PetscErrorCode DMPlexCheckInterfaceCones(DM dm) 167 { 168 PetscSF sf; 169 PetscInt nleaves, nroots; 170 PetscMPIInt nranks; 171 const PetscInt *mine, *roffset, *rmine, *rremote; 172 const PetscSFNode *remote; 173 const PetscMPIInt *ranks; 174 PetscSF msf, imsf; 175 PetscMPIInt niranks; 176 PetscInt nileaves; 177 const PetscMPIInt *iranks; 178 const PetscInt *iroffset, *irmine, *irremote; 179 PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */ 180 PetscInt *mine_orig_numbering; 181 Vec *sntCoordinatesPerRank; 182 Vec *refCoordinatesPerRank; 183 Vec *recCoordinatesPerRank = NULL; 184 PetscInt r; 185 PetscMPIInt size, rank; 186 PetscBool same; 187 PetscBool verbose = PETSC_FALSE; 188 MPI_Comm comm; 189 190 PetscFunctionBegin; 191 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 192 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 193 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 194 PetscCallMPI(MPI_Comm_size(comm, &size)); 195 if (size < 2) PetscFunctionReturn(PETSC_SUCCESS); 196 PetscCall(DMGetPointSF(dm, &sf)); 197 if (!sf) PetscFunctionReturn(PETSC_SUCCESS); 198 PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, &mine, &remote)); 199 if (nroots < 0) PetscFunctionReturn(PETSC_SUCCESS); 200 PetscCheck(dm->coordinates[0].x || dm->coordinates[0].xl, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM coordinates must be set"); 201 PetscCall(PetscSFSetUp(sf)); 202 PetscCall(PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote)); 203 204 /* Expand sent cones per rank */ 205 PetscCall(SortByRemote_Private(sf, &rmine1, &rremote1)); 206 PetscCall(GetRecursiveConeCoordinatesPerRank_Private(dm, sf, rmine1, &sntCoordinatesPerRank)); 207 208 /* Create inverse SF */ 209 PetscCall(PetscSFGetMultiSF(sf, &msf)); 210 PetscCall(PetscSFCreateInverseSF(msf, &imsf)); 211 PetscCall(PetscSFSetUp(imsf)); 212 PetscCall(PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL)); 213 PetscCall(PetscSFGetRootRanks(imsf, &niranks, &iranks, &iroffset, &irmine, &irremote)); 214 215 /* Compute original numbering of multi-roots (referenced points) */ 216 PetscCall(PetscSFComputeMultiRootOriginalNumberingByRank_Private(sf, imsf, &mine_orig_numbering)); 217 218 /* Expand coordinates of the referred cones per rank */ 219 PetscCall(GetRecursiveConeCoordinatesPerRank_Private(dm, imsf, mine_orig_numbering, &refCoordinatesPerRank)); 220 221 /* Send the coordinates */ 222 PetscCall(ExchangeVecByRank_Private((PetscObject)sf, nranks, ranks, sntCoordinatesPerRank, niranks, iranks, &recCoordinatesPerRank)); 223 224 /* verbose output */ 225 PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_check_cones_conform_on_interfaces_verbose", &verbose, NULL)); 226 if (verbose) { 227 PetscViewer sv, v = PETSC_VIEWER_STDOUT_WORLD; 228 PetscCall(PetscViewerASCIIPrintf(v, "============\nDMPlexCheckInterfaceCones output\n============\n")); 229 PetscCall(PetscViewerASCIIPushSynchronized(v)); 230 PetscCall(PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", rank)); 231 for (r = 0; r < size; r++) { 232 if (r < nranks) { 233 PetscCall(PetscViewerASCIISynchronizedPrintf(v, " r=%" PetscInt_FMT " ranks[r]=%d sntCoordinatesPerRank[r]:\n", r, ranks[r])); 234 PetscCall(PetscViewerASCIIPushTab(v)); 235 PetscCall(PetscViewerGetSubViewer(v, PETSC_COMM_SELF, &sv)); 236 PetscCall(VecView(sntCoordinatesPerRank[r], sv)); 237 PetscCall(PetscViewerRestoreSubViewer(v, PETSC_COMM_SELF, &sv)); 238 PetscCall(PetscViewerASCIIPopTab(v)); 239 } else { 240 PetscCall(PetscViewerGetSubViewer(v, PETSC_COMM_SELF, &sv)); 241 PetscCall(PetscViewerRestoreSubViewer(v, PETSC_COMM_SELF, &sv)); 242 } 243 } 244 PetscCall(PetscViewerASCIISynchronizedPrintf(v, " ----------\n")); 245 for (r = 0; r < size; r++) { 246 if (r < niranks) { 247 PetscCall(PetscViewerASCIISynchronizedPrintf(v, " r=%" PetscInt_FMT " iranks[r]=%d refCoordinatesPerRank[r]:\n", r, iranks[r])); 248 PetscCall(PetscViewerASCIIPushTab(v)); 249 PetscCall(PetscViewerGetSubViewer(v, PETSC_COMM_SELF, &sv)); 250 PetscCall(VecView(refCoordinatesPerRank[r], sv)); 251 PetscCall(PetscViewerRestoreSubViewer(v, PETSC_COMM_SELF, &sv)); 252 PetscCall(PetscViewerASCIIPopTab(v)); 253 } else { 254 PetscCall(PetscViewerGetSubViewer(v, PETSC_COMM_SELF, &sv)); 255 PetscCall(PetscViewerRestoreSubViewer(v, PETSC_COMM_SELF, &sv)); 256 } 257 } 258 PetscCall(PetscViewerASCIISynchronizedPrintf(v, " ----------\n")); 259 for (r = 0; r < size; r++) { 260 if (r < niranks) { 261 PetscCall(PetscViewerASCIISynchronizedPrintf(v, " r=%" PetscInt_FMT " iranks[r]=%d recCoordinatesPerRank[r]:\n", r, iranks[r])); 262 PetscCall(PetscViewerASCIIPushTab(v)); 263 PetscCall(PetscViewerGetSubViewer(v, PETSC_COMM_SELF, &sv)); 264 PetscCall(VecView(recCoordinatesPerRank[r], sv)); 265 PetscCall(PetscViewerRestoreSubViewer(v, PETSC_COMM_SELF, &sv)); 266 PetscCall(PetscViewerASCIIPopTab(v)); 267 } else { 268 PetscCall(PetscViewerGetSubViewer(v, PETSC_COMM_SELF, &sv)); 269 PetscCall(PetscViewerRestoreSubViewer(v, PETSC_COMM_SELF, &sv)); 270 } 271 } 272 PetscCall(PetscViewerASCIIPopSynchronized(v)); 273 } 274 275 /* Compare recCoordinatesPerRank with refCoordinatesPerRank */ 276 for (r = 0; r < niranks; r++) { 277 PetscCall(VecEqual(refCoordinatesPerRank[r], recCoordinatesPerRank[r], &same)); 278 PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_PLIB, "interface cones do not conform for remote rank %d", iranks[r]); 279 } 280 281 /* destroy sent stuff */ 282 for (r = 0; r < nranks; r++) PetscCall(VecDestroy(&sntCoordinatesPerRank[r])); 283 PetscCall(PetscFree(sntCoordinatesPerRank)); 284 PetscCall(PetscFree2(rmine1, rremote1)); 285 PetscCall(PetscSFDestroy(&imsf)); 286 287 /* destroy referenced stuff */ 288 for (r = 0; r < niranks; r++) PetscCall(VecDestroy(&refCoordinatesPerRank[r])); 289 PetscCall(PetscFree(refCoordinatesPerRank)); 290 PetscCall(PetscFree(mine_orig_numbering)); 291 292 /* destroy received stuff */ 293 for (r = 0; r < niranks; r++) PetscCall(VecDestroy(&recCoordinatesPerRank[r])); 294 PetscCall(PetscFree(recCoordinatesPerRank)); 295 PetscFunctionReturn(PETSC_SUCCESS); 296 } 297