Actual source code: ex106.c
1: static char help[] = "Test space-filling-curve reorder of a distributed cell list\n\n";
3: #include <petscdmplex.h>
4: #include <petscsf.h>
6: static PetscErrorCode SetGlobalCellTags(MPI_Comm comm, PetscInt numCells, PetscInt tags[])
7: {
8: PetscInt off = 0;
9: PetscMPIInt rank;
11: PetscFunctionBeginUser;
12: PetscCallMPI(MPI_Comm_rank(comm, &rank));
13: PetscCallMPI(MPI_Exscan(&numCells, &off, 1, MPIU_INT, MPI_SUM, comm));
14: if (!rank) off = 0;
15: for (PetscInt c = 0; c < numCells; ++c) tags[c] = off + c;
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
19: // Verify that migrationSF maps the input cells onto the output cells one-to-one. Each input cell
20: // carries a globally unique tag. After migration every tag must appear exactly once.
21: static PetscErrorCode CheckMigrationIsPermutation(MPI_Comm comm, PetscSF migrationSF, PetscInt numCells, PetscInt newNumCells, PetscInt NCells)
22: {
23: PetscInt *tags, *newtags, *hist;
25: PetscFunctionBeginUser;
26: PetscCall(PetscMalloc2(numCells, &tags, newNumCells, &newtags));
27: PetscCall(PetscCalloc1(NCells, &hist));
28: PetscCall(SetGlobalCellTags(comm, numCells, tags));
29: PetscCall(PetscSFBcastBegin(migrationSF, MPIU_INT, tags, newtags, MPI_REPLACE));
30: PetscCall(PetscSFBcastEnd(migrationSF, MPIU_INT, tags, newtags, MPI_REPLACE));
31: for (PetscInt c = 0; c < newNumCells; ++c) {
32: PetscCheck(newtags[c] >= 0 && newtags[c] < NCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Migrated tag %" PetscInt_FMT " out of range [0, %" PetscInt_FMT ")", newtags[c], NCells);
33: ++hist[newtags[c]];
34: }
35: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, hist, NCells, MPIU_INT, MPI_SUM, comm));
36: for (PetscInt c = 0; c < NCells; ++c) PetscCheck(hist[c] == 1, comm, PETSC_ERR_PLIB, "Cell tag %" PetscInt_FMT " appears %" PetscInt_FMT " times, not once", c, hist[c]);
37: PetscCall(PetscFree(hist));
38: PetscCall(PetscFree2(tags, newtags));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: // Independent Morton encoder, so the test does not reuse the implementation it checks. It must
43: // mirror the contract of DMPLEXCURVEMORTON: 21 bits per axis over the global bounding box, with
44: // axis 0 in the highest of each interleaved triple.
45: static uint64_t TestZEncode1(PetscInt t)
46: {
47: uint64_t z = (uint64_t)t & 0x1fffff;
49: z = (z | (z << 32)) & UINT64_C(0x1f00000000ffff);
50: z = (z | (z << 16)) & UINT64_C(0x1f0000ff0000ff);
51: z = (z | (z << 8)) & UINT64_C(0x100f00f00f00f00f);
52: z = (z | (z << 4)) & UINT64_C(0x10c30c30c30c30c3);
53: z = (z | (z << 2)) & UINT64_C(0x1249249249249249);
54: return z;
55: }
57: // Verify that the reordered cells form one ascending run of the curve. If useTags is true on every
58: // process, use tags as the tie-breaker for equal codes, matching the implementation's globally
59: // unique cell number. Empty processes may pass NULL for tags.
60: static PetscErrorCode CheckGloballyCurveSorted(MPI_Comm comm, PetscInt spaceDim, PetscInt n, const PetscReal centroids[], PetscBool useTags, const PetscInt tags[])
61: {
62: PetscReal lo[3], hi[3], span[3];
63: PetscInt64 *codes;
64: PetscInt64 *bounds;
65: PetscInt64 prevmax = PETSC_INT64_MIN;
66: PetscInt *counts;
67: PetscInt *tagbounds = NULL;
68: PetscInt prevtag = PETSC_INT_MIN;
69: const PetscInt maxidx = (1 << 21) - 1;
70: PetscMPIInt size;
72: PetscFunctionBeginUser;
73: PetscAssert(!useTags || !n || tags, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "useTags is set but tags is NULL on a process with %" PetscInt_FMT " cells", n);
74: PetscCallMPI(MPI_Comm_size(comm, &size));
75: for (PetscInt d = 0; d < 3; ++d) {
76: lo[d] = PETSC_MAX_REAL;
77: hi[d] = PETSC_MIN_REAL;
78: }
79: for (PetscInt c = 0; c < n; ++c) {
80: for (PetscInt d = 0; d < spaceDim; ++d) {
81: lo[d] = PetscMin(lo[d], centroids[c * spaceDim + d]);
82: hi[d] = PetscMax(hi[d], centroids[c * spaceDim + d]);
83: }
84: }
85: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, lo, 3, MPIU_REAL, MPIU_MIN, comm));
86: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, hi, 3, MPIU_REAL, MPIU_MAX, comm));
87: for (PetscInt d = 0; d < 3; ++d) {
88: if (lo[d] > hi[d]) {
89: lo[d] = 0.;
90: hi[d] = 0.;
91: }
92: }
93: for (PetscInt d = 0; d < 3; ++d) span[d] = hi[d] > lo[d] ? hi[d] - lo[d] : 1.;
95: PetscCall(PetscMalloc1(PetscMax(1, n), &codes));
96: for (PetscInt c = 0; c < n; ++c) {
97: PetscInt q[3] = {0, 0, 0};
99: for (PetscInt d = 0; d < spaceDim; ++d) {
100: const PetscReal t = (centroids[c * spaceDim + d] - lo[d]) / span[d];
102: q[d] = PetscMax(0, PetscMin(maxidx, (PetscInt)(t * (PetscReal)maxidx)));
103: }
104: codes[c] = (PetscInt64)((TestZEncode1(q[0]) << 2) | (TestZEncode1(q[1]) << 1) | TestZEncode1(q[2]));
105: }
106: for (PetscInt c = 1; c < n; ++c) PetscCheck(codes[c - 1] < codes[c] || (codes[c - 1] == codes[c] && (!useTags || tags[c - 1] < tags[c])), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local curve keys not ascending at %" PetscInt_FMT, c);
108: // Exchange each rank's key range. The counts identify empty ranks without reserving a sentinel
109: // value, because PETSC_INT64_MAX is itself a valid Morton code.
110: PetscCall(PetscMalloc1(2 * size, &bounds));
111: PetscCall(PetscMalloc1(size, &counts));
112: PetscCallMPI(MPI_Allgather(&n, 1, MPIU_INT, counts, 1, MPIU_INT, comm));
113: {
114: PetscInt64 mine[2];
116: mine[0] = n ? codes[0] : 0;
117: mine[1] = n ? codes[n - 1] : 0;
118: PetscCallMPI(MPI_Allgather(mine, 2, MPIU_INT64, bounds, 2, MPIU_INT64, comm));
119: }
120: if (useTags) {
121: PetscInt mine[2];
123: PetscCall(PetscMalloc1(2 * size, &tagbounds));
124: mine[0] = n ? tags[0] : 0;
125: mine[1] = n ? tags[n - 1] : 0;
126: PetscCallMPI(MPI_Allgather(mine, 2, MPIU_INT, tagbounds, 2, MPIU_INT, comm));
127: }
128: for (PetscMPIInt r = 0; r < size; ++r) {
129: PetscBool ordered;
131: if (!counts[r]) continue;
132: ordered = prevmax < bounds[2 * r] || (prevmax == bounds[2 * r] && (!useTags || prevtag < tagbounds[2 * r]));
133: PetscCheck(ordered, comm, PETSC_ERR_PLIB, "Curve order broken across ranks before rank %d", r);
134: prevmax = bounds[2 * r + 1];
135: if (useTags) prevtag = tagbounds[2 * r + 1];
136: }
137: PetscCall(PetscFree(tagbounds));
138: PetscCall(PetscFree(counts));
139: PetscCall(PetscFree(bounds));
140: PetscCall(PetscFree(codes));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: // Part A: reorder a lattice of synthetic 3-D centroids. The initial distribution is strided, so
145: // every rank starts with cells spread over the whole domain.
146: static PetscErrorCode TestFromCentroids(MPI_Comm comm, PetscInt N, PetscBool allOnRank0)
147: {
148: PetscSF sf;
149: PetscReal *centroids, *newcentroids;
150: PetscInt *tags, *newtags;
151: PetscInt numCells = 0, newNumCells, NCells = N * N * N, gnew;
152: PetscMPIInt size, rank;
154: PetscFunctionBeginUser;
155: PetscCallMPI(MPI_Comm_size(comm, &size));
156: PetscCallMPI(MPI_Comm_rank(comm, &rank));
157: // Two input distributions. The strided one gives every rank cells spread over the whole domain.
158: // The other puts every cell on rank 0, which is what reading a mesh serially and building in
159: // parallel produces, and which leaves every other rank with nothing to sample.
160: for (PetscInt g = 0; g < NCells; ++g)
161: if (allOnRank0 ? rank == 0 : g % size == rank) ++numCells;
162: PetscCall(PetscMalloc2(PetscMax(1, numCells) * 3, ¢roids, PetscMax(1, numCells), &tags));
163: {
164: PetscInt c = 0;
166: for (PetscInt g = 0; g < NCells; ++g) {
167: if (allOnRank0 ? rank != 0 : g % size != rank) continue;
168: centroids[c * 3 + 0] = (PetscReal)(g % N) + 0.5;
169: centroids[c * 3 + 1] = (PetscReal)((g / N) % N) + 0.5;
170: centroids[c * 3 + 2] = (PetscReal)(g / (N * N)) + 0.5;
171: ++c;
172: }
173: }
174: PetscCall(SetGlobalCellTags(comm, numCells, tags));
175: PetscCall(DMPlexReorderCellListByCurveFromCentroids(comm, DMPLEXCURVEMORTON, 3, numCells, centroids, &sf, &newNumCells));
176: gnew = newNumCells;
177: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gnew, 1, MPIU_INT, MPI_SUM, comm));
178: PetscCheck(gnew == NCells, comm, PETSC_ERR_PLIB, "Global cell count changed from %" PetscInt_FMT " to %" PetscInt_FMT, NCells, gnew);
179: PetscCall(CheckMigrationIsPermutation(comm, sf, numCells, newNumCells, NCells));
181: // Migrate the centroids themselves so the locality of the new distribution can be measured.
182: PetscCall(PetscMalloc2(newNumCells * 3, &newcentroids, newNumCells, &newtags));
183: {
184: MPI_Datatype ctype;
186: PetscCallMPI(MPI_Type_contiguous(3, MPIU_REAL, &ctype));
187: PetscCallMPI(MPI_Type_commit(&ctype));
188: PetscCall(PetscSFBcastBegin(sf, ctype, centroids, newcentroids, MPI_REPLACE));
189: PetscCall(PetscSFBcastEnd(sf, ctype, centroids, newcentroids, MPI_REPLACE));
190: PetscCallMPI(MPI_Type_free(&ctype));
191: }
192: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, tags, newtags, MPI_REPLACE));
193: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, tags, newtags, MPI_REPLACE));
194: PetscCall(CheckGloballyCurveSorted(comm, 3, newNumCells, newcentroids, PETSC_TRUE, newtags));
195: // The reorder must redistribute, not merely permute in place, and it must split the cells evenly.
196: // The second pass equidistributes from the global position of each cell, so the counts differ by
197: // at most one. The splitters alone only bound the busiest rank at twice the average.
198: {
199: PetscInt hi = newNumCells, lo = newNumCells;
201: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hi, 1, MPIU_INT, MPI_MAX, comm));
202: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &lo, 1, MPIU_INT, MPI_MIN, comm));
203: PetscCheck(hi - lo <= 1, comm, PETSC_ERR_PLIB, "Cells per rank run from %" PetscInt_FMT " to %" PetscInt_FMT ", which is not an exact split", lo, hi);
204: // Every rank must receive cells once there are enough to go round.
205: PetscCheck(NCells < (PetscInt)size || lo > 0, comm, PETSC_ERR_PLIB, "A rank received no cells");
206: }
207: PetscCall(PetscPrintf(comm, "FromCentroids: N=%" PetscInt_FMT " cells=%" PetscInt_FMT " permutation ok, globally curve sorted, balanced\n", N, NCells));
208: PetscCall(PetscSFDestroy(&sf));
209: PetscCall(PetscFree2(newcentroids, newtags));
210: PetscCall(PetscFree2(centroids, tags));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: // Part B: reorder the connectivity of an N x N quadrilateral grid, then build a DMPlex from the
215: // reordered list. The vertex distribution stays in natural order, as the reorder requires.
216: static PetscErrorCode TestCellList(MPI_Comm comm, PetscInt N)
217: {
218: DM dm;
219: PetscSF sf;
220: PetscInt *cells, *newcells, *cellsSaved = NULL;
221: PetscReal *coords;
222: PetscLayout vlayout;
223: PetscInt numCells = 0, newNumCells, NCells = N * N, NVertices = (N + 1) * (N + 1);
224: PetscInt numVertices, vStart, vEnd, cStart, cEnd, gcells;
225: PetscMPIInt size, rank;
227: PetscFunctionBeginUser;
228: PetscCallMPI(MPI_Comm_size(comm, &size));
229: PetscCallMPI(MPI_Comm_rank(comm, &rank));
230: for (PetscInt g = 0; g < NCells; ++g)
231: if (g % size == rank) ++numCells;
232: PetscCall(PetscMalloc1(numCells * 4, &cells));
233: {
234: PetscInt c = 0;
236: for (PetscInt g = 0; g < NCells; ++g) {
237: const PetscInt i = g % N, j = g / N;
239: if (g % size != rank) continue;
240: cells[c * 4 + 0] = j * (N + 1) + i;
241: cells[c * 4 + 1] = j * (N + 1) + i + 1;
242: cells[c * 4 + 2] = (j + 1) * (N + 1) + i + 1;
243: cells[c * 4 + 3] = (j + 1) * (N + 1) + i;
244: ++c;
245: }
246: }
247: // Own a contiguous slice of the vertices, matching what DMPlexCreateFromCellListParallelPetsc() expects.
248: PetscCall(PetscLayoutCreate(comm, &vlayout));
249: PetscCall(PetscLayoutSetSize(vlayout, NVertices));
250: PetscCall(PetscLayoutSetBlockSize(vlayout, 1));
251: PetscCall(PetscLayoutSetUp(vlayout));
252: PetscCall(PetscLayoutGetRange(vlayout, &vStart, &vEnd));
253: PetscCall(PetscLayoutDestroy(&vlayout));
254: numVertices = vEnd - vStart;
255: PetscCall(PetscMalloc1(numVertices * 2, &coords));
256: for (PetscInt v = vStart; v < vEnd; ++v) {
257: coords[(v - vStart) * 2 + 0] = (PetscReal)(v % (N + 1));
258: coords[(v - vStart) * 2 + 1] = (PetscReal)(v / (N + 1));
259: }
261: PetscCall(DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, numCells, 4, cells, 2, numVertices, NVertices, coords, &sf, &newNumCells, &newcells));
262: PetscCall(CheckMigrationIsPermutation(comm, sf, numCells, newNumCells, NCells));
263: for (PetscInt i = 0; i < newNumCells * 4; ++i) PetscCheck(newcells[i] >= 0 && newcells[i] < NVertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reordered connectivity entry %" PetscInt_FMT " out of range", newcells[i]);
265: // The reordered list must build a valid, interpolated, distributed plex.
266: PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, 2, newNumCells, numVertices, NVertices, 4, PETSC_TRUE, newcells, 2, coords, NULL, &cellsSaved, &dm));
267: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
268: gcells = cEnd - cStart;
269: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gcells, 1, MPIU_INT, MPI_SUM, comm));
270: PetscCheck(gcells == NCells, comm, PETSC_ERR_PLIB, "Built plex has %" PetscInt_FMT " cells, expected %" PetscInt_FMT, gcells, NCells);
271: PetscCall(DMPlexCheck(dm));
272: PetscCall(PetscPrintf(comm, "CellList: N=%" PetscInt_FMT " cells=%" PetscInt_FMT " plex built and checked\n", N, NCells));
273: PetscCall(PetscFree(cellsSaved));
274: PetscCall(DMDestroy(&dm));
275: PetscCall(PetscSFDestroy(&sf));
276: PetscCall(PetscFree(newcells));
277: PetscCall(PetscFree(coords));
278: PetscCall(PetscFree(cells));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: // Part C: the reorder exists to make parallel interpolation cheap. After interpolation the number
283: // of shared points measures how much of the mesh sits on rank boundaries. Reordering must not
284: // increase it.
285: static PetscErrorCode BuildAndCountSharedPoints(MPI_Comm comm, PetscInt numCells, const PetscInt cells[], PetscInt numVertices, PetscInt NVertices, const PetscReal coords[], PetscInt *nshared)
286: {
287: DM dm;
288: PetscSF pointSF;
289: PetscInt *cellsSaved = NULL;
290: PetscInt nleaves;
292: PetscFunctionBeginUser;
293: PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, 2, numCells, numVertices, NVertices, 4, PETSC_TRUE, cells, 2, coords, NULL, &cellsSaved, &dm));
294: PetscCall(DMGetPointSF(dm, &pointSF));
295: PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, NULL, NULL));
296: *nshared = nleaves;
297: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, nshared, 1, MPIU_INT, MPI_SUM, comm));
298: PetscCall(PetscFree(cellsSaved));
299: PetscCall(DMDestroy(&dm));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode TestLocalityImproves(MPI_Comm comm, PetscInt N)
304: {
305: PetscSF sf;
306: PetscInt *cells, *newcells;
307: PetscReal *coords;
308: PetscLayout vlayout;
309: PetscInt numCells = 0, newNumCells, NCells = N * N, NVertices = (N + 1) * (N + 1);
310: PetscInt numVertices, vStart, vEnd, sharedBefore, sharedAfter;
311: PetscMPIInt size, rank;
313: PetscFunctionBeginUser;
314: PetscCallMPI(MPI_Comm_size(comm, &size));
315: PetscCallMPI(MPI_Comm_rank(comm, &rank));
316: for (PetscInt g = 0; g < NCells; ++g)
317: if (g % size == rank) ++numCells;
318: PetscCall(PetscMalloc1(numCells * 4, &cells));
319: {
320: PetscInt c = 0;
322: for (PetscInt g = 0; g < NCells; ++g) {
323: const PetscInt i = g % N, j = g / N;
325: if (g % size != rank) continue;
326: cells[c * 4 + 0] = j * (N + 1) + i;
327: cells[c * 4 + 1] = j * (N + 1) + i + 1;
328: cells[c * 4 + 2] = (j + 1) * (N + 1) + i + 1;
329: cells[c * 4 + 3] = (j + 1) * (N + 1) + i;
330: ++c;
331: }
332: }
333: PetscCall(PetscLayoutCreate(comm, &vlayout));
334: PetscCall(PetscLayoutSetSize(vlayout, NVertices));
335: PetscCall(PetscLayoutSetBlockSize(vlayout, 1));
336: PetscCall(PetscLayoutSetUp(vlayout));
337: PetscCall(PetscLayoutGetRange(vlayout, &vStart, &vEnd));
338: PetscCall(PetscLayoutDestroy(&vlayout));
339: numVertices = vEnd - vStart;
340: PetscCall(PetscMalloc1(numVertices * 2, &coords));
341: for (PetscInt v = vStart; v < vEnd; ++v) {
342: coords[(v - vStart) * 2 + 0] = (PetscReal)(v % (N + 1));
343: coords[(v - vStart) * 2 + 1] = (PetscReal)(v / (N + 1));
344: }
346: PetscCall(BuildAndCountSharedPoints(comm, numCells, cells, numVertices, NVertices, coords, &sharedBefore));
347: PetscCall(DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, numCells, 4, cells, 2, numVertices, NVertices, coords, &sf, &newNumCells, &newcells));
348: PetscCall(BuildAndCountSharedPoints(comm, newNumCells, newcells, numVertices, NVertices, coords, &sharedAfter));
349: PetscCheck(sharedAfter <= sharedBefore, comm, PETSC_ERR_PLIB, "Shared points grew from %" PetscInt_FMT " to %" PetscInt_FMT, sharedBefore, sharedAfter);
350: // A strided input distribution shares almost every point. The reorder must cut that sharply, not
351: // merely avoid making it worse, so require at least a factor of two on more than one process.
352: // How much there is to gain depends on how many cells each process receives. With fewer than a
353: // couple of cells per process every cell touches a boundary and nothing can improve, so only
354: // require that the reorder does no harm. Above that, require a modest reduction, and once each
355: // process holds a real block require a factor of two. Measured values: a 2x2 grid over 8
356: // processes gives no change, an 8x8 grid over 8 processes gives 1.68, and a 16x16 grid gives
357: // more than 4.
358: if (size > 1 && NCells >= 8 * (PetscInt)size) {
359: PetscCheck(4 * sharedAfter <= 3 * sharedBefore, comm, PETSC_ERR_PLIB, "Shared points only fell from %" PetscInt_FMT " to %" PetscInt_FMT ", less than a quarter", sharedBefore, sharedAfter);
360: PetscCheck(NCells < 32 * (PetscInt)size || 2 * sharedAfter <= sharedBefore, comm, PETSC_ERR_PLIB, "Shared points only fell from %" PetscInt_FMT " to %" PetscInt_FMT ", less than a factor of two", sharedBefore, sharedAfter);
361: }
362: // The counts depend on the number of processes, so keep them out of the reference output.
363: PetscCall(PetscPrintf(comm, "Locality: reorder cuts shared points\n"));
364: PetscCall(PetscInfo(NULL, "shared points %" PetscInt_FMT " -> %" PetscInt_FMT "\n", sharedBefore, sharedAfter));
365: PetscCall(PetscSFDestroy(&sf));
366: PetscCall(PetscFree(newcells));
367: PetscCall(PetscFree(coords));
368: PetscCall(PetscFree(cells));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: // Part D: degenerate geometry. A curve code alone cannot order cells that quantize to the same grid
373: // point, and one distant node is enough to make the whole bulk of a mesh do that. The reorder must
374: // still balance, because a rank that receives every cell is worse than no reorder at all.
375: static PetscErrorCode TestDegenerateGeometry(MPI_Comm comm, PetscInt N)
376: {
377: const char *names[] = {"identical centroids", "outlier stretches box", "mostly coincident", "corner on every rank"};
378: const PetscInt nkinds = (PetscInt)(sizeof(names) / sizeof(names[0]));
379: PetscMPIInt size, rank;
381: PetscFunctionBeginUser;
382: PetscCallMPI(MPI_Comm_size(comm, &size));
383: PetscCallMPI(MPI_Comm_rank(comm, &rank));
384: for (PetscInt kind = 0; kind < nkinds; ++kind) {
385: PetscSF sf;
386: PetscReal *cent, *newcent;
387: PetscInt *tags, *newtags;
388: PetscInt NCells = kind == 3 ? 4095 : N * N * N, numCells = 0, newNumCells, lo, hi, tot;
390: for (PetscInt g = 0; g < NCells; ++g)
391: if (g % size == rank) ++numCells;
392: PetscCall(PetscMalloc2(PetscMax(1, numCells) * 3, ¢, PetscMax(1, numCells), &tags));
393: {
394: PetscInt c = 0;
396: for (PetscInt g = 0; g < NCells; ++g) {
397: if (g % size != rank) continue;
398: if (kind == 0) { // every centroid the same point
399: cent[c * 3 + 0] = 1.5;
400: cent[c * 3 + 1] = 2.5;
401: cent[c * 3 + 2] = 3.5;
402: } else if (kind == 1) { // a tight mesh plus one very distant node
403: cent[c * 3 + 0] = (PetscReal)(g % N) * 1e-3;
404: cent[c * 3 + 1] = (PetscReal)((g / N) % N) * 1e-3;
405: cent[c * 3 + 2] = (PetscReal)(g / (N * N)) * 1e-3;
406: if (g == 0) {
407: cent[c * 3 + 0] = 1e6;
408: cent[c * 3 + 1] = 1e6;
409: cent[c * 3 + 2] = 1e6;
410: }
411: } else if (kind == 3) {
412: // Every rank holds cells across the whole box, and its first cell sits on the corner, so
413: // every rank's own smallest curve code equals the global smallest. A mesh file in
414: // generator order looks like this. Too few samples then place a splitter on that shared
415: // smallest value, and nearly every cell lands on one rank. The cell count is deliberately
416: // indivisible so the per-rank sample counts differ after rounding.
417: if (c == 0) {
418: cent[c * 3 + 0] = 0.;
419: cent[c * 3 + 1] = 0.;
420: cent[c * 3 + 2] = 0.;
421: } else {
422: cent[c * 3 + 0] = (PetscReal)((g * 7919) % 64);
423: cent[c * 3 + 1] = (PetscReal)((g * 6271) % 64);
424: cent[c * 3 + 2] = (PetscReal)((g * 4643) % 64);
425: }
426: } else { // three quarters of the cells on one point
427: if (g % 4) {
428: cent[c * 3 + 0] = 5.;
429: cent[c * 3 + 1] = 5.;
430: cent[c * 3 + 2] = 5.;
431: } else {
432: cent[c * 3 + 0] = (PetscReal)(g % N);
433: cent[c * 3 + 1] = (PetscReal)((g / N) % N);
434: cent[c * 3 + 2] = (PetscReal)(g / (N * N));
435: }
436: }
437: ++c;
438: }
439: }
440: PetscCall(SetGlobalCellTags(comm, numCells, tags));
441: PetscCall(DMPlexReorderCellListByCurveFromCentroids(comm, DMPLEXCURVEMORTON, 3, numCells, cent, &sf, &newNumCells));
442: PetscCall(CheckMigrationIsPermutation(comm, sf, numCells, newNumCells, NCells));
443: PetscCall(PetscMalloc2(newNumCells * 3, &newcent, newNumCells, &newtags));
444: {
445: MPI_Datatype ctype;
447: PetscCallMPI(MPI_Type_contiguous(3, MPIU_REAL, &ctype));
448: PetscCallMPI(MPI_Type_commit(&ctype));
449: PetscCall(PetscSFBcastBegin(sf, ctype, cent, newcent, MPI_REPLACE));
450: PetscCall(PetscSFBcastEnd(sf, ctype, cent, newcent, MPI_REPLACE));
451: PetscCallMPI(MPI_Type_free(&ctype));
452: }
453: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, tags, newtags, MPI_REPLACE));
454: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, tags, newtags, MPI_REPLACE));
455: PetscCall(CheckGloballyCurveSorted(comm, 3, newNumCells, newcent, PETSC_TRUE, newtags));
456: lo = newNumCells;
457: hi = newNumCells;
458: tot = newNumCells;
459: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &lo, 1, MPIU_INT, MPI_MIN, comm));
460: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hi, 1, MPIU_INT, MPI_MAX, comm));
461: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &tot, 1, MPIU_INT, MPI_SUM, comm));
462: PetscCheck(tot == NCells, comm, PETSC_ERR_PLIB, "%s: cell count changed from %" PetscInt_FMT " to %" PetscInt_FMT, names[kind], NCells, tot);
463: // The split must be exact whatever the geometry does to the curve codes. The permutation and
464: // key-order checks above use this same degenerate input, because final balance alone cannot see
465: // a first-pass defect after the exact split.
466: PetscCheck(hi - lo <= 1, comm, PETSC_ERR_PLIB, "%s: cells per rank run from %" PetscInt_FMT " to %" PetscInt_FMT ", which is not an exact split", names[kind], lo, hi);
467: // Every rank must receive cells. The curve codes here are all equal or nearly so, so a split
468: // taken from the codes alone would leave one rank with every cell and the rest with none.
469: PetscCheck(NCells < (PetscInt)size || lo > 0, comm, PETSC_ERR_PLIB, "%s: a rank received no cells; the curve codes do not separate these centroids, so the split must come from the cell numbering", names[kind]);
470: PetscCall(PetscSFDestroy(&sf));
471: PetscCall(PetscFree2(newcent, newtags));
472: PetscCall(PetscFree2(cent, tags));
473: }
474: PetscCall(PetscPrintf(comm, "Degenerate geometry: %" PetscInt_FMT " cases permutation ok, globally key sorted, balanced\n", nkinds));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: // Part E: the curve is selected by name. An unregistered name must be rejected, because a silent
479: // fall back to a default curve would order the cells along a curve the caller did not ask for.
480: static PetscErrorCode TestUnknownCurveType(MPI_Comm comm)
481: {
482: const char *bad[] = {"hilbert", "morton ", "", NULL};
483: const PetscInt nbad = (PetscInt)(sizeof(bad) / sizeof(bad[0]));
484: PetscReal cent[3] = {0.5, 0.5, 0.5};
485: PetscSF sf = NULL;
486: PetscInt *newcells = NULL;
487: PetscInt newNumCells = -1;
488: PetscErrorCode ierr;
490: PetscFunctionBeginUser;
491: for (PetscInt i = 0; i < nbad; ++i) {
492: // Every process passes the same name, so the collective check fails on all of them together.
493: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
494: ierr = DMPlexReorderCellListByCurveFromCentroids(comm, bad[i], 3, 1, cent, &sf, &newNumCells);
495: PetscCall(PetscPopErrorHandler());
496: PetscCheck(ierr == PETSC_ERR_ARG_UNKNOWN_TYPE, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurveFromCentroids() returned %d for curve \"%s\", not PETSC_ERR_ARG_UNKNOWN_TYPE", (int)ierr, bad[i] ? bad[i] : "(null)");
497: // The connectivity interface checks the name through the centroid routine, so it gathers the
498: // corner coordinates first. Pass no cells, which keeps that gather empty.
499: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
500: ierr = DMPlexReorderCellListByCurve(comm, bad[i], 0, 4, NULL, 3, 0, 0, NULL, &sf, &newNumCells, &newcells);
501: PetscCall(PetscPopErrorHandler());
502: PetscCheck(ierr == PETSC_ERR_ARG_UNKNOWN_TYPE, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurve() returned %d for curve \"%s\", not PETSC_ERR_ARG_UNKNOWN_TYPE", (int)ierr, bad[i] ? bad[i] : "(null)");
503: }
504: // A rejected call must leave the output arguments alone, so the caller frees nothing.
505: PetscCheck(!sf && !newcells && newNumCells == -1, comm, PETSC_ERR_PLIB, "A rejected call wrote to its output arguments");
506: PetscCall(PetscPrintf(comm, "Unknown curve: %" PetscInt_FMT " names rejected by both interfaces\n", nbad));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: // Part F: one-dimensional coordinates. The Morton encoder interleaves three axes, and Parts A to D
511: // use two or three of them, so a single axis exercises a path of its own. In one dimension the
512: // curve is the axis itself, so the reordered cells must come out in ascending coordinate order.
513: static PetscErrorCode TestOneDimensional(MPI_Comm comm, PetscInt N)
514: {
515: DM dm;
516: PetscSF sf, sfCells;
517: PetscReal *cent, *newcent, *coords;
518: PetscInt *cells, *newcells, *cellsSaved = NULL;
519: PetscLayout vlayout;
520: PetscInt numCells = 0, nnewCent, nnewList, NCells = N * N * N, NVertices = N * N * N + 1;
521: PetscInt numVertices, vStart, vEnd, gnew, cStart, cEnd, gcells;
522: PetscMPIInt size, rank;
524: PetscFunctionBeginUser;
525: PetscCallMPI(MPI_Comm_size(comm, &size));
526: PetscCallMPI(MPI_Comm_rank(comm, &rank));
527: // A line of segments, distributed with a stride so that every process starts with cells spread
528: // over the whole line.
529: for (PetscInt g = 0; g < NCells; ++g)
530: if (g % size == rank) ++numCells;
531: PetscCall(PetscMalloc2(PetscMax(1, numCells), ¢, PetscMax(1, numCells) * 2, &cells));
532: {
533: PetscInt c = 0;
535: for (PetscInt g = 0; g < NCells; ++g) {
536: if (g % size != rank) continue;
537: cent[c] = (PetscReal)g + 0.5;
538: cells[c * 2 + 0] = g;
539: cells[c * 2 + 1] = g + 1;
540: ++c;
541: }
542: }
544: // The centroid interface in one dimension.
545: PetscCall(DMPlexReorderCellListByCurveFromCentroids(comm, DMPLEXCURVEMORTON, 1, numCells, cent, &sf, &nnewCent));
546: gnew = nnewCent;
547: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gnew, 1, MPIU_INT, MPI_SUM, comm));
548: PetscCheck(gnew == NCells, comm, PETSC_ERR_PLIB, "Global cell count changed from %" PetscInt_FMT " to %" PetscInt_FMT, NCells, gnew);
549: PetscCall(CheckMigrationIsPermutation(comm, sf, numCells, nnewCent, NCells));
550: PetscCall(PetscMalloc1(PetscMax(1, nnewCent), &newcent));
551: PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, cent, newcent, MPI_REPLACE));
552: PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, cent, newcent, MPI_REPLACE));
553: PetscCall(CheckGloballyCurveSorted(comm, 1, nnewCent, newcent, PETSC_FALSE, NULL));
554: // On one axis the centroids are distinct and the curve is the axis, so the order is strict.
555: for (PetscInt c = 1; c < nnewCent; ++c) PetscCheck(newcent[c - 1] < newcent[c], PETSC_COMM_SELF, PETSC_ERR_PLIB, "One-dimensional centroids not strictly ascending at %" PetscInt_FMT, c);
557: // The connectivity interface in one dimension. Own a contiguous slice of the vertices.
558: PetscCall(PetscLayoutCreate(comm, &vlayout));
559: PetscCall(PetscLayoutSetSize(vlayout, NVertices));
560: PetscCall(PetscLayoutSetBlockSize(vlayout, 1));
561: PetscCall(PetscLayoutSetUp(vlayout));
562: PetscCall(PetscLayoutGetRange(vlayout, &vStart, &vEnd));
563: PetscCall(PetscLayoutDestroy(&vlayout));
564: numVertices = vEnd - vStart;
565: PetscCall(PetscMalloc1(PetscMax(1, numVertices), &coords));
566: for (PetscInt v = vStart; v < vEnd; ++v) coords[v - vStart] = (PetscReal)v;
567: // Pass PETSC_DECIDE for the global vertex count, which the routine accepts and which Part B does
568: // not exercise. The layout then sums the local counts, which gives the same NVertices.
569: PetscCall(DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, numCells, 2, cells, 1, numVertices, PETSC_DECIDE, coords, &sfCells, &nnewList, &newcells));
570: PetscCheck(nnewList == nnewCent, comm, PETSC_ERR_PLIB, "The two interfaces split the same cells differently, %" PetscInt_FMT " against %" PetscInt_FMT, nnewList, nnewCent);
571: PetscCall(CheckMigrationIsPermutation(comm, sfCells, numCells, nnewList, NCells));
572: // Each segment must arrive whole, and the segments must arrive in ascending order.
573: for (PetscInt c = 0; c < nnewList; ++c) {
574: PetscCheck(newcells[c * 2 + 1] == newcells[c * 2] + 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Segment %" PetscInt_FMT " has vertices %" PetscInt_FMT " and %" PetscInt_FMT ", not consecutive", c, newcells[c * 2], newcells[c * 2 + 1]);
575: PetscCheck(!c || newcells[(c - 1) * 2] < newcells[c * 2], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Segments not ascending at %" PetscInt_FMT, c);
576: }
578: // The reordered list must build a valid one-dimensional plex.
579: PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, 1, nnewList, numVertices, NVertices, 2, PETSC_TRUE, newcells, 1, coords, NULL, &cellsSaved, &dm));
580: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
581: gcells = cEnd - cStart;
582: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gcells, 1, MPIU_INT, MPI_SUM, comm));
583: PetscCheck(gcells == NCells, comm, PETSC_ERR_PLIB, "Built plex has %" PetscInt_FMT " cells, expected %" PetscInt_FMT, gcells, NCells);
584: PetscCall(DMPlexCheck(dm));
585: PetscCall(PetscPrintf(comm, "OneDimensional: cells ascending on one axis, plex built and checked\n"));
586: PetscCall(PetscFree(cellsSaved));
587: PetscCall(DMDestroy(&dm));
588: PetscCall(PetscSFDestroy(&sfCells));
589: PetscCall(PetscSFDestroy(&sf));
590: PetscCall(PetscFree(newcells));
591: PetscCall(PetscFree(coords));
592: PetscCall(PetscFree(newcent));
593: PetscCall(PetscFree2(cent, cells));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: // Part G: argument validation. Every process passes the same invalid value, so each call fails
598: // collectively and the error handler returns the code instead of aborting.
599: static PetscErrorCode TestArgumentValidation(MPI_Comm comm)
600: {
601: const PetscInt baddim[] = {0, 4, -1};
602: const PetscInt ncorner[] = {0, -3};
603: const PetscInt ndim = (PetscInt)(sizeof(baddim) / sizeof(baddim[0]));
604: const PetscInt ncorn = (PetscInt)(sizeof(ncorner) / sizeof(ncorner[0]));
605: PetscReal cent[3] = {0.5, 0.5, 0.5};
606: PetscSF sf = NULL;
607: PetscInt *newcells = NULL;
608: PetscInt newNumCells = -1, nrejected = 0;
609: PetscErrorCode ierr;
611: PetscFunctionBeginUser;
612: // The centroid interface: spaceDim must be in [1, 3] and numCells must not be negative.
613: for (PetscInt i = 0; i < ndim; ++i) {
614: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
615: ierr = DMPlexReorderCellListByCurveFromCentroids(comm, DMPLEXCURVEMORTON, baddim[i], 1, cent, &sf, &newNumCells);
616: PetscCall(PetscPopErrorHandler());
617: PetscCheck(ierr == PETSC_ERR_ARG_OUTOFRANGE, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurveFromCentroids() returned %d for spaceDim %" PetscInt_FMT ", not PETSC_ERR_ARG_OUTOFRANGE", (int)ierr, baddim[i]);
618: ++nrejected;
619: }
620: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
621: ierr = DMPlexReorderCellListByCurveFromCentroids(comm, DMPLEXCURVEMORTON, 3, -1, cent, &sf, &newNumCells);
622: PetscCall(PetscPopErrorHandler());
623: PetscCheck(ierr == PETSC_ERR_ARG_OUTOFRANGE, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurveFromCentroids() returned %d for a negative numCells, not PETSC_ERR_ARG_OUTOFRANGE", (int)ierr);
624: ++nrejected;
626: // The connectivity interface: the same two ranges, and numCorners must be positive. Each check
627: // runs before the routine allocates, so a rejected call leaks nothing.
628: for (PetscInt i = 0; i < ndim; ++i) {
629: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
630: ierr = DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, 0, 4, NULL, baddim[i], 0, 0, NULL, &sf, &newNumCells, &newcells);
631: PetscCall(PetscPopErrorHandler());
632: PetscCheck(ierr == PETSC_ERR_ARG_OUTOFRANGE, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurve() returned %d for spaceDim %" PetscInt_FMT ", not PETSC_ERR_ARG_OUTOFRANGE", (int)ierr, baddim[i]);
633: ++nrejected;
634: }
635: for (PetscInt i = 0; i < ncorn; ++i) {
636: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
637: ierr = DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, 0, ncorner[i], NULL, 3, 0, 0, NULL, &sf, &newNumCells, &newcells);
638: PetscCall(PetscPopErrorHandler());
639: PetscCheck(ierr == PETSC_ERR_ARG_OUTOFRANGE, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurve() returned %d for numCorners %" PetscInt_FMT ", not PETSC_ERR_ARG_OUTOFRANGE", (int)ierr, ncorner[i]);
640: ++nrejected;
641: }
642: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
643: ierr = DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, -1, 4, NULL, 3, 0, 0, NULL, &sf, &newNumCells, &newcells);
644: PetscCall(PetscPopErrorHandler());
645: PetscCheck(ierr == PETSC_ERR_ARG_OUTOFRANGE, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurve() returned %d for a negative numCells, not PETSC_ERR_ARG_OUTOFRANGE", (int)ierr);
646: ++nrejected;
648: // The corner buffer of DMPlexReorderCellListByCurve() must fit in a PetscInt. The check runs
649: // before the routine reads the connectivity, so the array stays NULL and nothing is allocated.
650: // With 64-bit indices no reachable count passes PETSC_INT_MAX, so this case is 32-bit only.
651: if (!PetscDefined(USE_64BIT_INDICES)) {
652: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
653: ierr = DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, 300000000, 4, NULL, 2, 0, 0, NULL, &sf, &newNumCells, &newcells);
654: PetscCall(PetscPopErrorHandler());
655: PetscCheck(ierr == PETSC_ERR_SUP, comm, PETSC_ERR_PLIB, "DMPlexReorderCellListByCurve() returned %d for a corner buffer of 2.4e9 reals, not PETSC_ERR_SUP", (int)ierr);
656: ++nrejected;
657: }
659: PetscCheck(!sf && !newcells && newNumCells == -1, comm, PETSC_ERR_PLIB, "A rejected call wrote to its output arguments");
660: // Keep the count out of the message. It differs between a 32-bit and a 64-bit index build.
661: PetscCall(PetscPrintf(comm, "Validation: every invalid argument set rejected\n"));
662: PetscCall(PetscInfo(NULL, "invalid argument sets rejected %" PetscInt_FMT "\n", nrejected));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: // Part H: no cells anywhere. A reader that finds no cells of the requested type hands the reorder
667: // an empty list on every process. There is then no sample to take, so no splitter can come from the
668: // data, and both interfaces must return an empty migration rather than fail.
669: static PetscErrorCode TestNoCells(MPI_Comm comm)
670: {
671: PetscSF sf;
672: PetscInt *newcells = NULL;
673: PetscInt newNumCells = -1, nroots, nleaves;
675: PetscFunctionBeginUser;
676: PetscCall(DMPlexReorderCellListByCurveFromCentroids(comm, DMPLEXCURVEMORTON, 3, 0, NULL, &sf, &newNumCells));
677: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
678: PetscCheck(newNumCells == 0 && nroots == 0 && nleaves == 0, comm, PETSC_ERR_PLIB, "Empty input gave %" PetscInt_FMT " cells with %" PetscInt_FMT " roots and %" PetscInt_FMT " leaves", newNumCells, nroots, nleaves);
679: PetscCall(PetscSFDestroy(&sf));
681: newNumCells = -1;
682: PetscCall(DMPlexReorderCellListByCurve(comm, DMPLEXCURVEMORTON, 0, 4, NULL, 2, 0, PETSC_DECIDE, NULL, &sf, &newNumCells, &newcells));
683: PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
684: PetscCheck(newNumCells == 0 && nroots == 0 && nleaves == 0, comm, PETSC_ERR_PLIB, "Empty connectivity gave %" PetscInt_FMT " cells with %" PetscInt_FMT " roots and %" PetscInt_FMT " leaves", newNumCells, nroots, nleaves);
685: PetscCall(PetscSFDestroy(&sf));
686: PetscCall(PetscFree(newcells));
687: PetscCall(PetscPrintf(comm, "NoCells: empty input gives an empty migration\n"));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: // Part I: a mesh embedded in more than three dimensions. The curve interleaves three axes and the
692: // bounding box holds three values, so DMPlexGetOrdering() must reject a higher-dimensional mesh
693: // rather than write past that box. The check runs before the work arrays exist, so the rejected call
694: // frees everything, which every test proves because the suite runs with -malloc_dump.
695: static PetscErrorCode TestHighCoordinateDim(MPI_Comm comm)
696: {
697: DM dm;
698: IS perm = NULL;
699: PetscLayout vlayout;
700: PetscReal *coords;
701: PetscInt *cells, *cellsSaved = NULL;
702: PetscInt numCells = 0, numVertices, vStart, vEnd;
703: const PetscInt N = 2, NCells = 4, NVertices = 9;
704: PetscErrorCode ierr;
705: PetscMPIInt size, rank;
707: PetscFunctionBeginUser;
708: PetscCallMPI(MPI_Comm_size(comm, &size));
709: PetscCallMPI(MPI_Comm_rank(comm, &rank));
710: // A 2 x 2 quadrilateral grid whose vertices carry four coordinates each.
711: for (PetscInt g = 0; g < NCells; ++g)
712: if (g % size == rank) ++numCells;
713: PetscCall(PetscMalloc1(PetscMax(1, numCells) * 4, &cells));
714: {
715: PetscInt c = 0;
717: for (PetscInt g = 0; g < NCells; ++g) {
718: const PetscInt i = g % N, j = g / N;
720: if (g % size != rank) continue;
721: cells[c * 4 + 0] = j * (N + 1) + i;
722: cells[c * 4 + 1] = j * (N + 1) + i + 1;
723: cells[c * 4 + 2] = (j + 1) * (N + 1) + i + 1;
724: cells[c * 4 + 3] = (j + 1) * (N + 1) + i;
725: ++c;
726: }
727: }
728: PetscCall(PetscLayoutCreate(comm, &vlayout));
729: PetscCall(PetscLayoutSetSize(vlayout, NVertices));
730: PetscCall(PetscLayoutSetBlockSize(vlayout, 1));
731: PetscCall(PetscLayoutSetUp(vlayout));
732: PetscCall(PetscLayoutGetRange(vlayout, &vStart, &vEnd));
733: PetscCall(PetscLayoutDestroy(&vlayout));
734: numVertices = vEnd - vStart;
735: PetscCall(PetscMalloc1(PetscMax(1, numVertices) * 4, &coords));
736: for (PetscInt v = vStart; v < vEnd; ++v) {
737: coords[(v - vStart) * 4 + 0] = (PetscReal)(v % (N + 1));
738: coords[(v - vStart) * 4 + 1] = (PetscReal)(v / (N + 1));
739: coords[(v - vStart) * 4 + 2] = 0.;
740: coords[(v - vStart) * 4 + 3] = 0.;
741: }
742: PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, 2, numCells, numVertices, NVertices, 4, PETSC_TRUE, cells, 4, coords, NULL, &cellsSaved, &dm));
744: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
745: ierr = DMPlexGetOrdering(dm, DMPLEXCURVEMORTON, NULL, &perm);
746: PetscCall(PetscPopErrorHandler());
747: PetscCheck(ierr == PETSC_ERR_ARG_OUTOFRANGE, comm, PETSC_ERR_PLIB, "DMPlexGetOrdering() returned %d for a mesh in four dimensions, not PETSC_ERR_ARG_OUTOFRANGE", (int)ierr);
748: PetscCheck(!perm, comm, PETSC_ERR_PLIB, "The rejected call returned a permutation");
749: PetscCall(PetscPrintf(comm, "HighCoordinateDim: four coordinates per vertex rejected\n"));
750: PetscCall(PetscFree(cellsSaved));
751: PetscCall(DMDestroy(&dm));
752: PetscCall(PetscFree(coords));
753: PetscCall(PetscFree(cells));
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: // Part J: the sample-index numerator passes 2^31. DMPlexZCodeSelectSplitters() forms that product in
758: // 64 bits; a 32-bit product would wrap to a negative index, which its range check reports. Reaching
759: // that scale needs 32*size*NCells above 2^31, which is 4.3 million cells on 16 processes and 35
760: // million on two. That costs either processes or memory, so this part is off by default.
761: // The sample_overflow suffix turns it on. Run it by hand with, for example:
762: //
763: // mpiexec -n 32 ./ex106 -n 4 -sample_overflow -overflow_cells 2500000
764: static PetscErrorCode TestSampleOverflow(MPI_Comm comm, PetscInt NCells)
765: {
766: PetscSF sf;
767: PetscReal *cent, *newcent;
768: PetscInt64 maxProduct;
769: PetscInt *tags, *newtags;
770: PetscInt numCells, newNumCells, lo, hi, tot;
771: PetscMPIInt size, rank;
773: PetscFunctionBeginUser;
774: PetscCallMPI(MPI_Comm_size(comm, &size));
775: PetscCallMPI(MPI_Comm_rank(comm, &rank));
776: // Every cell starts on rank 0, the distribution that asks one rank for every sample. That rank
777: // then takes 32*size samples, and the largest sample-index numerator is that count minus one
778: // times the local cell count.
779: numCells = rank == 0 ? NCells : 0;
780: maxProduct = (32 * (PetscInt64)size - 1) * (PetscInt64)NCells;
781: PetscCheck(maxProduct > 2147483647, comm, PETSC_ERR_ARG_OUTOFRANGE, "The largest sample-index product here is %" PetscInt64_FMT ", which does not pass 2^31. Use more processes or more cells: 32*size*cells must pass 2^31", maxProduct);
782: PetscCall(PetscMalloc2(PetscMax(1, numCells), ¢, PetscMax(1, numCells), &tags));
783: // One axis, one cell per unit. The curve quantizes to 21 bits, so cells beyond 2^21 share a code
784: // and the global cell number separates them.
785: for (PetscInt c = 0; c < numCells; ++c) cent[c] = (PetscReal)c;
786: PetscCall(SetGlobalCellTags(comm, numCells, tags));
787: PetscCall(DMPlexReorderCellListByCurveFromCentroids(comm, DMPLEXCURVEMORTON, 1, numCells, cent, &sf, &newNumCells));
788: PetscCall(PetscMalloc2(newNumCells, &newcent, newNumCells, &newtags));
789: PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, cent, newcent, MPI_REPLACE));
790: PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, cent, newcent, MPI_REPLACE));
791: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, tags, newtags, MPI_REPLACE));
792: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, tags, newtags, MPI_REPLACE));
793: for (PetscInt c = 0; c < newNumCells; ++c) PetscCheck(newtags[c] >= 0 && newtags[c] < NCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Migrated tag %" PetscInt_FMT " out of range [0, %" PetscInt_FMT ")", newtags[c], NCells);
794: PetscCall(CheckGloballyCurveSorted(comm, 1, newNumCells, newcent, PETSC_TRUE, newtags));
795: lo = newNumCells;
796: hi = newNumCells;
797: tot = newNumCells;
798: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &lo, 1, MPIU_INT, MPI_MIN, comm));
799: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hi, 1, MPIU_INT, MPI_MAX, comm));
800: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &tot, 1, MPIU_INT, MPI_SUM, comm));
801: PetscCheck(tot == NCells, comm, PETSC_ERR_PLIB, "Cell count changed from %" PetscInt_FMT " to %" PetscInt_FMT, NCells, tot);
802: // The exact pass hides splitter quality from the counts, so also verify the complete
803: // lexicographic order above.
804: PetscCheck(lo > 0, comm, PETSC_ERR_PLIB, "A rank received no cells");
805: PetscCheck(hi - lo <= 1, comm, PETSC_ERR_PLIB, "Cells per rank run from %" PetscInt_FMT " to %" PetscInt_FMT ", which is not an exact split", lo, hi);
806: PetscCall(PetscPrintf(comm, "SampleOverflow: largest sample-index product %" PetscInt64_FMT " above 2^31, split balanced from %" PetscInt_FMT " to %" PetscInt_FMT " cells\n", maxProduct, lo, hi));
807: PetscCall(PetscSFDestroy(&sf));
808: PetscCall(PetscFree2(newcent, newtags));
809: PetscCall(PetscFree2(cent, tags));
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: // Part K: Morton ordering needs cell coordinates. A topology-only DMPLEX must be rejected before
814: // DMPlexGetCellCoordinates() reaches its coordinate-vector closure path.
815: static PetscErrorCode TestNoCoordinates(MPI_Comm comm)
816: {
817: DM dm;
818: IS perm = NULL;
819: PetscInt cone[4] = {1, 2, 3, 4};
820: PetscErrorCode ierr;
822: PetscFunctionBeginUser;
823: PetscCall(DMPlexCreate(comm, &dm));
824: PetscCall(DMSetDimension(dm, 2));
825: PetscCall(DMPlexSetChart(dm, 0, 5));
826: PetscCall(DMPlexSetConeSize(dm, 0, 4));
827: PetscCall(DMSetUp(dm));
828: PetscCall(DMPlexSetCone(dm, 0, cone));
829: PetscCall(DMPlexSymmetrize(dm));
830: PetscCall(DMPlexStratify(dm));
831: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
832: ierr = DMPlexGetOrdering(dm, DMPLEXCURVEMORTON, NULL, &perm);
833: PetscCall(PetscPopErrorHandler());
834: PetscCheck(ierr == PETSC_ERR_ARG_WRONGSTATE, comm, PETSC_ERR_PLIB, "DMPlexGetOrdering() returned %d for a mesh without coordinates, not PETSC_ERR_ARG_WRONGSTATE", (int)ierr);
835: PetscCheck(!perm, comm, PETSC_ERR_PLIB, "The rejected call returned a permutation");
836: PetscCall(DMDestroy(&dm));
837: PetscCall(PetscPrintf(comm, "NoCoordinates: Morton ordering rejected a mesh without coordinates\n"));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: int main(int argc, char **argv)
842: {
843: MPI_Comm comm;
844: PetscBool overflow = PETSC_FALSE;
845: PetscInt N = 4, overflowCells = 2500000;
847: PetscFunctionBeginUser;
848: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
849: comm = PETSC_COMM_WORLD;
850: PetscOptionsBegin(comm, "", "SFC cell-list reorder test options", "DMPLEX");
851: PetscCall(PetscOptionsInt("-n", "Cells per side", "ex106.c", N, &N, NULL));
852: PetscCall(PetscOptionsBool("-sample_overflow", "Run the sample-index case, which needs many processes and millions of cells", "ex106.c", overflow, &overflow, NULL));
853: PetscCall(PetscOptionsInt("-overflow_cells", "Cells for -sample_overflow", "ex106.c", overflowCells, &overflowCells, NULL));
854: PetscOptionsEnd();
855: PetscCheck(N > 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "-n must be positive");
856: PetscCheck(overflowCells > 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "-overflow_cells must be positive");
857: PetscCall(TestFromCentroids(comm, N, PETSC_FALSE));
858: PetscCall(TestFromCentroids(comm, N, PETSC_TRUE));
859: PetscCall(TestCellList(comm, N));
860: PetscCall(TestLocalityImproves(comm, N));
861: PetscCall(TestDegenerateGeometry(comm, N));
862: PetscCall(TestUnknownCurveType(comm));
863: PetscCall(TestOneDimensional(comm, N));
864: PetscCall(TestArgumentValidation(comm));
865: PetscCall(TestNoCells(comm));
866: PetscCall(TestHighCoordinateDim(comm));
867: PetscCall(TestNoCoordinates(comm));
868: if (overflow) PetscCall(TestSampleOverflow(comm, overflowCells));
869: PetscCall(PetscFinalize());
870: return 0;
871: }
873: /*TEST
875: test:
876: suffix: 0
877: nsize: {{1 2 3 4}}
878: args: -n 8
880: test:
881: suffix: 1
882: nsize: {{2 5 7}}
883: args: -n 16
885: # Small and empty-rank cases. At 2 the locality check has only two cells per process. At 8 the
886: # three-dimensional parts give every rank one cell. At 16 those parts leave half the ranks empty,
887: # which exercises collective key checks with locally NULL tag arrays.
888: test:
889: suffix: empty_ranks
890: nsize: {{2 8 16}}
891: args: -n 2
893: # Every cell starts on rank 0, the distribution a serial read produces. Check that the output is
894: # still a sorted, balanced permutation when the first exchange starts from one process.
895: test:
896: suffix: skewed_input
897: nsize: {{2 4 8}}
898: args: -n 8
900: # Degenerate geometry: coincident centroids, one distant node, and a corner cell on every rank.
901: # Apply the permutation and lexicographic key-order invariants directly to these inputs.
902: test:
903: suffix: degenerate
904: nsize: {{2 4 8}}
905: args: -n 6
907: # The sample-index product passes 2^31. This needs 32*size*cells above 2^31, so it costs either
908: # processes or memory: 16 processes with 4.3 million cells takes 0.3 s and 200 MB on the loaded
909: # process. A 32-bit product would wrap to a negative index at this scale, so this guards the 64-bit
910: # cast in DMPlexZCodeSelectSplitters(). A build with 64-bit indices cannot form a 32-bit product,
911: # so it skips the test and keeps the memory.
912: test:
913: suffix: sample_overflow
914: nsize: 16
915: requires: !defined(PETSC_USE_64BIT_INDICES)
916: args: -n 4 -sample_overflow -overflow_cells 4300000
918: TEST*/