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