Actual source code: plexpartition.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/partitionerimpl.h>
3: #include <petsc/private/matmetisimpl.h>
4: #include <petsc/private/matparmetisimpl.h>
5: #include <petsc/private/hashseti.h>
7: const char *const DMPlexCSRAlgorithms[] = {"mat", "graph", "overlap", "DMPlexCSRAlgorithm", "DM_PLEX_CSR_", NULL};
9: static inline PetscInt DMPlex_GlobalID(PetscInt point)
10: {
11: return point >= 0 ? point : -(point + 1);
12: }
14: static PetscErrorCode DMPlexCreatePartitionerGraph_Overlap(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
15: {
16: DM ovdm;
17: PetscSF sfPoint;
18: IS cellNumbering;
19: const PetscInt *cellNum;
20: PetscInt *adj = NULL, *vOffsets = NULL, *vAdj = NULL;
21: PetscBool useCone, useClosure;
22: PetscInt dim, depth, overlap, cStart, cEnd, c, v;
23: PetscMPIInt rank, size;
25: PetscFunctionBegin;
26: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
27: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
28: PetscCall(DMGetDimension(dm, &dim));
29: PetscCall(DMPlexGetDepth(dm, &depth));
30: if (dim != depth) {
31: /* We do not handle the uninterpolated case here */
32: PetscCall(DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency));
33: /* DMPlexCreateNeighborCSR does not make a numbering */
34: if (globalNumbering) PetscCall(DMPlexCreateCellNumbering(dm, PETSC_TRUE, globalNumbering));
35: /* Different behavior for empty graphs */
36: if (!*numVertices) {
37: PetscCall(PetscMalloc1(1, offsets));
38: (*offsets)[0] = 0;
39: }
40: /* Broken in parallel */
41: if (rank) PetscCheck(!*numVertices, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
44: /* Always use FVM adjacency to create partitioner graph */
45: PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
46: PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
47: /* Need overlap >= 1 */
48: PetscCall(DMPlexGetOverlap(dm, &overlap));
49: if (size && overlap < 1) {
50: PetscCall(DMPlexDistributeOverlap(dm, 1, NULL, &ovdm));
51: } else {
52: PetscCall(PetscObjectReference((PetscObject)dm));
53: ovdm = dm;
54: }
55: PetscCall(DMGetPointSF(ovdm, &sfPoint));
56: PetscCall(DMPlexGetHeightStratum(ovdm, height, &cStart, &cEnd));
57: PetscCall(DMPlexCreateNumbering_Plex(ovdm, cStart, cEnd, 0, NULL, sfPoint, &cellNumbering));
58: if (globalNumbering) {
59: PetscCall(PetscObjectReference((PetscObject)cellNumbering));
60: *globalNumbering = cellNumbering;
61: }
62: PetscCall(ISGetIndices(cellNumbering, &cellNum));
63: /* Determine sizes */
64: for (*numVertices = 0, c = cStart; c < cEnd; ++c) {
65: /* Skip non-owned cells in parallel (ParMETIS expects no overlap) */
66: if (cellNum[c - cStart] < 0) continue;
67: (*numVertices)++;
68: }
69: PetscCall(PetscCalloc1(*numVertices + 1, &vOffsets));
70: for (c = cStart, v = 0; c < cEnd; ++c) {
71: PetscInt adjSize = PETSC_DETERMINE, a, vsize = 0;
73: if (cellNum[c - cStart] < 0) continue;
74: PetscCall(DMPlexGetAdjacency(ovdm, c, &adjSize, &adj));
75: for (a = 0; a < adjSize; ++a) {
76: const PetscInt point = adj[a];
77: if (point != c && cStart <= point && point < cEnd) ++vsize;
78: }
79: vOffsets[v + 1] = vOffsets[v] + vsize;
80: ++v;
81: }
82: /* Determine adjacency */
83: PetscCall(PetscMalloc1(vOffsets[*numVertices], &vAdj));
84: for (c = cStart, v = 0; c < cEnd; ++c) {
85: PetscInt adjSize = PETSC_DETERMINE, a, off = vOffsets[v];
87: /* Skip non-owned cells in parallel (ParMETIS expects no overlap) */
88: if (cellNum[c - cStart] < 0) continue;
89: PetscCall(DMPlexGetAdjacency(ovdm, c, &adjSize, &adj));
90: for (a = 0; a < adjSize; ++a) {
91: const PetscInt point = adj[a];
92: if (point != c && cStart <= point && point < cEnd) vAdj[off++] = DMPlex_GlobalID(cellNum[point - cStart]);
93: }
94: PetscCheck(off == vOffsets[v + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Offsets %" PetscInt_FMT " should be %" PetscInt_FMT, off, vOffsets[v + 1]);
95: /* Sort adjacencies (not strictly necessary) */
96: PetscCall(PetscSortInt(off - vOffsets[v], &vAdj[vOffsets[v]]));
97: ++v;
98: }
99: PetscCall(PetscFree(adj));
100: PetscCall(ISRestoreIndices(cellNumbering, &cellNum));
101: PetscCall(ISDestroy(&cellNumbering));
102: PetscCall(DMSetBasicAdjacency(dm, useCone, useClosure));
103: PetscCall(DMDestroy(&ovdm));
104: if (offsets) {
105: *offsets = vOffsets;
106: } else PetscCall(PetscFree(vOffsets));
107: if (adjacency) {
108: *adjacency = vAdj;
109: } else PetscCall(PetscFree(vAdj));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
114: {
115: PetscInt dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
116: PetscInt *adj = NULL, *vOffsets = NULL, *graph = NULL;
117: IS cellNumbering;
118: const PetscInt *cellNum;
119: PetscBool useCone, useClosure;
120: PetscSection section;
121: PetscSegBuffer adjBuffer;
122: PetscSF sfPoint;
123: PetscInt *adjCells = NULL, *remoteCells = NULL;
124: const PetscInt *local;
125: PetscInt nroots, nleaves, l;
126: PetscMPIInt rank;
128: PetscFunctionBegin;
129: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
130: PetscCall(DMGetDimension(dm, &dim));
131: PetscCall(DMPlexGetDepth(dm, &depth));
132: if (dim != depth) {
133: /* We do not handle the uninterpolated case here */
134: PetscCall(DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency));
135: /* DMPlexCreateNeighborCSR does not make a numbering */
136: if (globalNumbering) PetscCall(DMPlexCreateCellNumbering(dm, PETSC_TRUE, globalNumbering));
137: /* Different behavior for empty graphs */
138: if (!*numVertices) {
139: PetscCall(PetscMalloc1(1, offsets));
140: (*offsets)[0] = 0;
141: }
142: /* Broken in parallel */
143: if (rank) PetscCheck(!*numVertices, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
146: PetscCall(DMGetPointSF(dm, &sfPoint));
147: PetscCall(DMPlexGetHeightStratum(dm, height, &pStart, &pEnd));
148: /* Build adjacency graph via a section/segbuffer */
149: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
150: PetscCall(PetscSectionSetChart(section, pStart, pEnd));
151: PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 1000, &adjBuffer));
152: /* Always use FVM adjacency to create partitioner graph */
153: PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
154: PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
155: PetscCall(DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering));
156: if (globalNumbering) {
157: PetscCall(PetscObjectReference((PetscObject)cellNumbering));
158: *globalNumbering = cellNumbering;
159: }
160: PetscCall(ISGetIndices(cellNumbering, &cellNum));
161: /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
162: Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
163: */
164: PetscCall(PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL));
165: if (nroots >= 0) {
166: PetscInt fStart, fEnd, f;
168: PetscCall(PetscCalloc2(nroots, &adjCells, nroots, &remoteCells));
169: PetscCall(DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd));
170: for (l = 0; l < nroots; ++l) adjCells[l] = -3;
171: for (f = fStart; f < fEnd; ++f) {
172: const PetscInt *support;
173: PetscInt supportSize;
175: PetscCall(DMPlexGetSupport(dm, f, &support));
176: PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
177: if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
178: else if (supportSize == 2) {
179: PetscCall(PetscFindInt(support[0], nleaves, local, &p));
180: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
181: PetscCall(PetscFindInt(support[1], nleaves, local, &p));
182: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
183: }
184: /* Handle non-conforming meshes */
185: if (supportSize > 2) {
186: PetscInt numChildren;
187: const PetscInt *children;
189: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, &children));
190: for (PetscInt i = 0; i < numChildren; ++i) {
191: const PetscInt child = children[i];
192: if (fStart <= child && child < fEnd) {
193: PetscCall(DMPlexGetSupport(dm, child, &support));
194: PetscCall(DMPlexGetSupportSize(dm, child, &supportSize));
195: if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
196: else if (supportSize == 2) {
197: PetscCall(PetscFindInt(support[0], nleaves, local, &p));
198: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
199: PetscCall(PetscFindInt(support[1], nleaves, local, &p));
200: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
201: }
202: }
203: }
204: }
205: }
206: for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
207: PetscCall(PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE));
208: PetscCall(PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE));
209: PetscCall(PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX));
210: PetscCall(PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX));
211: }
212: /* Combine local and global adjacencies */
213: for (*numVertices = 0, p = pStart; p < pEnd; p++) {
214: /* Skip non-owned cells in parallel (ParMETIS expects no overlap) */
215: if (nroots > 0) {
216: if (cellNum[p - pStart] < 0) continue;
217: }
218: /* Add remote cells */
219: if (remoteCells) {
220: const PetscInt gp = DMPlex_GlobalID(cellNum[p - pStart]);
221: PetscInt coneSize, numChildren, c, i;
222: const PetscInt *cone, *children;
224: PetscCall(DMPlexGetCone(dm, p, &cone));
225: PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
226: for (c = 0; c < coneSize; ++c) {
227: const PetscInt point = cone[c];
228: if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
229: PetscInt *PETSC_RESTRICT pBuf;
230: PetscCall(PetscSectionAddDof(section, p, 1));
231: PetscCall(PetscSegBufferGetInts(adjBuffer, 1, &pBuf));
232: *pBuf = remoteCells[point];
233: }
234: /* Handle non-conforming meshes */
235: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, &children));
236: for (i = 0; i < numChildren; ++i) {
237: const PetscInt child = children[i];
238: if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
239: PetscInt *PETSC_RESTRICT pBuf;
240: PetscCall(PetscSectionAddDof(section, p, 1));
241: PetscCall(PetscSegBufferGetInts(adjBuffer, 1, &pBuf));
242: *pBuf = remoteCells[child];
243: }
244: }
245: }
246: }
247: /* Add local cells */
248: adjSize = PETSC_DETERMINE;
249: PetscCall(DMPlexGetAdjacency(dm, p, &adjSize, &adj));
250: for (a = 0; a < adjSize; ++a) {
251: const PetscInt point = adj[a];
252: if (point != p && pStart <= point && point < pEnd) {
253: PetscInt *PETSC_RESTRICT pBuf;
254: PetscCall(PetscSectionAddDof(section, p, 1));
255: PetscCall(PetscSegBufferGetInts(adjBuffer, 1, &pBuf));
256: *pBuf = DMPlex_GlobalID(cellNum[point - pStart]);
257: }
258: }
259: (*numVertices)++;
260: }
261: PetscCall(PetscFree(adj));
262: PetscCall(PetscFree2(adjCells, remoteCells));
263: PetscCall(DMSetBasicAdjacency(dm, useCone, useClosure));
265: /* Derive CSR graph from section/segbuffer */
266: PetscCall(PetscSectionSetUp(section));
267: PetscCall(PetscSectionGetStorageSize(section, &size));
268: PetscCall(PetscMalloc1(*numVertices + 1, &vOffsets));
269: for (idx = 0, p = pStart; p < pEnd; p++) {
270: if (nroots > 0) {
271: if (cellNum[p - pStart] < 0) continue;
272: }
273: PetscCall(PetscSectionGetOffset(section, p, &vOffsets[idx++]));
274: }
275: vOffsets[*numVertices] = size;
276: PetscCall(PetscSegBufferExtractAlloc(adjBuffer, &graph));
278: if (nroots >= 0) {
279: /* Filter out duplicate edges using section/segbuffer */
280: PetscCall(PetscSectionReset(section));
281: PetscCall(PetscSectionSetChart(section, 0, *numVertices));
282: for (p = 0; p < *numVertices; p++) {
283: PetscInt start = vOffsets[p], end = vOffsets[p + 1];
284: PetscInt numEdges = end - start, *PETSC_RESTRICT edges;
285: PetscCall(PetscSortRemoveDupsInt(&numEdges, &graph[start]));
286: PetscCall(PetscSectionSetDof(section, p, numEdges));
287: PetscCall(PetscSegBufferGetInts(adjBuffer, numEdges, &edges));
288: PetscCall(PetscArraycpy(edges, &graph[start], numEdges));
289: }
290: PetscCall(PetscFree(vOffsets));
291: PetscCall(PetscFree(graph));
292: /* Derive CSR graph from section/segbuffer */
293: PetscCall(PetscSectionSetUp(section));
294: PetscCall(PetscSectionGetStorageSize(section, &size));
295: PetscCall(PetscMalloc1(*numVertices + 1, &vOffsets));
296: for (idx = 0, p = 0; p < *numVertices; p++) PetscCall(PetscSectionGetOffset(section, p, &vOffsets[idx++]));
297: vOffsets[*numVertices] = size;
298: PetscCall(PetscSegBufferExtractAlloc(adjBuffer, &graph));
299: } else {
300: /* Sort adjacencies (not strictly necessary) */
301: for (p = 0; p < *numVertices; p++) {
302: PetscInt start = vOffsets[p], end = vOffsets[p + 1];
303: PetscCall(PetscSortInt(end - start, &graph[start]));
304: }
305: }
307: if (offsets) *offsets = vOffsets;
308: if (adjacency) *adjacency = graph;
310: /* Cleanup */
311: PetscCall(PetscSegBufferDestroy(&adjBuffer));
312: PetscCall(PetscSectionDestroy(§ion));
313: PetscCall(ISRestoreIndices(cellNumbering, &cellNum));
314: PetscCall(ISDestroy(&cellNumbering));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
319: {
320: Mat conn, CSR;
321: IS fis, cis, cis_own;
322: PetscSF sfPoint;
323: const PetscInt *rows, *cols, *ii, *jj;
324: PetscInt *idxs, *idxs2;
325: PetscInt dim, depth, floc, cloc, i, M, N, c, m, cStart, cEnd, fStart, fEnd;
326: PetscMPIInt rank;
327: PetscBool flg;
329: PetscFunctionBegin;
330: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
331: PetscCall(DMGetDimension(dm, &dim));
332: PetscCall(DMPlexGetDepth(dm, &depth));
333: if (dim != depth) {
334: /* We do not handle the uninterpolated case here */
335: PetscCall(DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency));
336: /* DMPlexCreateNeighborCSR does not make a numbering */
337: if (globalNumbering) PetscCall(DMPlexCreateCellNumbering(dm, PETSC_TRUE, globalNumbering));
338: /* Different behavior for empty graphs */
339: if (!*numVertices) {
340: PetscCall(PetscMalloc1(1, offsets));
341: (*offsets)[0] = 0;
342: }
343: /* Broken in parallel */
344: if (rank) PetscCheck(!*numVertices, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
347: /* Interpolated and parallel case */
349: /* numbering */
350: PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sfPoint));
351: PetscCall(DMPlexGetHeightStratum(dm, height, &cStart, &cEnd));
352: PetscCall(DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd));
353: PetscCall(DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis));
354: PetscCall(DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis));
355: if (globalNumbering) PetscCall(ISDuplicate(cis, globalNumbering));
357: /* get positive global ids and local sizes for facets and cells */
358: PetscCall(ISGetLocalSize(fis, &m));
359: PetscCall(ISGetIndices(fis, &rows));
360: PetscCall(PetscMalloc1(m, &idxs));
361: for (i = 0, floc = 0; i < m; i++) {
362: const PetscInt p = rows[i];
364: if (p < 0) {
365: idxs[i] = -(p + 1);
366: } else {
367: idxs[i] = p;
368: floc += 1;
369: }
370: }
371: PetscCall(ISRestoreIndices(fis, &rows));
372: PetscCall(ISDestroy(&fis));
373: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis));
375: PetscCall(ISGetLocalSize(cis, &m));
376: PetscCall(ISGetIndices(cis, &cols));
377: PetscCall(PetscMalloc1(m, &idxs));
378: PetscCall(PetscMalloc1(m, &idxs2));
379: for (i = 0, cloc = 0; i < m; i++) {
380: const PetscInt p = cols[i];
382: if (p < 0) {
383: idxs[i] = -(p + 1);
384: } else {
385: idxs[i] = p;
386: idxs2[cloc++] = p;
387: }
388: }
389: PetscCall(ISRestoreIndices(cis, &cols));
390: PetscCall(ISDestroy(&cis));
391: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis));
392: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own));
394: /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
395: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &conn));
396: PetscCall(MatSetSizes(conn, floc, cloc, M, N));
397: PetscCall(MatSetType(conn, MATMPIAIJ));
398: PetscCall(DMPlexGetMaxSizes(dm, NULL, &m));
399: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
400: PetscCall(MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL));
402: /* Assemble matrix */
403: PetscCall(ISGetIndices(fis, &rows));
404: PetscCall(ISGetIndices(cis, &cols));
405: for (c = cStart; c < cEnd; c++) {
406: const PetscInt *cone;
407: PetscInt coneSize, row, col, f;
409: col = cols[c - cStart];
410: PetscCall(DMPlexGetCone(dm, c, &cone));
411: PetscCall(DMPlexGetConeSize(dm, c, &coneSize));
412: for (f = 0; f < coneSize; f++) {
413: const PetscScalar v = 1.0;
414: const PetscInt *children;
415: PetscInt numChildren;
417: row = rows[cone[f] - fStart];
418: PetscCall(MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES));
420: /* non-conforming meshes */
421: PetscCall(DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children));
422: for (PetscInt ch = 0; ch < numChildren; ch++) {
423: const PetscInt child = children[ch];
425: if (child < fStart || child >= fEnd) continue;
426: row = rows[child - fStart];
427: PetscCall(MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES));
428: }
429: }
430: }
431: PetscCall(ISRestoreIndices(fis, &rows));
432: PetscCall(ISRestoreIndices(cis, &cols));
433: PetscCall(ISDestroy(&fis));
434: PetscCall(ISDestroy(&cis));
435: PetscCall(MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY));
436: PetscCall(MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY));
438: /* Get parallel CSR by doing conn^T * conn */
439: PetscCall(MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CSR));
440: PetscCall(MatDestroy(&conn));
442: /* extract local part of the CSR */
443: PetscCall(MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn));
444: PetscCall(MatDestroy(&CSR));
445: PetscCall(MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
446: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
448: /* get back requested output */
449: if (numVertices) *numVertices = m;
450: if (offsets) {
451: PetscCall(PetscCalloc1(m + 1, &idxs));
452: for (i = 1; i < m + 1; i++) idxs[i] = ii[i] - i; /* ParMETIS does not like self-connectivity */
453: *offsets = idxs;
454: }
455: if (adjacency) {
456: PetscCall(PetscMalloc1(ii[m] - m, &idxs));
457: PetscCall(ISGetIndices(cis_own, &rows));
458: for (i = 0, c = 0; i < m; i++) {
459: PetscInt j, g = rows[i];
461: for (j = ii[i]; j < ii[i + 1]; j++) {
462: if (jj[j] == g) continue; /* again, self-connectivity */
463: idxs[c++] = jj[j];
464: }
465: }
466: PetscCheck(c == ii[m] - m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, c, ii[m] - m);
467: PetscCall(ISRestoreIndices(cis_own, &rows));
468: *adjacency = idxs;
469: }
471: /* cleanup */
472: PetscCall(ISDestroy(&cis_own));
473: PetscCall(MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg));
474: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
475: PetscCall(MatDestroy(&conn));
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: /*@
480: DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner
482: Collective
484: Input Parameters:
485: + dm - The mesh `DM`
486: - height - Height of the strata from which to construct the graph
488: Output Parameters:
489: + numVertices - Number of vertices in the graph
490: . offsets - Point offsets in the graph
491: . adjacency - Point connectivity in the graph
492: - globalNumbering - A map from the local cell numbering to the global numbering used in "adjacency". Negative indicates that the cell is a duplicate from another process.
494: Options Database Key:
495: . -dm_plex_csr_alg (mat|graph|overlap) - Choose the algorithm for computing the CSR graph
497: Level: developer
499: Note:
500: The user can control the definition of adjacency for the mesh using `DMSetAdjacency()`. They should choose the combination appropriate for the function
501: representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().
503: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCSRAlgorithm`, `PetscPartitionerGetType()`, `PetscPartitionerCreate()`, `DMSetAdjacency()`
504: @*/
505: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
506: {
507: DMPlexCSRAlgorithm alg = DM_PLEX_CSR_GRAPH;
509: PetscFunctionBegin;
510: PetscCall(PetscOptionsGetEnum(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_csr_alg", DMPlexCSRAlgorithms, (PetscEnum *)&alg, NULL));
511: switch (alg) {
512: case DM_PLEX_CSR_MAT:
513: PetscCall(DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering));
514: break;
515: case DM_PLEX_CSR_GRAPH:
516: PetscCall(DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering));
517: break;
518: case DM_PLEX_CSR_OVERLAP:
519: PetscCall(DMPlexCreatePartitionerGraph_Overlap(dm, height, numVertices, offsets, adjacency, globalNumbering));
520: break;
521: }
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /*@
526: DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.
528: Collective
530: Input Parameters:
531: + dm - The `DMPLEX`
532: - cellHeight - The height of mesh points to treat as cells (default should be 0)
534: Output Parameters:
535: + numVertices - The number of local vertices in the graph, or cells in the mesh.
536: . offsets - The offset to the adjacency list for each cell
537: - adjacency - The adjacency list for all cells
539: Level: advanced
541: Note:
542: This is suitable for input to a mesh partitioner like ParMETIS.
544: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreate()`
545: @*/
546: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
547: {
548: const PetscInt maxFaceCases = 30;
549: PetscInt numFaceCases = 0;
550: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
551: PetscInt *off, *adj;
552: PetscInt *neighborCells = NULL;
553: PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
555: PetscFunctionBegin;
556: /* For parallel partitioning, I think you have to communicate supports */
557: PetscCall(DMGetDimension(dm, &dim));
558: cellDim = dim - cellHeight;
559: PetscCall(DMPlexGetDepth(dm, &depth));
560: PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
561: if (cEnd - cStart == 0) {
562: if (numVertices) *numVertices = 0;
563: if (offsets) *offsets = NULL;
564: if (adjacency) *adjacency = NULL;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
567: numCells = cEnd - cStart;
568: faceDepth = depth - cellHeight;
569: if (dim == depth) {
570: PetscInt f, fStart, fEnd;
572: PetscCall(PetscCalloc1(numCells + 1, &off));
573: /* Count neighboring cells */
574: PetscCall(DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd));
575: for (f = fStart; f < fEnd; ++f) {
576: const PetscInt *support;
577: PetscInt supportSize;
578: PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
579: PetscCall(DMPlexGetSupport(dm, f, &support));
580: if (supportSize == 2) {
581: PetscInt numChildren;
583: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
584: if (!numChildren) {
585: ++off[support[0] - cStart + 1];
586: ++off[support[1] - cStart + 1];
587: }
588: }
589: }
590: /* Prefix sum */
591: for (c = 1; c <= numCells; ++c) off[c] += off[c - 1];
592: if (adjacency) {
593: PetscInt *tmp;
595: PetscCall(PetscMalloc1(off[numCells], &adj));
596: PetscCall(PetscMalloc1(numCells + 1, &tmp));
597: PetscCall(PetscArraycpy(tmp, off, numCells + 1));
598: /* Get neighboring cells */
599: for (f = fStart; f < fEnd; ++f) {
600: const PetscInt *support;
601: PetscInt supportSize;
602: PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
603: PetscCall(DMPlexGetSupport(dm, f, &support));
604: if (supportSize == 2) {
605: PetscInt numChildren;
607: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
608: if (!numChildren) {
609: adj[tmp[support[0] - cStart]++] = support[1];
610: adj[tmp[support[1] - cStart]++] = support[0];
611: }
612: }
613: }
614: for (c = 0; c < cEnd - cStart; ++c) PetscAssert(tmp[c] == off[c + 1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %" PetscInt_FMT " != %" PetscInt_FMT " for cell %" PetscInt_FMT, tmp[c], off[c], c + cStart);
615: PetscCall(PetscFree(tmp));
616: }
617: if (numVertices) *numVertices = numCells;
618: if (offsets) *offsets = off;
619: if (adjacency) *adjacency = adj;
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
622: /* Setup face recognition */
623: if (faceDepth == 1) {
624: PetscInt cornersSeen[30] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; /* Could use PetscBT */
626: for (c = cStart; c < cEnd; ++c) {
627: PetscInt corners;
629: PetscCall(DMPlexGetConeSize(dm, c, &corners));
630: if (!cornersSeen[corners]) {
631: PetscInt nFV;
633: PetscCheck(numFaceCases < maxFaceCases, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
634: cornersSeen[corners] = 1;
636: PetscCall(DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV));
638: numFaceVertices[numFaceCases++] = nFV;
639: }
640: }
641: }
642: PetscCall(PetscCalloc1(numCells + 1, &off));
643: /* Count neighboring cells */
644: for (cell = cStart; cell < cEnd; ++cell) {
645: PetscInt numNeighbors = PETSC_DETERMINE, n;
647: PetscCall(DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells));
648: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
649: for (n = 0; n < numNeighbors; ++n) {
650: PetscInt cellPair[2];
651: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
652: PetscInt meetSize = 0;
653: const PetscInt *meet = NULL;
655: cellPair[0] = cell;
656: cellPair[1] = neighborCells[n];
657: if (cellPair[0] == cellPair[1]) continue;
658: if (!found) {
659: PetscCall(DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet));
660: if (meetSize) {
661: for (PetscInt f = 0; f < numFaceCases; ++f) {
662: if (numFaceVertices[f] == meetSize) {
663: found = PETSC_TRUE;
664: break;
665: }
666: }
667: }
668: PetscCall(DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet));
669: }
670: if (found) ++off[cell - cStart + 1];
671: }
672: }
673: /* Prefix sum */
674: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell - 1];
676: if (adjacency) {
677: PetscCall(PetscMalloc1(off[numCells], &adj));
678: /* Get neighboring cells */
679: for (cell = cStart; cell < cEnd; ++cell) {
680: PetscInt numNeighbors = PETSC_DETERMINE, n;
681: PetscInt cellOffset = 0;
683: PetscCall(DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells));
684: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
685: for (n = 0; n < numNeighbors; ++n) {
686: PetscInt cellPair[2];
687: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
688: PetscInt meetSize = 0;
689: const PetscInt *meet = NULL;
691: cellPair[0] = cell;
692: cellPair[1] = neighborCells[n];
693: if (cellPair[0] == cellPair[1]) continue;
694: if (!found) {
695: PetscCall(DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet));
696: if (meetSize) {
697: for (PetscInt f = 0; f < numFaceCases; ++f) {
698: if (numFaceVertices[f] == meetSize) {
699: found = PETSC_TRUE;
700: break;
701: }
702: }
703: }
704: PetscCall(DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet));
705: }
706: if (found) {
707: adj[off[cell - cStart] + cellOffset] = neighborCells[n];
708: ++cellOffset;
709: }
710: }
711: }
712: }
713: PetscCall(PetscFree(neighborCells));
714: if (numVertices) *numVertices = numCells;
715: if (offsets) *offsets = off;
716: if (adjacency) *adjacency = adj;
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*@
721: PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh
723: Collective
725: Input Parameters:
726: + part - The `PetscPartitioner`
727: . targetSection - The `PetscSection` describing the absolute weight of each partition (can be `NULL`)
728: - dm - The mesh `DM`
730: Output Parameters:
731: + partSection - The `PetscSection` giving the division of points by partition
732: - partition - The list of points by partition
734: Level: developer
736: Note:
737: If the `DM` has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
738: by the section in the transitive closure of the point.
740: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `PetscPartitionerCreate()`, `PetscSectionCreate()`,
741: `PetscSectionSetChart()`, `PetscPartitionerPartition()`
742: @*/
743: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
744: {
745: PetscMPIInt size;
746: PetscBool isplex;
747: PetscSection vertSection = NULL, edgeSection = NULL;
749: PetscFunctionBegin;
754: PetscAssertPointer(partition, 5);
755: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
756: PetscCheck(isplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)dm)->type_name);
757: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)part), &size));
758: if (size == 1) {
759: PetscInt *points;
760: PetscInt cStart, cEnd, c;
762: PetscCall(DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd));
763: PetscCall(PetscSectionReset(partSection));
764: PetscCall(PetscSectionSetChart(partSection, 0, size));
765: PetscCall(PetscSectionSetDof(partSection, 0, cEnd - cStart));
766: PetscCall(PetscSectionSetUp(partSection));
767: PetscCall(PetscMalloc1(cEnd - cStart, &points));
768: for (c = cStart; c < cEnd; ++c) points[c] = c;
769: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), cEnd - cStart, points, PETSC_OWN_POINTER, partition));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
772: if (part->height == 0) {
773: PetscInt numVertices = 0;
774: PetscInt *start = NULL;
775: PetscInt *adjacency = NULL;
776: IS globalNumbering;
778: if (!part->noGraph || part->viewerGraph) {
779: PetscCall(DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering));
780: } else { /* only compute the number of owned local vertices */
781: const PetscInt *idxs;
782: PetscInt p, pStart, pEnd;
784: PetscCall(DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd));
785: PetscCall(DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering));
786: PetscCall(ISGetIndices(globalNumbering, &idxs));
787: for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
788: PetscCall(ISRestoreIndices(globalNumbering, &idxs));
789: }
790: if (part->usevwgt) {
791: PetscSection section = dm->localSection, clSection = NULL;
792: IS clPoints = NULL;
793: const PetscInt *gid, *clIdx;
794: PetscInt v, p, pStart, pEnd;
796: /* dm->localSection encodes degrees of freedom per point, not per cell. We need to get the closure index to properly specify cell weights (aka dofs) */
797: /* We do this only if the local section has been set */
798: if (section) {
799: PetscCall(PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL));
800: if (!clSection) PetscCall(DMPlexCreateClosureIndex(dm, NULL));
801: PetscCall(PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints));
802: PetscCall(ISGetIndices(clPoints, &clIdx));
803: }
804: PetscCall(DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd));
805: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &vertSection));
806: PetscCall(PetscSectionSetChart(vertSection, 0, numVertices));
807: if (globalNumbering) PetscCall(ISGetIndices(globalNumbering, &gid));
808: else gid = NULL;
809: for (p = pStart, v = 0; p < pEnd; ++p) {
810: PetscInt dof = 1;
812: /* skip cells in the overlap */
813: if (gid && gid[p - pStart] < 0) continue;
815: if (section) {
816: PetscInt cl, clSize, clOff;
818: dof = 0;
819: PetscCall(PetscSectionGetDof(clSection, p, &clSize));
820: PetscCall(PetscSectionGetOffset(clSection, p, &clOff));
821: for (cl = 0; cl < clSize; cl += 2) {
822: PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */
824: PetscCall(PetscSectionGetDof(section, clPoint, &clDof));
825: dof += clDof;
826: }
827: }
828: PetscCheck(dof, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of dofs for point %" PetscInt_FMT " in the local section should be positive", p);
829: PetscCall(PetscSectionSetDof(vertSection, v, dof));
830: v++;
831: }
832: if (globalNumbering) PetscCall(ISRestoreIndices(globalNumbering, &gid));
833: if (clPoints) PetscCall(ISRestoreIndices(clPoints, &clIdx));
834: PetscCall(PetscSectionSetUp(vertSection));
835: }
836: if (part->useewgt) {
837: const PetscInt numEdges = start[numVertices];
839: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &edgeSection));
840: PetscCall(PetscSectionSetChart(edgeSection, 0, numEdges));
841: for (PetscInt e = 0; e < start[numVertices]; ++e) PetscCall(PetscSectionSetDof(edgeSection, e, 1));
842: for (PetscInt v = 0; v < numVertices; ++v) {
843: DMPolytopeType ct;
845: // Assume v is the cell number
846: PetscCall(DMPlexGetCellType(dm, v, &ct));
847: if (ct != DM_POLYTOPE_POINT_PRISM_TENSOR && ct != DM_POLYTOPE_SEG_PRISM_TENSOR && ct != DM_POLYTOPE_TRI_PRISM_TENSOR && ct != DM_POLYTOPE_QUAD_PRISM_TENSOR) continue;
849: for (PetscInt e = start[v]; e < start[v + 1]; ++e) PetscCall(PetscSectionSetDof(edgeSection, e, 3));
850: }
851: PetscCall(PetscSectionSetUp(edgeSection));
852: }
853: PetscCall(PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, edgeSection, targetSection, partSection, partition));
854: PetscCall(PetscFree(start));
855: PetscCall(PetscFree(adjacency));
856: if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
857: const PetscInt *globalNum;
858: const PetscInt *partIdx;
859: PetscInt *map, cStart, cEnd;
860: PetscInt *adjusted, i, localSize, offset;
861: IS newPartition;
863: PetscCall(ISGetLocalSize(*partition, &localSize));
864: PetscCall(PetscMalloc1(localSize, &adjusted));
865: PetscCall(ISGetIndices(globalNumbering, &globalNum));
866: PetscCall(ISGetIndices(*partition, &partIdx));
867: PetscCall(PetscMalloc1(localSize, &map));
868: PetscCall(DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd));
869: for (i = cStart, offset = 0; i < cEnd; i++) {
870: if (globalNum[i - cStart] >= 0) map[offset++] = i;
871: }
872: for (i = 0; i < localSize; i++) adjusted[i] = map[partIdx[i]];
873: PetscCall(PetscFree(map));
874: PetscCall(ISRestoreIndices(*partition, &partIdx));
875: PetscCall(ISRestoreIndices(globalNumbering, &globalNum));
876: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, localSize, adjusted, PETSC_OWN_POINTER, &newPartition));
877: PetscCall(ISDestroy(&globalNumbering));
878: PetscCall(ISDestroy(partition));
879: *partition = newPartition;
880: }
881: } else SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %" PetscInt_FMT " for points to partition", part->height);
882: PetscCall(PetscSectionDestroy(&vertSection));
883: PetscCall(PetscSectionDestroy(&edgeSection));
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: /*@
888: DMPlexGetPartitioner - Get the mesh partitioner
890: Not Collective
892: Input Parameter:
893: . dm - The `DM`
895: Output Parameter:
896: . part - The `PetscPartitioner`
898: Level: developer
900: Note:
901: This gets a borrowed reference, so the user should not destroy this `PetscPartitioner`.
903: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `DMPlexSetPartitioner()`, `PetscPartitionerDMPlexPartition()`, `PetscPartitionerCreate()`
904: @*/
905: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
906: {
907: DM_Plex *mesh = (DM_Plex *)dm->data;
909: PetscFunctionBegin;
911: PetscAssertPointer(part, 2);
912: *part = mesh->partitioner;
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: /*@
917: DMPlexSetPartitioner - Set the mesh partitioner
919: logically Collective
921: Input Parameters:
922: + dm - The `DM`
923: - part - The partitioner
925: Level: developer
927: Note:
928: Any existing `PetscPartitioner` will be destroyed.
930: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `DMPlexDistribute()`, `DMPlexGetPartitioner()`, `PetscPartitionerCreate()`
931: @*/
932: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
933: {
934: DM_Plex *mesh = (DM_Plex *)dm->data;
936: PetscFunctionBegin;
939: PetscCall(PetscObjectReference((PetscObject)part));
940: PetscCall(PetscPartitionerDestroy(&mesh->partitioner));
941: mesh->partitioner = part;
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
946: {
947: const PetscInt *cone;
948: PetscInt coneSize, c;
949: PetscBool missing;
951: PetscFunctionBeginHot;
952: PetscCall(PetscHSetIQueryAdd(ht, point, &missing));
953: if (missing) {
954: PetscCall(DMPlexGetCone(dm, point, &cone));
955: PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
956: for (c = 0; c < coneSize; c++) PetscCall(DMPlexAddClosure_Private(dm, ht, cone[c]));
957: }
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
962: {
963: PetscFunctionBegin;
964: if (up) {
965: PetscInt parent;
967: PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
968: if (parent != point) {
969: PetscInt closureSize, *closure = NULL, i;
971: PetscCall(DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
972: for (i = 0; i < closureSize; i++) {
973: PetscInt cpoint = closure[2 * i];
975: PetscCall(PetscHSetIAdd(ht, cpoint));
976: PetscCall(DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_TRUE, PETSC_FALSE));
977: }
978: PetscCall(DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure));
979: }
980: }
981: if (down) {
982: PetscInt numChildren;
983: const PetscInt *children;
985: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, &children));
986: if (numChildren) {
987: for (PetscInt i = 0; i < numChildren; i++) {
988: PetscInt cpoint = children[i];
990: PetscCall(PetscHSetIAdd(ht, cpoint));
991: PetscCall(DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_FALSE, PETSC_TRUE));
992: }
993: }
994: }
995: PetscFunctionReturn(PETSC_SUCCESS);
996: }
998: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
999: {
1000: PetscInt parent;
1002: PetscFunctionBeginHot;
1003: PetscCall(DMPlexGetTreeParent(dm, point, &parent, NULL));
1004: if (point != parent) {
1005: const PetscInt *cone;
1006: PetscInt coneSize, c;
1008: PetscCall(DMPlexAddClosureTree_Up_Private(dm, ht, parent));
1009: PetscCall(DMPlexAddClosure_Private(dm, ht, parent));
1010: PetscCall(DMPlexGetCone(dm, parent, &cone));
1011: PetscCall(DMPlexGetConeSize(dm, parent, &coneSize));
1012: for (c = 0; c < coneSize; c++) {
1013: const PetscInt cp = cone[c];
1015: PetscCall(DMPlexAddClosureTree_Up_Private(dm, ht, cp));
1016: }
1017: }
1018: PetscFunctionReturn(PETSC_SUCCESS);
1019: }
1021: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
1022: {
1023: PetscInt numChildren;
1024: const PetscInt *children;
1026: PetscFunctionBeginHot;
1027: PetscCall(DMPlexGetTreeChildren(dm, point, &numChildren, &children));
1028: for (PetscInt i = 0; i < numChildren; i++) PetscCall(PetscHSetIAdd(ht, children[i]));
1029: PetscFunctionReturn(PETSC_SUCCESS);
1030: }
1032: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
1033: {
1034: const PetscInt *cone;
1035: PetscInt coneSize, c;
1037: PetscFunctionBeginHot;
1038: PetscCall(PetscHSetIAdd(ht, point));
1039: PetscCall(DMPlexAddClosureTree_Up_Private(dm, ht, point));
1040: PetscCall(DMPlexAddClosureTree_Down_Private(dm, ht, point));
1041: PetscCall(DMPlexGetCone(dm, point, &cone));
1042: PetscCall(DMPlexGetConeSize(dm, point, &coneSize));
1043: for (c = 0; c < coneSize; c++) PetscCall(DMPlexAddClosureTree_Private(dm, ht, cone[c]));
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
1048: {
1049: DM_Plex *mesh = (DM_Plex *)dm->data;
1050: const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
1051: PetscInt nelems, *elems, off = 0, p;
1052: PetscHSetI ht = NULL;
1054: PetscFunctionBegin;
1055: PetscCall(PetscHSetICreate(&ht));
1056: PetscCall(PetscHSetIResize(ht, numPoints * 16));
1057: if (!hasTree) {
1058: for (p = 0; p < numPoints; ++p) PetscCall(DMPlexAddClosure_Private(dm, ht, points[p]));
1059: } else {
1060: #if 1
1061: for (p = 0; p < numPoints; ++p) PetscCall(DMPlexAddClosureTree_Private(dm, ht, points[p]));
1062: #else
1063: PetscInt *closure = NULL, closureSize, c;
1064: for (p = 0; p < numPoints; ++p) {
1065: PetscCall(DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure));
1066: for (c = 0; c < closureSize * 2; c += 2) {
1067: PetscCall(PetscHSetIAdd(ht, closure[c]));
1068: if (hasTree) PetscCall(DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE));
1069: }
1070: }
1071: if (closure) PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure));
1072: #endif
1073: }
1074: PetscCall(PetscHSetIGetSize(ht, &nelems));
1075: PetscCall(PetscMalloc1(nelems, &elems));
1076: PetscCall(PetscHSetIGetElems(ht, &off, elems));
1077: PetscCall(PetscHSetIDestroy(&ht));
1078: PetscCall(PetscSortInt(nelems, elems));
1079: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: /*@
1084: DMPlexPartitionLabelClosure - Add the closure of all points to the partition label
1086: Input Parameters:
1087: + dm - The `DM`
1088: - label - `DMLabel` assigning ranks to remote roots
1090: Level: developer
1092: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`
1093: @*/
1094: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
1095: {
1096: IS rankIS, pointIS, closureIS;
1097: const PetscInt *ranks, *points;
1098: PetscInt numRanks, numPoints, r;
1100: PetscFunctionBegin;
1101: PetscCall(DMLabelGetValueIS(label, &rankIS));
1102: PetscCall(ISGetLocalSize(rankIS, &numRanks));
1103: PetscCall(ISGetIndices(rankIS, &ranks));
1104: for (r = 0; r < numRanks; ++r) {
1105: const PetscInt rank = ranks[r];
1106: PetscCall(DMLabelGetStratumIS(label, rank, &pointIS));
1107: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1108: PetscCall(ISGetIndices(pointIS, &points));
1109: PetscCall(DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS));
1110: PetscCall(ISRestoreIndices(pointIS, &points));
1111: PetscCall(ISDestroy(&pointIS));
1112: PetscCall(DMLabelSetStratumIS(label, rank, closureIS));
1113: PetscCall(ISDestroy(&closureIS));
1114: }
1115: PetscCall(ISRestoreIndices(rankIS, &ranks));
1116: PetscCall(ISDestroy(&rankIS));
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: /*@
1121: DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label
1123: Input Parameters:
1124: + dm - The `DM`
1125: - label - `DMLabel` assigning ranks to remote roots
1127: Level: developer
1129: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`
1130: @*/
1131: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
1132: {
1133: IS rankIS, pointIS;
1134: const PetscInt *ranks, *points;
1135: PetscInt numRanks, numPoints, r, p, a, adjSize;
1136: PetscInt *adj = NULL;
1138: PetscFunctionBegin;
1139: PetscCall(DMLabelGetValueIS(label, &rankIS));
1140: PetscCall(ISGetLocalSize(rankIS, &numRanks));
1141: PetscCall(ISGetIndices(rankIS, &ranks));
1142: for (r = 0; r < numRanks; ++r) {
1143: const PetscInt rank = ranks[r];
1145: PetscCall(DMLabelGetStratumIS(label, rank, &pointIS));
1146: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1147: PetscCall(ISGetIndices(pointIS, &points));
1148: for (p = 0; p < numPoints; ++p) {
1149: adjSize = PETSC_DETERMINE;
1150: PetscCall(DMPlexGetAdjacency(dm, points[p], &adjSize, &adj));
1151: for (a = 0; a < adjSize; ++a) PetscCall(DMLabelSetValue(label, adj[a], rank));
1152: }
1153: PetscCall(ISRestoreIndices(pointIS, &points));
1154: PetscCall(ISDestroy(&pointIS));
1155: }
1156: PetscCall(ISRestoreIndices(rankIS, &ranks));
1157: PetscCall(ISDestroy(&rankIS));
1158: PetscCall(PetscFree(adj));
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: /*@
1163: DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point `PetscSF`
1165: Input Parameters:
1166: + dm - The `DM`
1167: - label - `DMLabel` assigning ranks to remote roots
1169: Level: developer
1171: Note:
1172: This is required when generating multi-level overlaps to capture
1173: overlap points from non-neighboring partitions.
1175: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`
1176: @*/
1177: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
1178: {
1179: MPI_Comm comm;
1180: PetscMPIInt rank;
1181: PetscSF sfPoint;
1182: DMLabel lblRoots, lblLeaves;
1183: IS rankIS, pointIS;
1184: const PetscInt *ranks;
1185: PetscInt numRanks, r;
1187: PetscFunctionBegin;
1188: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1189: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1190: PetscCall(DMGetPointSF(dm, &sfPoint));
1191: /* Pull point contributions from remote leaves into local roots */
1192: PetscCall(DMLabelGather(label, sfPoint, &lblLeaves));
1193: PetscCall(DMLabelGetValueIS(lblLeaves, &rankIS));
1194: PetscCall(ISGetLocalSize(rankIS, &numRanks));
1195: PetscCall(ISGetIndices(rankIS, &ranks));
1196: for (r = 0; r < numRanks; ++r) {
1197: const PetscInt remoteRank = ranks[r];
1198: if (remoteRank == rank) continue;
1199: PetscCall(DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS));
1200: PetscCall(DMLabelInsertIS(label, pointIS, remoteRank));
1201: PetscCall(ISDestroy(&pointIS));
1202: }
1203: PetscCall(ISRestoreIndices(rankIS, &ranks));
1204: PetscCall(ISDestroy(&rankIS));
1205: PetscCall(DMLabelDestroy(&lblLeaves));
1206: /* Push point contributions from roots into remote leaves */
1207: PetscCall(DMLabelDistribute(label, sfPoint, &lblRoots));
1208: PetscCall(DMLabelGetValueIS(lblRoots, &rankIS));
1209: PetscCall(ISGetLocalSize(rankIS, &numRanks));
1210: PetscCall(ISGetIndices(rankIS, &ranks));
1211: for (r = 0; r < numRanks; ++r) {
1212: const PetscInt remoteRank = ranks[r];
1213: if (remoteRank == rank) continue;
1214: PetscCall(DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS));
1215: PetscCall(DMLabelInsertIS(label, pointIS, remoteRank));
1216: PetscCall(ISDestroy(&pointIS));
1217: }
1218: PetscCall(ISRestoreIndices(rankIS, &ranks));
1219: PetscCall(ISDestroy(&rankIS));
1220: PetscCall(DMLabelDestroy(&lblRoots));
1221: PetscFunctionReturn(PETSC_SUCCESS);
1222: }
1224: /*@
1225: DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label
1227: Input Parameters:
1228: + dm - The `DM`
1229: . rootLabel - `DMLabel` assigning ranks to local roots
1230: - processSF - A star forest mapping into the local index on each remote rank
1232: Output Parameter:
1233: . leafLabel - `DMLabel` assigning ranks to remote roots
1235: Level: developer
1237: Note:
1238: The rootLabel defines a send pattern by mapping local points to remote target ranks. The
1239: resulting leafLabel is a receiver mapping of remote roots to their parent rank.
1241: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`
1242: @*/
1243: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
1244: {
1245: MPI_Comm comm;
1246: PetscMPIInt rank, size, r;
1247: PetscInt p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
1248: PetscSF sfPoint;
1249: PetscSection rootSection;
1250: PetscSFNode *rootPoints, *leafPoints;
1251: const PetscSFNode *remote;
1252: const PetscInt *local, *neighbors;
1253: IS valueIS;
1254: PetscBool mpiOverflow = PETSC_FALSE;
1256: PetscFunctionBegin;
1257: PetscCall(PetscLogEventBegin(DMPLEX_PartLabelInvert, dm, 0, 0, 0));
1258: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1259: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1260: PetscCallMPI(MPI_Comm_size(comm, &size));
1261: PetscCall(DMGetPointSF(dm, &sfPoint));
1263: /* Convert to (point, rank) and use actual owners */
1264: PetscCall(PetscSectionCreate(comm, &rootSection));
1265: PetscCall(PetscSectionSetChart(rootSection, 0, size));
1266: PetscCall(DMLabelGetValueIS(rootLabel, &valueIS));
1267: PetscCall(ISGetLocalSize(valueIS, &numNeighbors));
1268: PetscCall(ISGetIndices(valueIS, &neighbors));
1269: for (n = 0; n < numNeighbors; ++n) {
1270: PetscCall(DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints));
1271: PetscCall(PetscSectionAddDof(rootSection, neighbors[n], numPoints));
1272: }
1273: PetscCall(PetscSectionSetUp(rootSection));
1274: PetscCall(PetscSectionGetStorageSize(rootSection, &rootSize));
1275: PetscCall(PetscMalloc1(rootSize, &rootPoints));
1276: PetscCall(PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote));
1277: for (n = 0; n < numNeighbors; ++n) {
1278: IS pointIS;
1279: const PetscInt *points;
1281: PetscCall(PetscSectionGetOffset(rootSection, neighbors[n], &off));
1282: PetscCall(DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS));
1283: PetscCall(ISGetLocalSize(pointIS, &numPoints));
1284: PetscCall(ISGetIndices(pointIS, &points));
1285: for (p = 0; p < numPoints; ++p) {
1286: if (local) PetscCall(PetscFindInt(points[p], nleaves, local, &l));
1287: else l = -1;
1288: if (l >= 0) {
1289: rootPoints[off + p] = remote[l];
1290: } else {
1291: rootPoints[off + p].index = points[p];
1292: rootPoints[off + p].rank = rank;
1293: }
1294: }
1295: PetscCall(ISRestoreIndices(pointIS, &points));
1296: PetscCall(ISDestroy(&pointIS));
1297: }
1299: /* Try to communicate overlap using All-to-All */
1300: if (!processSF) {
1301: PetscCount counter = 0;
1302: PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;
1304: PetscCall(PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls));
1305: for (n = 0; n < numNeighbors; ++n) {
1306: PetscCall(PetscSectionGetDof(rootSection, neighbors[n], &dof));
1307: PetscCall(PetscSectionGetOffset(rootSection, neighbors[n], &off));
1308: #if PetscDefined(USE_64BIT_INDICES)
1309: if (dof > PETSC_MPI_INT_MAX) {
1310: mpiOverflow = PETSC_TRUE;
1311: break;
1312: }
1313: if (off > PETSC_MPI_INT_MAX) {
1314: mpiOverflow = PETSC_TRUE;
1315: break;
1316: }
1317: #endif
1318: PetscCall(PetscMPIIntCast(dof, &scounts[neighbors[n]]));
1319: PetscCall(PetscMPIIntCast(off, &sdispls[neighbors[n]]));
1320: }
1321: PetscCallMPI(MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm));
1322: for (r = 0; r < size; ++r) {
1323: PetscCall(PetscMPIIntCast(counter, &rdispls[r]));
1324: counter += rcounts[r];
1325: }
1326: if (counter > PETSC_MPI_INT_MAX) mpiOverflow = PETSC_TRUE;
1327: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mpiOverflow, 1, MPI_C_BOOL, MPI_LOR, comm));
1328: if (!mpiOverflow) {
1329: PetscCall(PetscInfo(dm, "Using MPI_Alltoallv() for mesh distribution\n"));
1330: PetscCall(PetscIntCast(counter, &leafSize));
1331: PetscCall(PetscMalloc1(leafSize, &leafPoints));
1332: PetscCallMPI(MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_SF_NODE, leafPoints, rcounts, rdispls, MPIU_SF_NODE, comm));
1333: }
1334: PetscCall(PetscFree4(scounts, sdispls, rcounts, rdispls));
1335: }
1337: /* Communicate overlap using process star forest */
1338: if (processSF || mpiOverflow) {
1339: PetscSF procSF;
1340: PetscSection leafSection;
1342: if (processSF) {
1343: PetscCall(PetscInfo(dm, "Using processSF for mesh distribution\n"));
1344: PetscCall(PetscObjectReference((PetscObject)processSF));
1345: procSF = processSF;
1346: } else {
1347: PetscCall(PetscInfo(dm, "Using processSF for mesh distribution (MPI overflow)\n"));
1348: PetscCall(PetscSFCreate(comm, &procSF));
1349: PetscCall(PetscSFSetGraphWithPattern(procSF, NULL, PETSCSF_PATTERN_ALLTOALL));
1350: }
1352: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection));
1353: PetscCall(DMPlexDistributeData(dm, procSF, rootSection, MPIU_SF_NODE, rootPoints, leafSection, (void **)&leafPoints));
1354: PetscCall(PetscSectionGetStorageSize(leafSection, &leafSize));
1355: PetscCall(PetscSectionDestroy(&leafSection));
1356: PetscCall(PetscSFDestroy(&procSF));
1357: }
1359: for (p = 0; p < leafSize; p++) PetscCall(DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank));
1361: PetscCall(ISRestoreIndices(valueIS, &neighbors));
1362: PetscCall(ISDestroy(&valueIS));
1363: PetscCall(PetscSectionDestroy(&rootSection));
1364: PetscCall(PetscFree(rootPoints));
1365: PetscCall(PetscFree(leafPoints));
1366: PetscCall(PetscLogEventEnd(DMPLEX_PartLabelInvert, dm, 0, 0, 0));
1367: PetscFunctionReturn(PETSC_SUCCESS);
1368: }
1370: /*@
1371: DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points
1373: Input Parameters:
1374: + dm - The `DM`
1375: . label - `DMLabel` assigning ranks to remote roots
1376: - sortRanks - Whether or not to sort the `PetscSF` leaves by rank
1378: Output Parameter:
1379: . sf - The star forest communication context encapsulating the defined mapping
1381: Level: developer
1383: Note:
1384: The incoming label is a receiver mapping of remote points to their parent rank.
1386: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`
1387: @*/
1388: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscBool sortRanks, PetscSF *sf)
1389: {
1390: PetscMPIInt rank;
1391: PetscInt n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
1392: PetscSFNode *remotePoints;
1393: IS remoteRootIS, neighborsIS;
1394: const PetscInt *remoteRoots, *neighbors;
1395: PetscMPIInt myRank;
1397: PetscFunctionBegin;
1398: PetscCall(PetscLogEventBegin(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0));
1399: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1400: PetscCall(DMLabelGetValueIS(label, &neighborsIS));
1402: if (sortRanks) {
1403: IS is;
1405: PetscCall(ISDuplicate(neighborsIS, &is));
1406: PetscCall(ISSort(is));
1407: PetscCall(ISDestroy(&neighborsIS));
1408: neighborsIS = is;
1409: }
1410: myRank = sortRanks ? -1 : rank;
1412: PetscCall(ISGetLocalSize(neighborsIS, &nNeighbors));
1413: PetscCall(ISGetIndices(neighborsIS, &neighbors));
1414: for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
1415: PetscCall(DMLabelGetStratumSize(label, neighbors[n], &numPoints));
1416: numRemote += numPoints;
1417: }
1418: PetscCall(PetscMalloc1(numRemote, &remotePoints));
1420: /* Put owned points first if not sorting the ranks */
1421: if (!sortRanks) {
1422: PetscCall(DMLabelGetStratumSize(label, rank, &numPoints));
1423: if (numPoints > 0) {
1424: PetscCall(DMLabelGetStratumIS(label, rank, &remoteRootIS));
1425: PetscCall(ISGetIndices(remoteRootIS, &remoteRoots));
1426: for (p = 0; p < numPoints; p++) {
1427: remotePoints[idx].index = remoteRoots[p];
1428: remotePoints[idx].rank = rank;
1429: idx++;
1430: }
1431: PetscCall(ISRestoreIndices(remoteRootIS, &remoteRoots));
1432: PetscCall(ISDestroy(&remoteRootIS));
1433: }
1434: }
1436: /* Now add remote points */
1437: for (n = 0; n < nNeighbors; ++n) {
1438: PetscMPIInt nn;
1440: PetscCall(PetscMPIIntCast(neighbors[n], &nn));
1441: PetscCall(DMLabelGetStratumSize(label, nn, &numPoints));
1442: if (nn == myRank || numPoints <= 0) continue;
1443: PetscCall(DMLabelGetStratumIS(label, nn, &remoteRootIS));
1444: PetscCall(ISGetIndices(remoteRootIS, &remoteRoots));
1445: for (p = 0; p < numPoints; p++) {
1446: remotePoints[idx].index = remoteRoots[p];
1447: remotePoints[idx].rank = nn;
1448: idx++;
1449: }
1450: PetscCall(ISRestoreIndices(remoteRootIS, &remoteRoots));
1451: PetscCall(ISDestroy(&remoteRootIS));
1452: }
1454: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), sf));
1455: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1456: PetscCall(PetscSFSetGraph(*sf, pEnd - pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER));
1457: PetscCall(PetscSFSetUp(*sf));
1458: PetscCall(ISDestroy(&neighborsIS));
1459: PetscCall(PetscLogEventEnd(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0));
1460: PetscFunctionReturn(PETSC_SUCCESS);
1461: }
1463: #if PetscDefined(HAVE_PARMETIS)
1464: #include <parmetis.h>
1465: #endif
1467: /*
1468: DMPlexRewriteSF - Rewrites the ownership of the `PetscSF` of a `DM` (in place).
1470: Input parameters:b
1471: + dm - The `DMPLEX` object.
1472: . n - The number of points.
1473: . pointsToRewrite - The points in the `PetscSF` whose ownership will change.
1474: . targetOwners - New owner for each element in pointsToRewrite.
1475: - degrees - Degrees of the points in the `PetscSF` as obtained by `PetscSFComputeDegreeBegin()`/`PetscSFComputeDegreeEnd()`.
1477: Level: developer
1479: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`
1480: */
1481: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
1482: {
1483: PetscInt pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
1484: PetscInt *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
1485: PetscSFNode *leafLocationsNew;
1486: const PetscSFNode *iremote;
1487: const PetscInt *ilocal;
1488: PetscBool *isLeaf;
1489: PetscSF sf;
1490: MPI_Comm comm;
1491: PetscMPIInt rank, size;
1493: PetscFunctionBegin;
1494: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1495: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1496: PetscCallMPI(MPI_Comm_size(comm, &size));
1497: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1499: PetscCall(DMGetPointSF(dm, &sf));
1500: PetscCall(PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote));
1501: PetscCall(PetscMalloc1(pEnd - pStart, &isLeaf));
1502: for (i = 0; i < pEnd - pStart; i++) isLeaf[i] = PETSC_FALSE;
1503: for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;
1505: PetscCall(PetscMalloc1(pEnd - pStart + 1, &cumSumDegrees));
1506: cumSumDegrees[0] = 0;
1507: for (i = 1; i <= pEnd - pStart; i++) cumSumDegrees[i] = cumSumDegrees[i - 1] + degrees[i - 1];
1508: sumDegrees = cumSumDegrees[pEnd - pStart];
1509: /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */
1511: PetscCall(PetscMalloc1(sumDegrees, &locationsOfLeafs));
1512: PetscCall(PetscMalloc1(pEnd - pStart, &rankOnLeafs));
1513: for (i = 0; i < pEnd - pStart; i++) rankOnLeafs[i] = rank;
1514: PetscCall(PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs));
1515: PetscCall(PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs));
1516: PetscCall(PetscFree(rankOnLeafs));
1518: /* get the remote local points of my leaves */
1519: PetscCall(PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs));
1520: PetscCall(PetscMalloc1(pEnd - pStart, &points));
1521: for (i = 0; i < pEnd - pStart; i++) points[i] = pStart + i;
1522: PetscCall(PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs));
1523: PetscCall(PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs));
1524: PetscCall(PetscFree(points));
1525: /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
1526: PetscCall(PetscMalloc1(pEnd - pStart, &newOwners));
1527: PetscCall(PetscMalloc1(pEnd - pStart, &newNumbers));
1528: for (i = 0; i < pEnd - pStart; i++) {
1529: newOwners[i] = -1;
1530: newNumbers[i] = -1;
1531: }
1532: {
1533: PetscInt oldNumber, newNumber, oldOwner, newOwner;
1534: for (i = 0; i < n; i++) {
1535: oldNumber = pointsToRewrite[i];
1536: newNumber = -1;
1537: oldOwner = rank;
1538: newOwner = targetOwners[i];
1539: if (oldOwner == newOwner) {
1540: newNumber = oldNumber;
1541: } else {
1542: for (j = 0; j < degrees[oldNumber]; j++) {
1543: if (locationsOfLeafs[cumSumDegrees[oldNumber] + j] == newOwner) {
1544: newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber] + j];
1545: break;
1546: }
1547: }
1548: }
1549: PetscCheck(newNumber != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Couldn't find the new owner of vertex.");
1551: newOwners[oldNumber] = newOwner;
1552: newNumbers[oldNumber] = newNumber;
1553: }
1554: }
1555: PetscCall(PetscFree(cumSumDegrees));
1556: PetscCall(PetscFree(locationsOfLeafs));
1557: PetscCall(PetscFree(remoteLocalPointOfLeafs));
1559: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE));
1560: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE));
1561: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE));
1562: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE));
1564: /* Now count how many leafs we have on each processor. */
1565: leafCounter = 0;
1566: for (i = 0; i < pEnd - pStart; i++) {
1567: if (newOwners[i] >= 0) {
1568: if (newOwners[i] != rank) leafCounter++;
1569: } else {
1570: if (isLeaf[i]) leafCounter++;
1571: }
1572: }
1574: /* Now set up the new sf by creating the leaf arrays */
1575: PetscCall(PetscMalloc1(leafCounter, &leafsNew));
1576: PetscCall(PetscMalloc1(leafCounter, &leafLocationsNew));
1578: leafCounter = 0;
1579: counter = 0;
1580: for (i = 0; i < pEnd - pStart; i++) {
1581: if (newOwners[i] >= 0) {
1582: if (newOwners[i] != rank) {
1583: leafsNew[leafCounter] = i;
1584: leafLocationsNew[leafCounter].rank = newOwners[i];
1585: leafLocationsNew[leafCounter].index = newNumbers[i];
1586: leafCounter++;
1587: }
1588: } else {
1589: if (isLeaf[i]) {
1590: leafsNew[leafCounter] = i;
1591: leafLocationsNew[leafCounter].rank = iremote[counter].rank;
1592: leafLocationsNew[leafCounter].index = iremote[counter].index;
1593: leafCounter++;
1594: }
1595: }
1596: if (isLeaf[i]) counter++;
1597: }
1599: PetscCall(PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER));
1600: PetscCall(PetscFree(newOwners));
1601: PetscCall(PetscFree(newNumbers));
1602: PetscCall(PetscFree(isLeaf));
1603: PetscFunctionReturn(PETSC_SUCCESS);
1604: }
1606: /* The function below is used by DMPlexRebalanceSharedPoints which errors
1607: * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
1608: * this out in that case. */
1609: #if PetscDefined(HAVE_PARMETIS)
1610: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
1611: {
1612: PetscInt *distribution, min, max, sum;
1613: PetscMPIInt rank, size;
1615: PetscFunctionBegin;
1616: PetscCallMPI(MPI_Comm_size(comm, &size));
1617: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1618: PetscCall(PetscCalloc1(size, &distribution));
1619: for (PetscInt i = 0; i < n; i++) {
1620: if (part) distribution[part[i]] += vtxwgt[skip * i];
1621: else distribution[rank] += vtxwgt[skip * i];
1622: }
1623: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm));
1624: min = distribution[0];
1625: max = distribution[0];
1626: sum = distribution[0];
1627: for (PetscInt i = 1; i < size; i++) {
1628: if (distribution[i] < min) min = distribution[i];
1629: if (distribution[i] > max) max = distribution[i];
1630: sum += distribution[i];
1631: }
1632: PetscCall(PetscViewerASCIIPrintf(viewer, "Min: %" PetscInt_FMT ", Avg: %" PetscInt_FMT ", Max: %" PetscInt_FMT ", Balance: %f\n", min, sum / size, max, (max * 1. * size) / sum));
1633: PetscCall(PetscFree(distribution));
1634: PetscFunctionReturn(PETSC_SUCCESS);
1635: }
1636: #endif
1638: /*@
1639: DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the `PointSF` of the `DM` inplace.
1641: Input Parameters:
1642: + dm - The `DMPLEX` object.
1643: . entityDepth - depth of the entity to balance (0 -> balance vertices).
1644: . useInitialGuess - whether to use the current distribution as initial guess (only used by ParMETIS).
1645: - parallel - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.
1647: Output Parameter:
1648: . success - whether the graph partitioning was successful or not, optional. Unsuccessful simply means no change to the partitioning
1650: Options Database Keys:
1651: + -dm_plex_rebalance_shared_points_parmetis - Use ParMETIS instead of METIS for the partitioner
1652: . -dm_plex_rebalance_shared_points_use_initial_guess - Use current partition to bootstrap ParMETIS partition
1653: . -dm_plex_rebalance_shared_points_use_mat_partitioning - Use the MatPartitioning object to perform the partition, the prefix for those operations is -dm_plex_rebalance_shared_points_
1654: - -dm_plex_rebalance_shared_points_monitor - Monitor the shared points rebalance process
1656: Level: intermediate
1658: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexDistribute()`
1659: @*/
1660: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
1661: {
1662: #if PetscDefined(HAVE_PARMETIS)
1663: PetscSF sf;
1664: PetscInt i, j, idx, jdx;
1665: PetscInt eBegin, eEnd, nroots, nleafs, pStart, pEnd;
1666: const PetscInt *degrees, *ilocal;
1667: const PetscSFNode *iremote;
1668: PetscBool *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
1669: PetscInt numExclusivelyOwned, numNonExclusivelyOwned;
1670: PetscMPIInt rank, size;
1671: PetscInt *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
1672: const PetscInt *cumSumVertices;
1673: PetscInt offset, counter;
1674: PetscInt *vtxwgt;
1675: const PetscInt *xadj, *adjncy;
1676: PetscInt *part, *options;
1677: PetscInt nparts, wgtflag, numflag, ncon, edgecut;
1678: real_t *ubvec;
1679: PetscInt *firstVertices, *renumbering;
1680: PetscInt failed;
1681: MPI_Comm comm;
1682: Mat A;
1683: PetscViewer viewer;
1684: PetscViewerFormat format;
1685: PetscLayout layout;
1686: real_t *tpwgts;
1687: PetscMPIInt *counts, *mpiCumSumVertices;
1688: PetscInt *pointsToRewrite;
1689: PetscInt numRows;
1690: PetscBool done, usematpartitioning = PETSC_FALSE;
1691: IS ispart = NULL;
1692: MatPartitioning mp;
1693: const char *prefix;
1695: PetscFunctionBegin;
1696: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1697: PetscCallMPI(MPI_Comm_size(comm, &size));
1698: if (size == 1) {
1699: if (success) *success = PETSC_TRUE;
1700: PetscFunctionReturn(PETSC_SUCCESS);
1701: }
1702: if (success) *success = PETSC_FALSE;
1703: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1705: parallel = PETSC_FALSE;
1706: useInitialGuess = PETSC_FALSE;
1707: PetscObjectOptionsBegin((PetscObject)dm);
1708: PetscCall(PetscOptionsName("-dm_plex_rebalance_shared_points_parmetis", "Use ParMETIS instead of METIS for the partitioner", "DMPlexRebalanceSharedPoints", ¶llel));
1709: PetscCall(PetscOptionsBool("-dm_plex_rebalance_shared_points_use_initial_guess", "Use current partition to bootstrap ParMETIS partition", "DMPlexRebalanceSharedPoints", useInitialGuess, &useInitialGuess, NULL));
1710: PetscCall(PetscOptionsBool("-dm_plex_rebalance_shared_points_use_mat_partitioning", "Use the MatPartitioning object to partition", "DMPlexRebalanceSharedPoints", usematpartitioning, &usematpartitioning, NULL));
1711: PetscCall(PetscOptionsViewer("-dm_plex_rebalance_shared_points_monitor", "Monitor the shared points rebalance process", "DMPlexRebalanceSharedPoints", &viewer, &format, NULL));
1712: PetscOptionsEnd();
1713: if (viewer) PetscCall(PetscViewerPushFormat(viewer, format));
1715: PetscCall(PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0));
1717: PetscCall(DMGetOptionsPrefix(dm, &prefix));
1718: PetscCall(PetscOptionsCreateViewer(comm, ((PetscObject)dm)->options, prefix, "-dm_rebalance_partition_view", &viewer, &format, NULL));
1719: if (viewer) PetscCall(PetscViewerPushFormat(viewer, format));
1721: /* Figure out all points in the plex that we are interested in balancing. */
1722: PetscCall(DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd));
1723: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1724: PetscCall(PetscMalloc1(pEnd - pStart, &toBalance));
1725: for (i = 0; i < pEnd - pStart; i++) toBalance[i] = (PetscBool)(i >= eBegin && i < eEnd);
1727: /* There are three types of points:
1728: * exclusivelyOwned: points that are owned by this process and only seen by this process
1729: * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
1730: * leaf: a point that is seen by this process but owned by a different process
1731: */
1732: PetscCall(DMGetPointSF(dm, &sf));
1733: PetscCall(PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote));
1734: PetscCall(PetscMalloc1(pEnd - pStart, &isLeaf));
1735: PetscCall(PetscMalloc1(pEnd - pStart, &isNonExclusivelyOwned));
1736: PetscCall(PetscMalloc1(pEnd - pStart, &isExclusivelyOwned));
1737: for (i = 0; i < pEnd - pStart; i++) {
1738: isNonExclusivelyOwned[i] = PETSC_FALSE;
1739: isExclusivelyOwned[i] = PETSC_FALSE;
1740: isLeaf[i] = PETSC_FALSE;
1741: }
1743: /* mark all the leafs */
1744: for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;
1746: /* for an owned point, we can figure out whether another processor sees it or
1747: * not by calculating its degree */
1748: PetscCall(PetscSFComputeDegreeBegin(sf, °rees));
1749: PetscCall(PetscSFComputeDegreeEnd(sf, °rees));
1750: numExclusivelyOwned = 0;
1751: numNonExclusivelyOwned = 0;
1752: for (i = 0; i < pEnd - pStart; i++) {
1753: if (toBalance[i]) {
1754: if (degrees[i] > 0) {
1755: isNonExclusivelyOwned[i] = PETSC_TRUE;
1756: numNonExclusivelyOwned += 1;
1757: } else {
1758: if (!isLeaf[i]) {
1759: isExclusivelyOwned[i] = PETSC_TRUE;
1760: numExclusivelyOwned += 1;
1761: }
1762: }
1763: }
1764: }
1766: /* Build a graph with one vertex per core representing the
1767: * exclusively owned points and then one vertex per nonExclusively owned
1768: * point. */
1769: PetscCall(PetscLayoutCreate(comm, &layout));
1770: PetscCall(PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned));
1771: PetscCall(PetscLayoutSetUp(layout));
1772: PetscCall(PetscLayoutGetRanges(layout, &cumSumVertices));
1773: PetscCall(PetscMalloc1(pEnd - pStart, &globalNumbersOfLocalOwnedVertices));
1774: for (i = 0; i < pEnd - pStart; i++) globalNumbersOfLocalOwnedVertices[i] = pStart - 1;
1775: offset = cumSumVertices[rank];
1776: counter = 0;
1777: for (i = 0; i < pEnd - pStart; i++) {
1778: if (toBalance[i]) {
1779: if (degrees[i] > 0) {
1780: globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
1781: counter++;
1782: }
1783: }
1784: }
1786: /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
1787: PetscCall(PetscMalloc1(pEnd - pStart, &leafGlobalNumbers));
1788: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE));
1789: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE));
1791: /* Build the graph for partitioning */
1792: numRows = 1 + numNonExclusivelyOwned;
1793: PetscCall(PetscLogEventBegin(DMPLEX_RebalBuildGraph, dm, 0, 0, 0));
1794: PetscCall(MatCreate(comm, &A));
1795: PetscCall(MatSetType(A, MATMPIADJ));
1796: PetscCall(MatSetSizes(A, numRows, numRows, cumSumVertices[size], cumSumVertices[size]));
1797: idx = cumSumVertices[rank];
1798: for (i = 0; i < pEnd - pStart; i++) {
1799: if (toBalance[i]) {
1800: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1801: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1802: else continue;
1803: PetscCall(MatSetValue(A, idx, jdx, 1, INSERT_VALUES));
1804: PetscCall(MatSetValue(A, jdx, idx, 1, INSERT_VALUES));
1805: }
1806: }
1807: PetscCall(PetscFree(globalNumbersOfLocalOwnedVertices));
1808: PetscCall(PetscFree(leafGlobalNumbers));
1809: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1810: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1811: PetscCall(PetscLogEventEnd(DMPLEX_RebalBuildGraph, dm, 0, 0, 0));
1813: nparts = size;
1814: ncon = 1;
1815: PetscCall(PetscMalloc1(ncon * nparts, &tpwgts));
1816: for (i = 0; i < ncon * nparts; i++) tpwgts[i] = (real_t)(1. / (nparts));
1817: PetscCall(PetscMalloc1(ncon, &ubvec));
1818: for (i = 0; i < ncon; i++) ubvec[i] = (real_t)1.05;
1820: PetscCall(PetscMalloc1(ncon * (1 + numNonExclusivelyOwned), &vtxwgt));
1821: if (ncon == 2) {
1822: vtxwgt[0] = numExclusivelyOwned;
1823: vtxwgt[1] = 1;
1824: for (i = 0; i < numNonExclusivelyOwned; i++) {
1825: vtxwgt[ncon * (i + 1)] = 1;
1826: vtxwgt[ncon * (i + 1) + 1] = 0;
1827: }
1828: } else {
1829: PetscInt base, ms;
1830: PetscCallMPI(MPIU_Allreduce(&numExclusivelyOwned, &base, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
1831: PetscCall(MatGetSize(A, &ms, NULL));
1832: ms -= size;
1833: base = PetscMax(base, ms);
1834: vtxwgt[0] = base + numExclusivelyOwned;
1835: for (i = 0; i < numNonExclusivelyOwned; i++) vtxwgt[i + 1] = 1;
1836: }
1838: if (viewer) {
1839: PetscCall(PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %" PetscInt_FMT " on interface of mesh distribution.\n", entityDepth));
1840: PetscCall(PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %" PetscInt_FMT "\n", cumSumVertices[size]));
1841: }
1842: /* TODO: Drop the parallel/sequential choice here and just use MatPartioner for much more flexibility */
1843: if (usematpartitioning) {
1844: const char *prefix;
1846: PetscCall(MatPartitioningCreate(PetscObjectComm((PetscObject)dm), &mp));
1847: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mp, "dm_plex_rebalance_shared_points_"));
1848: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
1849: PetscCall(PetscObjectPrependOptionsPrefix((PetscObject)mp, prefix));
1850: PetscCall(MatPartitioningSetAdjacency(mp, A));
1851: PetscCall(MatPartitioningSetNumberVertexWeights(mp, ncon));
1852: PetscCall(MatPartitioningSetVertexWeights(mp, vtxwgt));
1853: PetscCall(MatPartitioningSetFromOptions(mp));
1854: PetscCall(MatPartitioningApply(mp, &ispart));
1855: PetscCall(ISGetIndices(ispart, (const PetscInt **)&part));
1856: } else if (parallel) {
1857: if (viewer) PetscCall(PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n"));
1858: PetscCall(PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part));
1859: PetscCall(MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done));
1860: PetscCheck(done, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not get adjacency information");
1861: PetscCall(PetscMalloc1(4, &options));
1862: options[0] = 1;
1863: options[1] = 0; /* Verbosity */
1864: if (viewer) options[1] = 1;
1865: options[2] = 0; /* Seed */
1866: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
1867: wgtflag = 2;
1868: numflag = 0;
1869: if (useInitialGuess) {
1870: PetscCall(PetscViewerASCIIPrintf(viewer, "THIS DOES NOT WORK! I don't know why. Using current distribution of points as initial guess.\n"));
1871: for (i = 0; i < numRows; i++) part[i] = rank;
1872: if (viewer) PetscCall(PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n"));
1873: PetscCall(PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0));
1874: PetscCallParMETIS(ParMETIS_V3_RefineKway, (PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1875: PetscCall(PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0));
1876: } else {
1877: PetscCall(PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0));
1878: PetscCallParMETIS(ParMETIS_V3_PartKway, (PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1879: PetscCall(PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0));
1880: }
1881: PetscCall(MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done));
1882: PetscCall(PetscFree(options));
1883: } else {
1884: if (viewer) PetscCall(PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n"));
1885: Mat As;
1886: PetscInt *partGlobal;
1887: PetscInt *numExclusivelyOwnedAll;
1889: PetscCall(PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part));
1890: PetscCall(MatGetSize(A, &numRows, NULL));
1891: PetscCall(PetscLogEventBegin(DMPLEX_RebalGatherGraph, dm, 0, 0, 0));
1892: PetscCall(MatMPIAdjToSeqRankZero(A, &As));
1893: PetscCall(PetscLogEventEnd(DMPLEX_RebalGatherGraph, dm, 0, 0, 0));
1895: PetscCall(PetscMalloc1(size, &numExclusivelyOwnedAll));
1896: numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
1897: PetscCallMPI(MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, numExclusivelyOwnedAll, 1, MPIU_INT, comm));
1899: PetscCall(PetscMalloc1(numRows, &partGlobal));
1900: PetscCall(PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0));
1901: if (rank == 0) {
1902: PetscInt *vtxwgt_g, numRows_g;
1904: PetscCall(MatGetRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done));
1905: PetscCall(PetscMalloc1(2 * numRows_g, &vtxwgt_g));
1906: for (i = 0; i < size; i++) {
1907: vtxwgt_g[ncon * cumSumVertices[i]] = numExclusivelyOwnedAll[i];
1908: if (ncon > 1) vtxwgt_g[ncon * cumSumVertices[i] + 1] = 1;
1909: for (j = cumSumVertices[i] + 1; j < cumSumVertices[i + 1]; j++) {
1910: vtxwgt_g[ncon * j] = 1;
1911: if (ncon > 1) vtxwgt_g[2 * j + 1] = 0;
1912: }
1913: }
1915: PetscCall(PetscMalloc1(64, &options));
1916: PetscCallMETIS(METIS_SetDefaultOptions, options);
1917: options[METIS_OPTION_CONTIG] = 1;
1918: PetscCallMETIS(METIS_PartGraphKway, &numRows_g, &ncon, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
1919: PetscCall(PetscFree(options));
1920: PetscCall(PetscFree(vtxwgt_g));
1921: PetscCall(MatRestoreRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done));
1922: PetscCall(MatDestroy(&As));
1923: }
1924: PetscCall(PetscBarrier((PetscObject)dm));
1925: PetscCall(PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0));
1926: PetscCall(PetscFree(numExclusivelyOwnedAll));
1928: /* scatter the partitioning information to ranks */
1929: PetscCall(PetscLogEventBegin(DMPLEX_RebalScatterPart, 0, 0, 0, 0));
1930: PetscCall(PetscMalloc1(size, &counts));
1931: PetscCall(PetscMalloc1(size + 1, &mpiCumSumVertices));
1932: for (i = 0; i < size; i++) PetscCall(PetscMPIIntCast(cumSumVertices[i + 1] - cumSumVertices[i], &counts[i]));
1933: for (i = 0; i <= size; i++) PetscCall(PetscMPIIntCast(cumSumVertices[i], &mpiCumSumVertices[i]));
1934: PetscCallMPI(MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm));
1935: PetscCall(PetscFree(counts));
1936: PetscCall(PetscFree(mpiCumSumVertices));
1937: PetscCall(PetscFree(partGlobal));
1938: PetscCall(PetscLogEventEnd(DMPLEX_RebalScatterPart, 0, 0, 0, 0));
1939: }
1940: PetscCall(PetscFree(ubvec));
1941: PetscCall(PetscFree(tpwgts));
1943: /* Rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
1944: PetscCall(PetscMalloc2(size, &firstVertices, size, &renumbering));
1945: firstVertices[rank] = part[0];
1946: PetscCallMPI(MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, firstVertices, 1, MPIU_INT, comm));
1947: for (i = 0; i < size; i++) renumbering[firstVertices[i]] = i;
1948: for (i = 0; i < cumSumVertices[rank + 1] - cumSumVertices[rank]; i++) part[i] = renumbering[part[i]];
1949: PetscCall(PetscFree2(firstVertices, renumbering));
1951: /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
1952: failed = (PetscInt)(part[0] != rank);
1953: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &failed, 1, MPIU_INT, MPI_SUM, comm));
1954: if (failed > 0) {
1955: PetscCheck(failed <= 0, comm, PETSC_ERR_LIB, "METIS/ParMETIS returned a bad partition");
1956: PetscCall(PetscFree(vtxwgt));
1957: PetscCall(PetscFree(toBalance));
1958: PetscCall(PetscFree(isLeaf));
1959: PetscCall(PetscFree(isNonExclusivelyOwned));
1960: PetscCall(PetscFree(isExclusivelyOwned));
1961: if (usematpartitioning) {
1962: PetscCall(ISRestoreIndices(ispart, (const PetscInt **)&part));
1963: PetscCall(ISDestroy(&ispart));
1964: } else PetscCall(PetscFree(part));
1965: if (viewer) {
1966: PetscCall(PetscViewerPopFormat(viewer));
1967: PetscCall(PetscViewerDestroy(&viewer));
1968: }
1969: PetscCall(PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0));
1970: PetscFunctionReturn(PETSC_SUCCESS);
1971: }
1973: /* Check how well we did distributing points*/
1974: if (viewer) {
1975: PetscCall(PetscViewerASCIIPrintf(viewer, "Number of owned entities of depth %" PetscInt_FMT ".\n", entityDepth));
1976: PetscCall(PetscViewerASCIIPrintf(viewer, "Initial "));
1977: PetscCall(DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, NULL, viewer));
1978: PetscCall(PetscViewerASCIIPrintf(viewer, "Rebalanced "));
1979: PetscCall(DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer));
1980: }
1982: /* Check that every vertex is owned by a process that it is actually connected to. */
1983: PetscCall(MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
1984: for (i = 1; i <= numNonExclusivelyOwned; i++) {
1985: PetscInt loc = 0;
1986: PetscCall(PetscFindInt(cumSumVertices[part[i]], xadj[i + 1] - xadj[i], &adjncy[xadj[i]], &loc));
1987: /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
1988: if (loc < 0) part[i] = rank;
1989: }
1990: PetscCall(MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done));
1991: PetscCall(MatDestroy(&A));
1993: /* See how significant the influences of the previous fixing up step was.*/
1994: if (viewer) {
1995: PetscCall(PetscViewerASCIIPrintf(viewer, "After fix "));
1996: PetscCall(DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer));
1997: }
1998: if (!usematpartitioning) PetscCall(PetscFree(vtxwgt));
1999: else PetscCall(MatPartitioningDestroy(&mp));
2001: PetscCall(PetscLayoutDestroy(&layout));
2003: PetscCall(PetscLogEventBegin(DMPLEX_RebalRewriteSF, dm, 0, 0, 0));
2004: /* Rewrite the SF to reflect the new ownership. */
2005: PetscCall(PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite));
2006: counter = 0;
2007: for (i = 0; i < pEnd - pStart; i++) {
2008: if (toBalance[i]) {
2009: if (isNonExclusivelyOwned[i]) {
2010: pointsToRewrite[counter] = i + pStart;
2011: counter++;
2012: }
2013: }
2014: }
2015: PetscCall(DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part + 1, degrees));
2016: PetscCall(PetscFree(pointsToRewrite));
2017: PetscCall(PetscLogEventEnd(DMPLEX_RebalRewriteSF, dm, 0, 0, 0));
2019: PetscCall(PetscFree(toBalance));
2020: PetscCall(PetscFree(isLeaf));
2021: PetscCall(PetscFree(isNonExclusivelyOwned));
2022: PetscCall(PetscFree(isExclusivelyOwned));
2023: if (usematpartitioning) {
2024: PetscCall(ISRestoreIndices(ispart, (const PetscInt **)&part));
2025: PetscCall(ISDestroy(&ispart));
2026: } else PetscCall(PetscFree(part));
2027: if (viewer) {
2028: PetscCall(PetscViewerPopFormat(viewer));
2029: PetscCall(PetscViewerDestroy(&viewer));
2030: }
2031: if (success) *success = PETSC_TRUE;
2032: PetscCall(PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0));
2033: PetscFunctionReturn(PETSC_SUCCESS);
2034: #else
2035: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2036: #endif
2037: }
2039: // If the point is in the closure of a label cell, set the owner to this process
2040: static PetscErrorCode CheckLabelPoint_Private(DM plex, DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt point, PetscInt *owner)
2041: {
2042: PetscInt starSize, *star = NULL;
2043: PetscMPIInt rank;
2045: PetscFunctionBegin;
2046: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)plex), &rank));
2047: PetscCall(DMPlexGetTransitiveClosure(plex, point, PETSC_FALSE, &starSize, &star));
2048: for (PetscInt s = 0; s < starSize; ++s) {
2049: const PetscInt spoint = star[s * 2];
2050: PetscInt depth, val;
2052: PetscCall(DMPlexGetPointDepth(plex, spoint, &depth));
2053: PetscCall(DMLabelGetValue(label, spoint, &val));
2054: for (PetscInt v = 0; v < Nv; ++v) {
2055: if (val == values[v]) {
2056: *owner = rank;
2057: s = starSize;
2058: break;
2059: }
2060: }
2061: }
2062: PetscCall(DMPlexRestoreTransitiveClosure(plex, point, PETSC_FALSE, &starSize, &star));
2063: PetscFunctionReturn(PETSC_SUCCESS);
2064: }
2066: /*@
2067: DMPlexRebalanceSharedLabelPoints - Change ownership of labeled points in the Plex so that processes owning shared label points also own a cell that contains them. This routine updates the `PointSF` of the `DM` inplace.
2069: Input Parameters:
2070: + dm - the `DMPLEX` object
2071: . label - the `DMLabel` object
2072: . Nv - the number of label values
2073: . values - the array of label values
2074: - cellDepth - depth of the cells in the label
2076: Level: intermediate
2078: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexDistribute()`, `DMPlexRebalanceSharedPoints()`
2079: @*/
2080: PetscErrorCode DMPlexRebalanceSharedLabelPoints(DM dm, DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt cellDepth)
2081: {
2082: DM plex;
2083: PetscSF sf;
2084: const PetscInt *degrees, *leaves;
2085: PetscInt *lowner, *gowner, *newPoints, *newOwner;
2086: PetscInt Nr, Nl, numNewOwners = 0, tmp = 0;
2087: PetscMPIInt rank, size;
2088: MPI_Comm comm;
2089: PetscInt debug = ((DM_Plex *)dm->data)->printCohesive;
2091: PetscFunctionBegin;
2092: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2093: PetscCallMPI(MPI_Comm_size(comm, &size));
2094: if (size == 1) {
2095: PetscFunctionReturn(PETSC_SUCCESS);
2096: }
2097: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2098: PetscCall(DMConvert(dm, DMPLEX, &plex));
2099: PetscCall(DMGetPointSF(dm, &sf));
2100: PetscCall(PetscSFGetGraph(sf, &Nr, &Nl, &leaves, NULL));
2101: PetscCall(PetscSFComputeDegreeBegin(sf, °rees));
2102: PetscCall(PetscSFComputeDegreeEnd(sf, °rees));
2103: PetscCall(PetscMalloc2(Nr, &lowner, Nr, &gowner));
2104: // Loop over all shared points
2105: // If a shared point is in the label and connected to a label cell, then vote for it
2106: // -2: Not a labeled point
2107: // -1: Labeled point that cannot be owned by this process
2108: // rank: Labeled point that could be owned by this process
2109: for (PetscInt l = 0; l < Nl; ++l) {
2110: const PetscInt leaf = leaves ? leaves[l] : l;
2111: PetscInt val;
2113: lowner[leaf] = -2;
2114: PetscCall(DMLabelGetValue(label, leaf, &val));
2115: for (PetscInt v = 0; v < Nv; ++v) {
2116: if (val == values[v]) {
2117: lowner[leaf] = -1;
2118: PetscCall(CheckLabelPoint_Private(plex, label, Nv, values, leaf, &lowner[leaf]));
2119: if (debug > 3) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Marked leaf point %" PetscInt_FMT " as %" PetscInt_FMT "\n", rank, leaf, lowner[leaf]));
2120: break;
2121: }
2122: }
2123: }
2124: for (PetscInt root = 0; root < Nr; ++root) {
2125: gowner[root] = -2;
2126: if (degrees[root] > 0) {
2127: PetscInt val;
2129: PetscCall(DMLabelGetValue(label, root, &val));
2130: for (PetscInt v = 0; v < Nv; ++v) {
2131: if (val == values[v]) {
2132: gowner[root] = -1;
2133: PetscCall(CheckLabelPoint_Private(plex, label, Nv, values, root, &gowner[root]));
2134: if (debug > 3) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Marked root point %" PetscInt_FMT " as %" PetscInt_FMT "\n", rank, root, gowner[root]));
2135: break;
2136: }
2137: }
2138: }
2139: }
2140: // Process votes
2141: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, lowner, gowner, MPI_MAX));
2142: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, lowner, gowner, MPI_MAX));
2143: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, gowner, lowner, MPI_MAX));
2144: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, gowner, lowner, MPI_MAX));
2145: // Figure out which points changed ownership and rewrite the SF
2146: for (PetscInt l = 0; l < Nl; ++l) {
2147: const PetscInt leaf = leaves ? leaves[l] : l;
2149: if (lowner[leaf] == rank) ++numNewOwners;
2150: }
2151: for (PetscInt root = 0; root < Nr; ++root) {
2152: PetscCheck(!degrees[root] || gowner[root] != -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "[%d]Shared point %" PetscInt_FMT " was not assigned an owner", rank, root);
2153: if (degrees[root] > 0 && gowner[root] != -2 && gowner[root] != rank) ++numNewOwners;
2154: }
2155: PetscCall(PetscMalloc2(numNewOwners, &newPoints, numNewOwners, &newOwner));
2156: for (PetscInt l = 0; l < Nl; ++l) {
2157: const PetscInt leaf = leaves ? leaves[l] : l;
2159: if (lowner[leaf] == rank) {
2160: newPoints[tmp] = leaf;
2161: newOwner[tmp++] = rank;
2162: if (debug > 3) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Changed leaf point %" PetscInt_FMT " to owner %d\n", rank, leaf, rank));
2163: }
2164: }
2165: for (PetscInt root = 0; root < Nr; ++root) {
2166: if (degrees[root] > 0 && gowner[root] != -2 && gowner[root] != rank) {
2167: newPoints[tmp] = root;
2168: newOwner[tmp++] = gowner[root];
2169: if (debug > 3) PetscCall(PetscSynchronizedPrintf(comm, "[%d] Changed root point %" PetscInt_FMT " to owner %" PetscInt_FMT "\n", rank, root, gowner[root]));
2170: }
2171: }
2172: PetscCheck(tmp == numNewOwners, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of new owners %" PetscInt_FMT " != %" PetscInt_FMT, tmp, numNewOwners);
2173: if (debug > 3) {
2174: PetscCall(PetscSynchronizedPrintf(comm, "[%d] Rewrote %" PetscInt_FMT " points\n", rank, numNewOwners));
2175: for (PetscInt n = 0; n < numNewOwners; ++n) {
2176: PetscCall(PetscSynchronizedPrintf(comm, "[%d] point %" PetscInt_FMT " now owned by %" PetscInt_FMT "\n", rank, newPoints[n], newOwner[n]));
2177: }
2178: PetscCall(PetscSynchronizedFlush(comm, NULL));
2179: }
2180: PetscCall(DMPlexRewriteSF(dm, numNewOwners, newPoints, newOwner, degrees));
2181: PetscCall(PetscFree2(newPoints, newOwner));
2182: PetscCall(PetscFree2(lowner, gowner));
2183: PetscCall(DMDestroy(&plex));
2184: PetscFunctionReturn(PETSC_SUCCESS);
2185: }