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