42 # include <parmetis.h>
63 Foam::label Foam::parMetisDecomp::decompose
68 Field<int>& cellWeights,
69 Field<int>& faceWeights,
70 const List<int>& options,
72 List<int>& finalDecomp
82 if (cellCentres.size() != xadj.size()-1)
85 <<
"cellCentres:" << cellCentres.size()
86 <<
" xadj:" << xadj.size()
92 List<int> nLocalCells(Pstream::nProcs());
93 nLocalCells[Pstream::myProcNo()] = xadj.size()-1;
94 Pstream::gatherList(nLocalCells);
95 Pstream::scatterList(nLocalCells);
98 List<int> cellOffsets(Pstream::nProcs()+1);
100 forAll(nLocalCells, procI)
102 cellOffsets[procI] = nGlobalCells;
103 nGlobalCells += nLocalCells[procI];
105 cellOffsets[Pstream::nProcs()] = nGlobalCells;
108 Field<floatScalar> xyz(3*cellCentres.size());
110 forAll(cellCentres, cellI)
112 const point& cc = cellCentres[cellI];
113 xyz[compI++] = float(cc.x());
114 xyz[compI++] = float(cc.y());
115 xyz[compI++] = float(cc.z());
127 List<int> nSendCells(Pstream::nProcs(), 0);
129 for (label procI = nLocalCells.size()-1; procI >=1; procI--)
131 if (nLocalCells[procI]-nSendCells[procI] < 1)
133 nSendCells[procI-1] = nSendCells[procI]-nLocalCells[procI]+1;
139 if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
142 IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
144 Field<int> prevXadj(fromPrevProc);
145 Field<int> prevAdjncy(fromPrevProc);
146 Field<floatScalar> prevXyz(fromPrevProc);
147 Field<int> prevCellWeights(fromPrevProc);
148 Field<int> prevFaceWeights(fromPrevProc);
150 if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1])
153 <<
"Expected from processor " << Pstream::myProcNo()-1
154 <<
" connectivity for " << nSendCells[Pstream::myProcNo()-1]
155 <<
" nCells but only received " << prevXadj.size()
160 prepend(prevAdjncy, adjncy);
162 xadj += prevAdjncy.size();
163 prepend(prevXadj, xadj);
165 prepend(prevXyz, xyz);
167 prepend(prevCellWeights, cellWeights);
168 prepend(prevFaceWeights, faceWeights);
174 if (nSendCells[Pstream::myProcNo()] > 0)
177 OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1);
179 int nCells = nSendCells[Pstream::myProcNo()];
180 int startCell = xadj.size()-1 - nCells;
181 int startFace = xadj[startCell];
182 int nFaces = adjncy.size()-startFace;
187 << Field<int>::subField(xadj, nCells, startCell)-startFace
188 << Field<int>::subField(adjncy, nFaces, startFace)
189 << SubField<floatScalar>(xyz, nDims*nCells, nDims*startCell)
193 ?
static_cast<const Field<int>&
>
195 Field<int>::subField(cellWeights, nCells, startCell)
202 ?
static_cast<const Field<int>&
>
204 Field<int>::subField(faceWeights, nFaces, startFace)
210 if (faceWeights.size())
212 faceWeights.setSize(faceWeights.size()-nFaces);
214 if (cellWeights.size())
216 cellWeights.setSize(cellWeights.size()-nCells);
218 xyz.setSize(xyz.size()-nDims*nCells);
219 adjncy.setSize(adjncy.size()-nFaces);
220 xadj.setSize(xadj.size() - nCells);
229 nLocalCells[procI] -= nSendCells[procI];
234 nLocalCells[procI] += nSendCells[procI-1];
239 forAll(nLocalCells, procI)
241 cellOffsets[procI] = nGlobalCells;
242 nGlobalCells += nLocalCells[procI];
246 if (nLocalCells[Pstream::myProcNo()] != (xadj.size()-1))
249 <<
"Have connectivity for " << xadj.size()-1
250 <<
" cells but nLocalCells:" << nLocalCells[Pstream::myProcNo()]
257 int* adjwgtPtr = NULL;
259 if (cellWeights.size())
261 vwgtPtr = cellWeights.begin();
264 if (faceWeights.size())
266 adjwgtPtr = faceWeights.begin();
274 Field<floatScalar> tpwgts(nCon*nProcessors_, 1./nProcessors_);
276 Field<floatScalar> ubvec(nCon, 1.02);
277 if (nProcessors_ == 1)
283 MPI_Comm comm = MPI_COMM_WORLD;
286 finalDecomp.setSize(nLocalCells[Pstream::myProcNo()]);
292 ParMETIS_V3_PartGeomKway
307 const_cast<List<int>&
>(options).begin(),
318 if (nSendCells[Pstream::myProcNo()] > 0)
320 IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1);
322 List<int> nextFinalDecomp(fromNextProc);
324 if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()])
327 <<
"Expected from processor " << Pstream::myProcNo()+1
328 <<
" decomposition for " << nSendCells[Pstream::myProcNo()]
329 <<
" nCells but only received " << nextFinalDecomp.size()
333 append(nextFinalDecomp, finalDecomp);
337 if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
339 OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
341 int nToPrevious = nSendCells[Pstream::myProcNo()-1];
348 finalDecomp.size()-nToPrevious
352 finalDecomp.setSize(finalDecomp.size()-nToPrevious);
361 Foam::parMetisDecomp::parMetisDecomp
363 const dictionary& decompositionDict,
367 decompositionMethod(decompositionDict),
376 "parMetisDecomp::parMetisDecomp()"
377 ) <<
"Failed to load the library 'metisDecomp.so'" <<
endl
392 if (cc.size() != mesh_.nCells())
396 "parMetisDecomp::decompose"
397 "(const pointField&, const scalarField&)"
398 ) <<
"Can use this decomposition method only for the whole mesh"
400 <<
"and supply one coordinate (cellCentre) for every cell." <<
endl
401 <<
"The number of coordinates " << cc.size() <<
endl
402 <<
"The number of cells in the mesh " << mesh_.nCells()
409 return metisDecomp(decompositionDict_, mesh_).decompose
421 calcMetisDistributedCSR
430 List<int> options(3, 0);
436 Field<int> cellWeights;
439 Field<int> faceWeights;
443 scalar minWeights =
gMin(cWeights);
444 if (cWeights.size() > 0)
450 "metisDecomp::decompose"
451 "(const pointField&, const scalarField&)"
452 ) <<
"Illegal minimum weight " << minWeights
456 if (cWeights.size() != mesh_.nCells())
460 "parMetisDecomp::decompose"
461 "(const pointField&, const scalarField&)"
462 ) <<
"Number of cell weights " << cWeights.size()
463 <<
" does not equal number of cells " << mesh_.nCells()
468 cellWeights.setSize(cWeights.size());
471 cellWeights[i] = int(cWeights[i]/minWeights);
477 if (decompositionDict_.found(
"metisCoeffs"))
479 const dictionary& metisCoeffs =
480 decompositionDict_.subDict(
"metisCoeffs");
483 if (metisCoeffs.readIfPresent(
"cellWeightsFile", weightsFile))
485 Info<<
"parMetisDecomp : Using cell-based weights read from "
486 << weightsFile <<
endl;
493 mesh_.time().timeName(),
499 cellWeights.transfer(cellIOWeights);
501 if (cellWeights.size() != mesh_.nCells())
505 "parMetisDecomp::decompose"
506 "(const pointField&, const scalarField&)"
507 ) <<
"Number of cell weights " << cellWeights.size()
508 <<
" read from " << cellIOWeights.objectPath()
509 <<
" does not equal number of cells " << mesh_.nCells()
514 if (metisCoeffs.readIfPresent(
"faceWeightsFile", weightsFile))
516 Info<<
"parMetisDecomp : Using face-based weights read from "
517 << weightsFile <<
endl;
524 mesh_.time().timeName(),
531 if (weights.size() != mesh_.nFaces())
535 "parMetisDecomp::decompose"
536 "(const pointField&, const scalarField&)"
537 ) <<
"Number of face weights " << weights.size()
538 <<
" does not equal number of internal and boundary faces "
543 faceWeights.setSize(adjncy.size());
546 List<int> nFacesPerCell(mesh_.nCells(), 0);
549 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
551 label w = weights[faceI];
553 label own = mesh_.faceOwner()[faceI];
554 label nei = mesh_.faceNeighbour()[faceI];
556 faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
557 faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
560 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
564 const polyPatch& pp = patches[patchI];
568 label faceI = pp.start();
572 label w = weights[faceI];
573 label own = mesh_.faceOwner()[faceI];
574 faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
581 if (metisCoeffs.readIfPresent(
"options", options))
583 Info<<
"Using Metis options " << options
586 if (options.size() != 3)
590 "parMetisDecomp::decompose"
591 "(const pointField&, const scalarField&)"
592 ) <<
"Number of options " << options.size()
600 List<int> finalDecomp;
617 decomp[i] = finalDecomp[i];
630 const labelList& faceOwner = mesh_.faceOwner();
631 const labelList& faceNeighbour = mesh_.faceNeighbour();
632 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
634 if (cellToRegion.size() != mesh_.nCells())
638 "parMetisDecomp::decompose(const labelList&, const pointField&)"
639 ) <<
"Size of cell-to-coarse map " << cellToRegion.size()
640 <<
" differs from number of cells in mesh " << mesh_.nCells()
646 globalIndex globalRegions(regionPoints.size());
652 List<int> globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
656 const polyPatch& pp = patches[patchI];
660 label faceI = pp.start();
661 label bFaceI = pp.start() - mesh_.nInternalFaces();
665 label ownRegion = cellToRegion[faceOwner[faceI]];
666 globalNeighbour[bFaceI++] = globalRegions.toGlobal(ownRegion);
681 List<DynamicList<label> > dynRegionRegions(regionPoints.size());
684 forAll(faceNeighbour, faceI)
686 label ownRegion = cellToRegion[faceOwner[faceI]];
687 label neiRegion = cellToRegion[faceNeighbour[faceI]];
689 if (ownRegion != neiRegion)
691 label globalOwn = globalRegions.toGlobal(ownRegion);
692 label globalNei = globalRegions.toGlobal(neiRegion);
694 if (
findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
696 dynRegionRegions[ownRegion].append(globalNei);
698 if (
findIndex(dynRegionRegions[neiRegion], globalOwn) == -1)
700 dynRegionRegions[neiRegion].append(globalOwn);
708 const polyPatch& pp = patches[patchI];
712 label faceI = pp.start();
713 label bFaceI = pp.start() - mesh_.nInternalFaces();
717 label ownRegion = cellToRegion[faceOwner[faceI]];
718 label globalNei = globalNeighbour[bFaceI++];
721 if (
findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
723 dynRegionRegions[ownRegion].append(globalNei);
729 globalRegionRegions.setSize(dynRegionRegions.size());
730 forAll(dynRegionRegions, i)
732 globalRegionRegions[i].transfer(dynRegionRegions[i]);
747 labelList cellDistribution(cellToRegion.size());
749 forAll(cellDistribution, cellI)
751 cellDistribution[cellI] = regionDecomp[cellToRegion[cellI]];
753 return cellDistribution;
764 if (cellCentres.size() != globalCellCells.size())
768 "parMetisDecomp::decompose(const labelListList&"
769 ", const pointField&, const scalarField&)"
770 ) <<
"Inconsistent number of cells (" << globalCellCells.size()
771 <<
") and number of cell centres (" << cellCentres.size()
772 <<
") or weights (" << cWeights.size() <<
")." <<
exit(
FatalError);
778 return metisDecomp(decompositionDict_, mesh_)
779 .decompose(globalCellCells, cellCentres, cWeights);
792 List<int> options(3, 0);
798 Field<int> cellWeights;
801 Field<int> faceWeights;
805 scalar minWeights =
gMin(cWeights);
806 if (cWeights.size() > 0)
812 "parMetisDecomp::decompose(const labelListList&"
813 ", const pointField&, const scalarField&)"
814 ) <<
"Illegal minimum weight " << minWeights
818 if (cWeights.size() != globalCellCells.size())
822 "parMetisDecomp::decompose(const labelListList&"
823 ", const pointField&, const scalarField&)"
824 ) <<
"Number of cell weights " << cWeights.size()
825 <<
" does not equal number of cells " << globalCellCells.size()
830 cellWeights.setSize(cWeights.size());
833 cellWeights[i] = int(cWeights[i]/minWeights);
839 if (decompositionDict_.found(
"metisCoeffs"))
841 const dictionary& metisCoeffs =
842 decompositionDict_.subDict(
"metisCoeffs");
844 if (metisCoeffs.readIfPresent(
"options", options))
846 Info<<
"Using Metis options " << options
849 if (options.size() != 3)
853 "parMetisDecomp::decompose(const labelListList&"
854 ", const pointField&, const scalarField&)"
855 ) <<
"Number of options " << options.size()
863 List<int> finalDecomp;
880 decomp[i] = finalDecomp[i];
886 void Foam::parMetisDecomp::calcMetisDistributedCSR
888 const polyMesh& mesh,
896 globalIndex globalCells(mesh.nCells());
906 const labelList& faceOwner = mesh.faceOwner();
907 const labelList& faceNeighbour = mesh.faceNeighbour();
908 const polyBoundaryMesh& patches = mesh.boundaryMesh();
914 List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
918 const polyPatch& pp = patches[patchI];
922 label faceI = pp.start();
923 label bFaceI = pp.start() - mesh.nInternalFaces();
927 globalNeighbour[bFaceI++] = globalCells.toGlobal
943 List<int> nFacesPerCell(mesh.nCells(), 0);
946 label nCoupledFaces = 0;
948 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
950 nFacesPerCell[faceOwner[faceI]]++;
951 nFacesPerCell[faceNeighbour[faceI]]++;
954 HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
958 const polyPatch& pp = patches[patchI];
962 label faceI = pp.start();
966 label own = faceOwner[faceI];
967 label globalNei = globalNeighbour[faceI-mesh.nInternalFaces()];
969 if (cellPair.insert(edge(own, globalNei)))
972 nFacesPerCell[faceOwner[faceI]]++;
983 xadj.setSize(mesh.nCells()+1);
987 for (label cellI = 0; cellI < mesh.nCells(); cellI++)
989 xadj[cellI] = freeAdj;
991 freeAdj += nFacesPerCell[cellI];
993 xadj[mesh.nCells()] = freeAdj;
1000 adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
1005 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
1007 label own = faceOwner[faceI];
1008 label nei = faceNeighbour[faceI];
1010 adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
1011 adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
1017 const polyPatch& pp = patches[patchI];
1021 label faceI = pp.start();
1022 label bFaceI = pp.start()-mesh.nInternalFaces();
1026 label own = faceOwner[faceI];
1027 label globalNei = globalNeighbour[bFaceI];
1029 if (cellPair.insert(edge(own, globalNei)))
1031 adjncy[xadj[own] + nFacesPerCell[own]++] =
1032 globalNeighbour[bFaceI];