1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2
3 /* TODO PetscArrayExchangeBegin/End */
4 /* TODO blocksize */
5 /* TODO move to API ? */
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[])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 ? */
ExchangeVecByRank_Private(PetscObject obj,PetscMPIInt nsranks,const PetscMPIInt sranks[],Vec svecs[],PetscMPIInt nrranks,const PetscMPIInt rranks[],Vec * rvecs[])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
SortByRemote_Private(PetscSF sf,PetscInt * rmine1[],PetscInt * rremote1[])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
GetRecursiveConeCoordinatesPerRank_Private(DM dm,PetscSF sf,PetscInt rmine[],Vec * coordinatesPerRank[])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
PetscSFComputeMultiRootOriginalNumberingByRank_Private(PetscSF sf,PetscSF imsf,PetscInt * irmine1[])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 @*/
DMPlexCheckInterfaceCones(DM dm)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