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, 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;
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: }
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;
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: }
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;
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: }
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;
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: }
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;
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: }
145: /*@
146: DMPlexCheckInterfaceCones - Check that points on inter-partition interfaces have conforming order of cone points.
148: Input Parameter:
149: . dm - The `DMPLEX` object
151: Level: developer
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.
157: This is mainly intended for debugging/testing purposes. Does not check cone orientation, for this purpose use `DMPlexCheckFaces()`.
159: For the complete list of DMPlexCheck* functions, see `DMSetFromOptions()`.
161: Developer Notes:
162: Interface cones are expanded into vertices and then their coordinates are compared.
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;
190: PetscFunctionBegin;
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));
204: /* Expand sent cones per rank */
205: PetscCall(SortByRemote_Private(sf, &rmine1, &rremote1));
206: PetscCall(GetRecursiveConeCoordinatesPerRank_Private(dm, sf, rmine1, &sntCoordinatesPerRank));
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));
215: /* Compute original numbering of multi-roots (referenced points) */
216: PetscCall(PetscSFComputeMultiRootOriginalNumberingByRank_Private(sf, imsf, &mine_orig_numbering));
218: /* Expand coordinates of the referred cones per rank */
219: PetscCall(GetRecursiveConeCoordinatesPerRank_Private(dm, imsf, mine_orig_numbering, &refCoordinatesPerRank));
221: /* Send the coordinates */
222: PetscCall(ExchangeVecByRank_Private((PetscObject)sf, nranks, ranks, sntCoordinatesPerRank, niranks, iranks, &recCoordinatesPerRank));
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: }
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: }
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));
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));
292: /* destroy received stuff */
293: for (r = 0; r < niranks; r++) PetscCall(VecDestroy(&recCoordinatesPerRank[r]));
294: PetscCall(PetscFree(recCoordinatesPerRank));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }