Actual source code: plexcheckinterface.c
petsc-3.12.5 2020-03-29
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 defined(PETSC_USE_DEBUG)
144: if (PetscUnlikely(nileaves != iroffset[niranks])) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nileaves != iroffset[niranks])");
145: #endif
146: PetscSFComputeDegreeBegin(sf, °ree);
147: PetscSFComputeDegreeEnd(sf, °ree);
148: PetscSFComputeMultiRootOriginalNumbering(sf, degree, NULL, &mRootsOrigNumbering);
149: PetscMalloc1(nileaves, irmine1);
150: for (r=0; r<niranks; r++) {
151: o = iroffset[r];
152: n = iroffset[r+1] - o;
153: for (i=0; i<n; i++) (*irmine1)[o+i] = mRootsOrigNumbering[irmine[o+i]];
154: }
155: PetscFree(mRootsOrigNumbering);
156: return(0);
157: }
159: /*@
160: DMPlexCheckConesConformOnInterfaces - Check that points on inter-partition interfaces have conforming order of cone points.
161: 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),
162: 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.
164: Input Parameters:
165: . dm - The DMPlex object
167: Note: This is mainly intended for debugging/testing purposes. Does not check cone orientation, for this purpose use DMPlexCheckFaces().
169: Developer Note: Interface cones are expanded into vertices and then their coordinates are compared.
171: Level: developer
173: .seealso: DMPlexGetCone(), DMPlexGetConeSize(), DMGetPointSF(), DMGetCoordinates(), DMPlexCheckFaces(), DMPlexCheckPointSF(), DMPlexCheckSymmetry(), DMPlexCheckSkeleton()
174: @*/
175: PetscErrorCode DMPlexCheckConesConformOnInterfaces(DM dm)
176: {
177: PetscSF sf;
178: PetscInt nleaves, nranks, nroots;
179: const PetscInt *mine, *roffset, *rmine, *rremote;
180: const PetscSFNode *remote;
181: const PetscMPIInt *ranks;
182: PetscSF msf, imsf;
183: PetscInt nileaves, niranks;
184: const PetscMPIInt *iranks;
185: const PetscInt *iroffset, *irmine, *irremote;
186: PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */
187: PetscInt *mine_orig_numbering;
188: Vec *sntCoordinatesPerRank;
189: Vec *refCoordinatesPerRank;
190: Vec *recCoordinatesPerRank=0;
191: PetscInt r;
192: PetscMPIInt commsize, myrank;
193: PetscBool same;
194: PetscBool verbose=PETSC_FALSE;
195: MPI_Comm comm;
196: PetscErrorCode ierr;
200: PetscObjectGetComm((PetscObject)dm, &comm);
201: MPI_Comm_rank(comm, &myrank);
202: MPI_Comm_size(comm, &commsize);
203: if (commsize < 2) return(0);
204: DMGetPointSF(dm, &sf);
205: if (!sf) return(0);
206: PetscSFGetGraph(sf, &nroots, &nleaves, &mine, &remote);
207: if (nroots < 0) return(0);
208: if (!dm->coordinates && !dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM coordinates must be set");
209: PetscSFSetUp(sf);
210: PetscSFGetRootRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);
212: /* Expand sent cones per rank */
213: SortByRemote_Private(sf, &rmine1, &rremote1);
214: GetRecursiveConeCoordinatesPerRank_Private(dm, sf, rmine1, &sntCoordinatesPerRank);
216: /* Create inverse SF */
217: PetscSFGetMultiSF(sf,&msf);
218: PetscSFCreateInverseSF(msf,&imsf);
219: PetscSFSetUp(imsf);
220: PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL);
221: PetscSFGetRootRanks(imsf, &niranks, &iranks, &iroffset, &irmine, &irremote);
223: /* Compute original numbering of multi-roots (referenced points) */
224: PetscSFComputeMultiRootOriginalNumberingByRank_Private(sf, imsf, &mine_orig_numbering);
226: /* Expand coordinates of the referred cones per rank */
227: GetRecursiveConeCoordinatesPerRank_Private(dm, imsf, mine_orig_numbering, &refCoordinatesPerRank);
229: /* Send the coordinates */
230: ExchangeVecByRank_Private((PetscObject)sf, nranks, ranks, sntCoordinatesPerRank, niranks, iranks, &recCoordinatesPerRank);
232: /* verbose output */
233: PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_check_cones_conform_on_interfaces_verbose", &verbose, NULL);
234: if (verbose) {
235: PetscViewer sv, v = PETSC_VIEWER_STDOUT_WORLD;
236: PetscViewerASCIIPrintf(v, "============\nDMPlexCheckConesConformOnInterfaces output\n============\n");
237: PetscViewerASCIIPushSynchronized(v);
238: PetscViewerASCIISynchronizedPrintf(v, "[%d] --------\n", myrank);
239: for (r=0; r<nranks; r++) {
240: PetscViewerASCIISynchronizedPrintf(v, " r=%D ranks[r]=%d sntCoordinatesPerRank[r]:\n", r, ranks[r]);
241: PetscViewerASCIIPushTab(v);
242: PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);
243: VecView(sntCoordinatesPerRank[r], sv);
244: PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);
245: PetscViewerASCIIPopTab(v);
246: }
247: PetscViewerASCIISynchronizedPrintf(v, " ----------\n");
248: for (r=0; r<niranks; r++) {
249: PetscViewerASCIISynchronizedPrintf(v, " r=%D iranks[r]=%d refCoordinatesPerRank[r]:\n", r, iranks[r]);
250: PetscViewerASCIIPushTab(v);
251: PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);
252: VecView(refCoordinatesPerRank[r], sv);
253: PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);
254: PetscViewerASCIIPopTab(v);
255: }
256: PetscViewerASCIISynchronizedPrintf(v, " ----------\n");
257: for (r=0; r<niranks; r++) {
258: PetscViewerASCIISynchronizedPrintf(v, " r=%D iranks[r]=%d recCoordinatesPerRank[r]:\n", r, iranks[r]);
259: PetscViewerASCIIPushTab(v);
260: PetscViewerGetSubViewer(v,PETSC_COMM_SELF,&sv);
261: VecView(recCoordinatesPerRank[r], sv);
262: PetscViewerRestoreSubViewer(v,PETSC_COMM_SELF,&sv);
263: PetscViewerASCIIPopTab(v);
264: }
265: }
267: /* Compare recCoordinatesPerRank with refCoordinatesPerRank */
268: for (r=0; r<niranks; r++) {
269: VecEqual(refCoordinatesPerRank[r], recCoordinatesPerRank[r], &same);
270: if (!same) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interface cones do not conform for remote rank %d", iranks[r]);
271: }
273: /* destroy sent stuff */
274: for (r=0; r<nranks; r++) {
275: VecDestroy(&sntCoordinatesPerRank[r]);
276: }
277: PetscFree(sntCoordinatesPerRank);
278: PetscFree2(rmine1, rremote1);
279: PetscSFDestroy(&imsf);
281: /* destroy referenced stuff */
282: for (r=0; r<niranks; r++) {
283: VecDestroy(&refCoordinatesPerRank[r]);
284: }
285: PetscFree(refCoordinatesPerRank);
286: PetscFree(mine_orig_numbering);
288: /* destroy received stuff */
289: for (r=0; r<niranks; r++) {
290: VecDestroy(&recCoordinatesPerRank[r]);
291: }
292: PetscFree(recCoordinatesPerRank);
293: return(0);
294: }