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 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 ierr = MPI_Type_size(dt, &unitsize);CHKERRMPI(ierr); 18 ierr = PetscObjectGetComm(obj, &comm);CHKERRQ(ierr); 19 ierr = PetscMalloc2(nrranks, &rsize, nrranks, &rarr);CHKERRQ(ierr); 20 ierr = PetscMalloc2(nrranks, &rreq, nsranks, &sreq);CHKERRQ(ierr); 21 /* exchange array size */ 22 ierr = PetscObjectGetNewTag(obj,&tag);CHKERRQ(ierr); 23 for (r=0; r<nrranks; r++) { 24 ierr = MPI_Irecv(&rsize[r], 1, MPIU_INT, rranks[r], tag, comm, &rreq[r]);CHKERRMPI(ierr); 25 } 26 for (r=0; r<nsranks; r++) { 27 ierr = MPI_Isend(&ssize[r], 1, MPIU_INT, sranks[r], tag, comm, &sreq[r]);CHKERRMPI(ierr); 28 } 29 ierr = MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE);CHKERRMPI(ierr); 30 ierr = MPI_Waitall(nsranks, sreq, MPI_STATUSES_IGNORE);CHKERRMPI(ierr); 31 /* exchange array */ 32 ierr = PetscObjectGetNewTag(obj,&tag);CHKERRQ(ierr); 33 for (r=0; r<nrranks; r++) { 34 ierr = PetscMalloc(rsize[r]*unitsize, &rarr[r]);CHKERRQ(ierr); 35 ierr = MPI_Irecv(rarr[r], rsize[r], dt, rranks[r], tag, comm, &rreq[r]);CHKERRMPI(ierr); 36 } 37 for (r=0; r<nsranks; r++) { 38 ierr = MPI_Isend(sarr[r], ssize[r], dt, sranks[r], tag, comm, &sreq[r]);CHKERRMPI(ierr); 39 } 40 ierr = MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE);CHKERRMPI(ierr); 41 ierr = MPI_Waitall(nsranks, sreq, MPI_STATUSES_IGNORE);CHKERRMPI(ierr); 42 ierr = PetscFree2(rreq, sreq);CHKERRQ(ierr); 43 *rsize_out = rsize; 44 *rarr_out = rarr; 45 PetscFunctionReturn(0); 46 } 47 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; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 ierr = PetscMalloc4(nsranks, &ssize, nsranks, &sarr, nrranks, &rreq, nsranks, &sreq);CHKERRQ(ierr); 62 for (r=0; r<nsranks; r++) { 63 ierr = VecGetLocalSize(svecs[r], &ssize[r]);CHKERRQ(ierr); 64 ierr = VecGetArrayRead(svecs[r], &sarr[r]);CHKERRQ(ierr); 65 } 66 ierr = ExchangeArrayByRank_Private(obj, MPIU_SCALAR, nsranks, sranks, ssize, (const void**)sarr, nrranks, rranks, &rsize, (void***)&rarr);CHKERRQ(ierr); 67 ierr = PetscMalloc1(nrranks, &rvecs_);CHKERRQ(ierr); 68 for (r=0; r<nrranks; r++) { 69 /* set array in two steps to mimic PETSC_OWN_POINTER */ 70 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, 1, rsize[r], NULL, &rvecs_[r]);CHKERRQ(ierr); 71 ierr = VecReplaceArray(rvecs_[r], rarr[r]);CHKERRQ(ierr); 72 } 73 for (r=0; r<nsranks; r++) { 74 ierr = VecRestoreArrayRead(svecs[r], &sarr[r]);CHKERRQ(ierr); 75 } 76 ierr = PetscFree2(rsize, rarr);CHKERRQ(ierr); 77 ierr = PetscFree4(ssize, sarr, rreq, sreq);CHKERRQ(ierr); 78 *rvecs = rvecs_; 79 PetscFunctionReturn(0); 80 } 81 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; 90 91 PetscFunctionBegin; 92 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr); 93 nleaves = roffset[nranks]; 94 ierr = PetscMalloc2(nleaves, rmine1, nleaves, rremote1);CHKERRQ(ierr); 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 ierr = PetscArraycpy(&(*rmine1)[o], &rmine[o], n);CHKERRQ(ierr); 101 ierr = PetscArraycpy(&(*rremote1)[o], &rremote[o], n);CHKERRQ(ierr); 102 ierr = PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);CHKERRQ(ierr); 103 } 104 PetscFunctionReturn(0); 105 } 106 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; 115 116 PetscFunctionBegin; 117 ierr = DMGetCoordinatesLocalSetUp(dm);CHKERRQ(ierr); 118 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);CHKERRQ(ierr); 119 ierr = PetscMalloc1(nranks, coordinatesPerRank);CHKERRQ(ierr); 120 for (r=0; r<nranks; r++) { 121 o = roffset[r]; 122 n = roffset[r+1] - o; 123 ierr = ISCreateGeneral(PETSC_COMM_SELF, n, &rmine[o], PETSC_USE_POINTER, &pointsPerRank);CHKERRQ(ierr); 124 ierr = DMPlexGetConeRecursiveVertices(dm, pointsPerRank, &conesPerRank);CHKERRQ(ierr); 125 ierr = DMGetCoordinatesLocalTuple(dm, conesPerRank, NULL, &(*coordinatesPerRank)[r]);CHKERRQ(ierr); 126 ierr = ISDestroy(&pointsPerRank);CHKERRQ(ierr); 127 ierr = ISDestroy(&conesPerRank);CHKERRQ(ierr); 128 } 129 PetscFunctionReturn(0); 130 } 131 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; 139 140 PetscFunctionBegin; 141 ierr = PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL);CHKERRQ(ierr); 142 ierr = PetscSFGetRootRanks(imsf, &niranks, NULL, &iroffset, &irmine, NULL);CHKERRQ(ierr); 143 if (PetscUnlikely(nileaves != iroffset[niranks])) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nileaves != iroffset[niranks])"); 144 ierr = PetscSFComputeDegreeBegin(sf, °ree);CHKERRQ(ierr); 145 ierr = PetscSFComputeDegreeEnd(sf, °ree);CHKERRQ(ierr); 146 ierr = PetscSFComputeMultiRootOriginalNumbering(sf, degree, NULL, &mRootsOrigNumbering);CHKERRQ(ierr); 147 ierr = PetscMalloc1(nileaves, irmine1);CHKERRQ(ierr); 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 ierr = PetscFree(mRootsOrigNumbering);CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 /*@ 158 DMPlexCheckInterfaceCones - Check that points on inter-partition interfaces have conforming order of cone points. 159 160 Input Parameters: 161 . dm - The DMPlex object 162 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. 166 167 This is mainly intended for debugging/testing purposes. Does not check cone orientation, for this purpose use DMPlexCheckFaces(). 168 169 For the complete list of DMPlexCheck* functions, see DMSetFromOptions(). 170 171 Developer Note: 172 Interface cones are expanded into vertices and then their coordinates are compared. 173 174 Level: developer 175 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; 200 201 PetscFunctionBegin; 202 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 203 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 204 ierr = MPI_Comm_rank(comm, &myrank);CHKERRMPI(ierr); 205 ierr = MPI_Comm_size(comm, &commsize);CHKERRMPI(ierr); 206 if (commsize < 2) PetscFunctionReturn(0); 207 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 208 if (!sf) PetscFunctionReturn(0); 209 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, &mine, &remote);CHKERRQ(ierr); 210 if (nroots < 0) PetscFunctionReturn(0); 211 if (!dm->coordinates && !dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM coordinates must be set"); 212 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 213 ierr = PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr); 214 215 /* Expand sent cones per rank */ 216 ierr = SortByRemote_Private(sf, &rmine1, &rremote1);CHKERRQ(ierr); 217 ierr = GetRecursiveConeCoordinatesPerRank_Private(dm, sf, rmine1, &sntCoordinatesPerRank);CHKERRQ(ierr); 218 219 /* Create inverse SF */ 220 ierr = PetscSFGetMultiSF(sf,&msf);CHKERRQ(ierr); 221 ierr = PetscSFCreateInverseSF(msf,&imsf);CHKERRQ(ierr); 222 ierr = PetscSFSetUp(imsf);CHKERRQ(ierr); 223 ierr = PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL);CHKERRQ(ierr); 224 ierr = PetscSFGetRootRanks(imsf, &niranks, &iranks, &iroffset, &irmine, &irremote);CHKERRQ(ierr); 225 226 /* Compute original numbering of multi-roots (referenced points) */ 227 ierr = PetscSFComputeMultiRootOriginalNumberingByRank_Private(sf, imsf, &mine_orig_numbering);CHKERRQ(ierr); 228 229 /* Expand coordinates of the referred cones per rank */ 230 ierr = GetRecursiveConeCoordinatesPerRank_Private(dm, imsf, mine_orig_numbering, &refCoordinatesPerRank);CHKERRQ(ierr); 231 232 /* Send the coordinates */ 233 ierr = ExchangeVecByRank_Private((PetscObject)sf, nranks, ranks, sntCoordinatesPerRank, niranks, iranks, &recCoordinatesPerRank);CHKERRQ(ierr); 234 235 /* verbose output */ 236 ierr = PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_check_cones_conform_on_interfaces_verbose", &verbose, NULL);CHKERRQ(ierr); 237 if (verbose) { 238 PetscViewer sv, v = PETSC_VIEWER_STDOUT_WORLD; 239 ierr = PetscViewerASCIIPrintf(v, "============\nDMPlexCheckInterfaceCones output\n============\n");CHKERRQ(ierr); 240 ierr = PetscViewerASCIIPushSynchronized(v);CHKERRQ(ierr); 241 ierr = PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", myrank);CHKERRQ(ierr); 242 for (r=0; r<nranks; r++) { 243 ierr = PetscViewerASCIISynchronizedPrintf(v, " r=%D ranks[r]=%d sntCoordinatesPerRank[r]:\n", r, ranks[r]);CHKERRQ(ierr); 244 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 245 ierr = PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);CHKERRQ(ierr); 246 ierr = VecView(sntCoordinatesPerRank[r], sv);CHKERRQ(ierr); 247 ierr = PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);CHKERRQ(ierr); 248 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 249 } 250 ierr = PetscViewerASCIISynchronizedPrintf(v, " ----------\n");CHKERRQ(ierr); 251 for (r=0; r<niranks; r++) { 252 ierr = PetscViewerASCIISynchronizedPrintf(v, " r=%D iranks[r]=%d refCoordinatesPerRank[r]:\n", r, iranks[r]);CHKERRQ(ierr); 253 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 254 ierr = PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);CHKERRQ(ierr); 255 ierr = VecView(refCoordinatesPerRank[r], sv);CHKERRQ(ierr); 256 ierr = PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);CHKERRQ(ierr); 257 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 258 } 259 ierr = PetscViewerASCIISynchronizedPrintf(v, " ----------\n");CHKERRQ(ierr); 260 for (r=0; r<niranks; r++) { 261 ierr = PetscViewerASCIISynchronizedPrintf(v, " r=%D iranks[r]=%d recCoordinatesPerRank[r]:\n", r, iranks[r]);CHKERRQ(ierr); 262 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr); 263 ierr = PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);CHKERRQ(ierr); 264 ierr = VecView(recCoordinatesPerRank[r], sv);CHKERRQ(ierr); 265 ierr = PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);CHKERRQ(ierr); 266 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr); 267 } 268 ierr = PetscViewerFlush(v);CHKERRQ(ierr); 269 ierr = PetscViewerASCIIPopSynchronized(v);CHKERRQ(ierr); 270 } 271 272 /* Compare recCoordinatesPerRank with refCoordinatesPerRank */ 273 for (r=0; r<niranks; r++) { 274 ierr = VecEqual(refCoordinatesPerRank[r], recCoordinatesPerRank[r], &same);CHKERRQ(ierr); 275 if (!same) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interface cones do not conform for remote rank %d", iranks[r]); 276 } 277 278 /* destroy sent stuff */ 279 for (r=0; r<nranks; r++) { 280 ierr = VecDestroy(&sntCoordinatesPerRank[r]);CHKERRQ(ierr); 281 } 282 ierr = PetscFree(sntCoordinatesPerRank);CHKERRQ(ierr); 283 ierr = PetscFree2(rmine1, rremote1);CHKERRQ(ierr); 284 ierr = PetscSFDestroy(&imsf);CHKERRQ(ierr); 285 286 /* destroy referenced stuff */ 287 for (r=0; r<niranks; r++) { 288 ierr = VecDestroy(&refCoordinatesPerRank[r]);CHKERRQ(ierr); 289 } 290 ierr = PetscFree(refCoordinatesPerRank);CHKERRQ(ierr); 291 ierr = PetscFree(mine_orig_numbering);CHKERRQ(ierr); 292 293 /* destroy received stuff */ 294 for (r=0; r<niranks; r++) { 295 ierr = VecDestroy(&recCoordinatesPerRank[r]);CHKERRQ(ierr); 296 } 297 ierr = PetscFree(recCoordinatesPerRank);CHKERRQ(ierr); 298 PetscFunctionReturn(0); 299 } 300