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