Merge branch 'master' into particleInteractions

Conflicts:
	src/parallel/decompose/decompositionMethods/Make/options
This commit is contained in:
graham
2010-03-23 11:18:53 +00:00
115 changed files with 5324 additions and 2119 deletions

View File

@ -9,7 +9,7 @@ wmake libso metisDecomp
if [ -d "$FOAM_MPI_LIBBIN" ]
then
( WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB; wmake libso parMetisDecomp )
( WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB; wmake libso ptscotchDecomp )
fi
wmake libso decompositionMethods

View File

@ -4,7 +4,7 @@ set -x
wmakeLnInclude decompositionMethods
wmakeLnInclude metisDecomp
wmakeLnInclude parMetisDecomp
wmakeLnInclude scotchDecomp
wmakeLnInclude ptscotchDecomp
# ----------------------------------------------------------------- end-of-file

View File

@ -1,10 +1,4 @@
EXE_INC =
LIB_LIBS = \
-L$(FOAM_LIBBIN)/dummy \
-L$(FOAM_MPI_LIBBIN) \
-lscotchDecomp \
-lmetisDecomp \
-lparMetisDecomp
LIB_LIBS =

View File

@ -28,6 +28,9 @@ InClass
\*---------------------------------------------------------------------------*/
#include "decompositionMethod.H"
#include "globalIndex.H"
#include "cyclicPolyPatch.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -156,6 +159,18 @@ Foam::labelList Foam::decompositionMethod::decompose
}
Foam::labelList Foam::decompositionMethod::decompose
(
const labelListList& globalCellCells,
const pointField& cc
)
{
scalarField cWeights(0);
return decompose(globalCellCells, cc, cWeights);
}
void Foam::decompositionMethod::calcCellCells
(
const polyMesh& mesh,
@ -201,15 +216,284 @@ void Foam::decompositionMethod::calcCellCells
}
Foam::labelList Foam::decompositionMethod::decompose
void Foam::decompositionMethod::calcCSR
(
const labelListList& globalCellCells,
const pointField& cc
const polyMesh& mesh,
List<int>& adjncy,
List<int>& xadj
)
{
scalarField cWeights(0);
// Make Metis CSR (Compressed Storage Format) storage
// adjncy : contains neighbours (= edges in graph)
// xadj(celli) : start of information in adjncy for celli
return decompose(globalCellCells, cc, cWeights);
xadj.setSize(mesh.nCells()+1);
// Initialise the number of internal faces of the cells to twice the
// number of internal faces
label nInternalFaces = 2*mesh.nInternalFaces();
// Check the boundary for coupled patches and add to the number of
// internal faces
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
forAll(pbm, patchi)
{
if (isA<cyclicPolyPatch>(pbm[patchi]))
{
nInternalFaces += pbm[patchi].size();
}
}
// Create the adjncy array the size of the total number of internal and
// coupled faces
adjncy.setSize(nInternalFaces);
// Fill in xadj
// ~~~~~~~~~~~~
label freeAdj = 0;
for (label cellI = 0; cellI < mesh.nCells(); cellI++)
{
xadj[cellI] = freeAdj;
const labelList& cFaces = mesh.cells()[cellI];
forAll(cFaces, i)
{
label faceI = cFaces[i];
if
(
mesh.isInternalFace(faceI)
|| isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
)
{
freeAdj++;
}
}
}
xadj[mesh.nCells()] = freeAdj;
// Fill in adjncy
// ~~~~~~~~~~~~~~
labelList nFacesPerCell(mesh.nCells(), 0);
// Internal faces
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = mesh.faceOwner()[faceI];
label nei = mesh.faceNeighbour()[faceI];
adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
}
// Coupled faces. Only cyclics done.
forAll(pbm, patchi)
{
if (isA<cyclicPolyPatch>(pbm[patchi]))
{
const unallocLabelList& faceCells = pbm[patchi].faceCells();
label sizeby2 = faceCells.size()/2;
for (label facei=0; facei<sizeby2; facei++)
{
label own = faceCells[facei];
label nei = faceCells[facei + sizeby2];
adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
}
}
}
}
// From cell-cell connections to Metis format (like CompactListList)
void Foam::decompositionMethod::calcCSR
(
const labelListList& cellCells,
List<int>& adjncy,
List<int>& xadj
)
{
// Count number of internal faces
label nConnections = 0;
forAll(cellCells, coarseI)
{
nConnections += cellCells[coarseI].size();
}
// Create the adjncy array as twice the size of the total number of
// internal faces
adjncy.setSize(nConnections);
xadj.setSize(cellCells.size()+1);
// Fill in xadj
// ~~~~~~~~~~~~
label freeAdj = 0;
forAll(cellCells, coarseI)
{
xadj[coarseI] = freeAdj;
const labelList& cCells = cellCells[coarseI];
forAll(cCells, i)
{
adjncy[freeAdj++] = cCells[i];
}
}
xadj[cellCells.size()] = freeAdj;
}
void Foam::decompositionMethod::calcDistributedCSR
(
const polyMesh& mesh,
List<int>& adjncy,
List<int>& xadj
)
{
// Create global cell numbers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
globalIndex globalCells(mesh.nCells());
//
// Make Metis Distributed CSR (Compressed Storage Format) storage
// adjncy : contains cellCells (= edges in graph)
// xadj(celli) : start of information in adjncy for celli
//
const labelList& faceOwner = mesh.faceOwner();
const labelList& faceNeighbour = mesh.faceNeighbour();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Get renumbered owner on other side of coupled faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
label bFaceI = pp.start() - mesh.nInternalFaces();
forAll(pp, i)
{
globalNeighbour[bFaceI++] = globalCells.toGlobal
(
faceOwner[faceI++]
);
}
}
}
// Get the cell on the other side of coupled patches
syncTools::swapBoundaryFaceList(mesh, globalNeighbour, false);
// Count number of faces (internal + coupled)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Number of faces per cell
List<int> nFacesPerCell(mesh.nCells(), 0);
// Number of coupled faces
label nCoupledFaces = 0;
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
nFacesPerCell[faceOwner[faceI]]++;
nFacesPerCell[faceNeighbour[faceI]]++;
}
// Handle coupled faces
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
forAll(pp, i)
{
nCoupledFaces++;
nFacesPerCell[faceOwner[faceI++]]++;
}
}
}
// Fill in xadj
// ~~~~~~~~~~~~
xadj.setSize(mesh.nCells()+1);
int freeAdj = 0;
for (label cellI = 0; cellI < mesh.nCells(); cellI++)
{
xadj[cellI] = freeAdj;
freeAdj += nFacesPerCell[cellI];
}
xadj[mesh.nCells()] = freeAdj;
// Fill in adjncy
// ~~~~~~~~~~~~~~
adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
nFacesPerCell = 0;
// For internal faces is just offsetted owner and neighbour
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = faceOwner[faceI];
label nei = faceNeighbour[faceI];
adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
}
// For boundary faces is offsetted coupled neighbour
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
label bFaceI = pp.start()-mesh.nInternalFaces();
forAll(pp, i)
{
label own = faceOwner[faceI];
adjncy[xadj[own] + nFacesPerCell[own]++] =
globalNeighbour[bFaceI];
faceI++;
bFaceI++;
}
}
}
}

View File

@ -66,6 +66,31 @@ protected:
labelListList& cellCells
);
// From mesh to compact row storage format
// (like CompactListList)
static void calcCSR
(
const polyMesh& mesh,
List<int>& adjncy,
List<int>& xadj
);
// From cell-cell connections to compact row storage format
// (like CompactListList)
static void calcCSR
(
const labelListList& cellCells,
List<int>& adjncy,
List<int>& xadj
);
static void calcDistributedCSR
(
const polyMesh& mesh,
List<int>& adjncy,
List<int>& xadj
);
private:
// Private Member Functions

View File

@ -1,7 +1,6 @@
EXE_INC = \
-I$(WM_THIRD_PARTY_DIR)/metis-5.0pre2/include \
-I../decompositionMethods/lnInclude \
-I../scotchDecomp/lnInclude
-I../decompositionMethods/lnInclude
LIB_LIBS = \
-lmetis \

View File

@ -28,7 +28,6 @@ License
#include "addToRunTimeSelectionTable.H"
#include "floatScalar.H"
#include "Time.H"
#include "scotchDecomp.H"
extern "C"
{
@ -340,7 +339,7 @@ Foam::labelList Foam::metisDecomp::decompose
List<int> adjncy;
List<int> xadj;
scotchDecomp::calcCSR(mesh_, adjncy, xadj);
calcCSR(mesh_, adjncy, xadj);
// Decompose using default weights
List<int> finalDecomp;
@ -390,7 +389,7 @@ Foam::labelList Foam::metisDecomp::decompose
cellCells
);
scotchDecomp::calcCSR(cellCells, adjncy, xadj);
calcCSR(cellCells, adjncy, xadj);
}
// Decompose using default weights
@ -435,7 +434,7 @@ Foam::labelList Foam::metisDecomp::decompose
List<int> adjncy;
List<int> xadj;
scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
calcCSR(globalCellCells, adjncy, xadj);
// Decompose using default weights

View File

@ -1,3 +0,0 @@
parMetisDecomp.C
LIB = $(FOAM_MPI_LIBBIN)/libparMetisDecomp

View File

@ -1,14 +0,0 @@
include $(RULES)/mplib$(WM_MPLIB)
EXE_INC = \
$(PFLAGS) $(PINC) \
-I$(WM_THIRD_PARTY_DIR)/ParMetis-3.1/ParMETISLib \
-I$(WM_THIRD_PARTY_DIR)/ParMetis-3.1 \
-I../decompositionMethods/lnInclude \
-I../metisDecomp/lnInclude \
-I../scotchDecomp/lnInclude
LIB_LIBS = \
-L$(FOAM_MPI_LIBBIN) \
-lmetis-parmetis \
-lparmetis

File diff suppressed because it is too large Load Diff

View File

@ -1,76 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "parMetisDecomp.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Insert at front of list
template<class Type>
void Foam::parMetisDecomp::prepend
(
const UList<Type>& extraLst,
List<Type>& lst
)
{
label nExtra = extraLst.size();
// Make space for initial elements
lst.setSize(lst.size() + nExtra);
for (label i = lst.size()-1; i >= nExtra; i--)
{
lst[i] = lst[i-nExtra];
}
// Insert at front
forAll(extraLst, i)
{
lst[i] = extraLst[i];
}
}
// Insert at back of list
template<class Type>
void Foam::parMetisDecomp::append
(
const UList<Type>& extraLst,
List<Type>& lst
)
{
label sz = lst.size();
// Make space for initial elements
lst.setSize(sz + extraLst.size());
// Insert at back
forAll(extraLst, i)
{
lst[sz++] = extraLst[i];
}
}
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
ptscotchDecomp.C
LIB = $(FOAM_MPI_LIBBIN)/libptscotchDecomp

View File

@ -0,0 +1,10 @@
include $(RULES)/mplib$(WM_MPLIB)
EXE_INC = \
$(PFLAGS) $(PINC) \
-I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \
-I/usr/include/scotch \
-I../decompositionMethods/lnInclude
LIB_LIBS = \
-L$(FOAM_MPI_LIBBIN) -lptscotch -lptscotcherrexit

View File

@ -0,0 +1,592 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
From scotch forum:
By: Francois PELLEGRINI RE: Graph mapping 'strategy' string [ reply ]
2008-08-22 10:09 Strategy handling in Scotch is a bit tricky. In order
not to be confused, you must have a clear view of how they are built.
Here are some rules:
1- Strategies are made up of "methods" which are combined by means of
"operators".
2- A method is of the form "m{param=value,param=value,...}", where "m"
is a single character (this is your first error: "f" is a method name,
not a parameter name).
3- There exist different sort of strategies : bipartitioning strategies,
mapping strategies, ordering strategies, which cannot be mixed. For
instance, you cannot build a bipartitioning strategy and feed it to a
mapping method (this is your second error).
To use the "mapCompute" routine, you must create a mapping strategy, not
a bipartitioning one, and so use stratGraphMap() and not
stratGraphBipart(). Your mapping strategy should however be based on the
"recursive bipartitioning" method ("b"). For instance, a simple (and
hence not very efficient) mapping strategy can be :
"b{sep=f}"
which computes mappings with the recursive bipartitioning method "b",
this latter using the Fiduccia-Mattheyses method "f" to compute its
separators.
If you want an exact partition (see your previous post), try
"b{sep=fx}".
However, these strategies are not the most efficient, as they do not
make use of the multi-level framework.
To use the multi-level framework, try for instance:
"b{sep=m{vert=100,low=h,asc=f}x}"
The current default mapping strategy in Scotch can be seen by using the
"-vs" option of program gmap. It is, to date:
b
{
job=t,
map=t,
poli=S,
sep=
(
m
{
asc=b
{
bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
org=f{move=80,pass=-1,bal=0.005},
width=3
},
low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
type=h,
vert=80,
rat=0.8
}
| m
{
asc=b
{
bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
org=f{move=80,pass=-1,bal=0.005},
width=3
},
low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
type=h,
vert=80,
rat=0.8
}
)
}
\*---------------------------------------------------------------------------*/
#include "ptscotchDecomp.H"
#include "addToRunTimeSelectionTable.H"
#include "Time.H"
#include "OFstream.H"
extern "C"
{
#include <stdio.h>
#include "mpi.h"
#include "ptscotch.h"
}
// Hack: scotch generates floating point errors so need to switch of error
// trapping!
#if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
# define LINUX
#endif
#if defined(LINUX) && defined(__GNUC__)
# define LINUX_GNUC
#endif
#ifdef LINUX_GNUC
# ifndef __USE_GNU
# define __USE_GNU
# endif
# include <fenv.h>
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(ptscotchDecomp, 0);
addToRunTimeSelectionTable
(
decompositionMethod,
ptscotchDecomp,
dictionaryMesh
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::ptscotchDecomp::check(const int retVal, const char* str)
{
if (retVal)
{
FatalErrorIn("ptscotchDecomp::decompose(..)")
<< "Call to scotch routine " << str << " failed."
<< exit(FatalError);
}
}
// Call scotch with options from dictionary.
Foam::label Foam::ptscotchDecomp::decompose
(
List<int>& adjncy,
List<int>& xadj,
const scalarField& cWeights,
List<int>& finalDecomp
)
{
// // Dump graph
// if (decompositionDict_.found("ptscotchCoeffs"))
// {
// const dictionary& scotchCoeffs =
// decompositionDict_.subDict("ptscotchCoeffs");
//
// if (scotchCoeffs.found("writeGraph"))
// {
// Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
//
// if (writeGraph)
// {
// OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
//
// Info<< "Dumping Scotch graph file to " << str.name() << endl
// << "Use this in combination with gpart." << endl;
//
// label version = 0;
// str << version << nl;
// // Numer of vertices
// str << xadj.size()-1 << ' ' << adjncy.size() << nl;
// // Numbering starts from 0
// label baseval = 0;
// // Has weights?
// label hasEdgeWeights = 0;
// label hasVertexWeights = 0;
// label numericflag = 10*hasEdgeWeights+hasVertexWeights;
// str << baseval << ' ' << numericflag << nl;
// for (label cellI = 0; cellI < xadj.size()-1; cellI++)
// {
// label start = xadj[cellI];
// label end = xadj[cellI+1];
// str << end-start;
//
// for (label i = start; i < end; i++)
// {
// str << ' ' << adjncy[i];
// }
// str << nl;
// }
// }
// }
// }
// Strategy
// ~~~~~~~~
// Default.
SCOTCH_Strat stradat;
check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
if (decompositionDict_.found("scotchCoeffs"))
{
const dictionary& scotchCoeffs =
decompositionDict_.subDict("scotchCoeffs");
string strategy;
if (scotchCoeffs.readIfPresent("strategy", strategy))
{
if (debug)
{
Info<< "ptscotchDecomp : Using strategy " << strategy << endl;
}
SCOTCH_stratDgraphMap(&stradat, strategy.c_str());
//fprintf(stdout, "S\tStrat=");
//SCOTCH_stratSave(&stradat, stdout);
//fprintf(stdout, "\n");
}
}
// Graph
// ~~~~~
List<int> velotab;
// Check for externally provided cellweights and if so initialise weights
scalar minWeights = gMin(cWeights);
if (cWeights.size() > 0)
{
if (minWeights <= 0)
{
WarningIn
(
"ptscotchDecomp::decompose"
"(const pointField&, const scalarField&)"
) << "Illegal minimum weight " << minWeights
<< endl;
}
if (cWeights.size() != xadj.size()-1)
{
FatalErrorIn
(
"ptscotchDecomp::decompose"
"(const pointField&, const scalarField&)"
) << "Number of cell weights " << cWeights.size()
<< " does not equal number of cells " << xadj.size()-1
<< exit(FatalError);
}
// Convert to integers.
velotab.setSize(cWeights.size());
forAll(velotab, i)
{
velotab[i] = int(cWeights[i]/minWeights);
}
}
SCOTCH_Dgraph grafdat;
check(SCOTCH_dgraphInit(&grafdat, MPI_COMM_WORLD), "SCOTCH_dgraphInit");
check
(
SCOTCH_dgraphBuild
(
&grafdat, // grafdat
0, // baseval, c-style numbering
xadj.size()-1, // vertlocnbr, nCells
xadj.size()-1, // vertlocmax
const_cast<SCOTCH_Num*>(xadj.begin()), // vertloctab, start index per cell into
// adjncy
&xadj[1], // vendloctab, end index ,,
velotab.begin(), // veloloctab, vertex weights
NULL, // vlblloctab
adjncy.size(), // edgelocnbr, number of arcs
adjncy.size(), // edgelocsiz
adjncy.begin(), // edgeloctab
NULL, // edgegsttab
NULL // edlotab, edge weights
),
"SCOTCH_dgraphBuild"
);
check(SCOTCH_dgraphCheck(&grafdat), "SCOTCH_dgraphCheck");
// Architecture
// ~~~~~~~~~~~~
// (fully connected network topology since using switch)
SCOTCH_Arch archdat;
check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
List<label> processorWeights;
if (decompositionDict_.found("scotchCoeffs"))
{
const dictionary& scotchCoeffs =
decompositionDict_.subDict("scotchCoeffs");
scotchCoeffs.readIfPresent("processorWeights", processorWeights);
}
if (processorWeights.size())
{
if (debug)
{
Info<< "ptscotchDecomp : Using procesor weights " << processorWeights
<< endl;
}
check
(
SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
"SCOTCH_archCmpltw"
);
}
else
{
check
(
SCOTCH_archCmplt(&archdat, nProcessors_),
"SCOTCH_archCmplt"
);
}
//SCOTCH_Mapping mapdat;
//SCOTCH_dgraphMapInit(&grafdat, &mapdat, &archdat, NULL);
//SCOTCH_dgraphMapCompute(&grafdat, &mapdat, &stradat); /*Perform mapping*/
//SCOTCHdgraphMapExit(&grafdat, &mapdat);
// Hack:switch off fpu error trapping
# ifdef LINUX_GNUC
int oldExcepts = fedisableexcept
(
FE_DIVBYZERO
| FE_INVALID
| FE_OVERFLOW
);
# endif
finalDecomp.setSize(xadj.size()-1);
finalDecomp = 0;
check
(
SCOTCH_dgraphMap
(
&grafdat,
&archdat,
&stradat, // const SCOTCH_Strat *
finalDecomp.begin() // parttab
),
"SCOTCH_graphMap"
);
# ifdef LINUX_GNUC
feenableexcept(oldExcepts);
# endif
//finalDecomp.setSize(xadj.size()-1);
//check
//(
// SCOTCH_dgraphPart
// (
// &grafdat,
// nProcessors_, // partnbr
// &stradat, // const SCOTCH_Strat *
// finalDecomp.begin() // parttab
// ),
// "SCOTCH_graphPart"
//);
// Release storage for graph
SCOTCH_dgraphExit(&grafdat);
// Release storage for strategy
SCOTCH_stratExit(&stradat);
// Release storage for network topology
SCOTCH_archExit(&archdat);
return 0;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::ptscotchDecomp::ptscotchDecomp
(
const dictionary& decompositionDict,
const polyMesh& mesh
)
:
decompositionMethod(decompositionDict),
mesh_(mesh)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::labelList Foam::ptscotchDecomp::decompose
(
const pointField& points,
const scalarField& pointWeights
)
{
if (points.size() != mesh_.nCells())
{
FatalErrorIn
(
"ptscotchDecomp::decompose(const pointField&, const scalarField&)"
)
<< "Can use this decomposition method only for the whole mesh"
<< endl
<< "and supply one coordinate (cellCentre) for every cell." << endl
<< "The number of coordinates " << points.size() << endl
<< "The number of cells in the mesh " << mesh_.nCells()
<< exit(FatalError);
}
// // For running sequential ...
// if (Pstream::nProcs() <= 1)
// {
// return scotchDecomp(decompositionDict_, mesh_)
// .decompose(points, pointWeights);
// }
// Make Metis CSR (Compressed Storage Format) storage
// adjncy : contains neighbours (= edges in graph)
// xadj(celli) : start of information in adjncy for celli
// Connections
Field<int> adjncy;
// Offsets into adjncy
Field<int> xadj;
calcDistributedCSR
(
mesh_,
adjncy,
xadj
);
// Decompose using default weights
List<int> finalDecomp;
decompose(adjncy, xadj, pointWeights, finalDecomp);
// Copy back to labelList
labelList decomp(finalDecomp.size());
forAll(decomp, i)
{
decomp[i] = finalDecomp[i];
}
return decomp;
}
Foam::labelList Foam::ptscotchDecomp::decompose
(
const labelList& agglom,
const pointField& agglomPoints,
const scalarField& pointWeights
)
{
if (agglom.size() != mesh_.nCells())
{
FatalErrorIn
(
"ptscotchDecomp::decompose(const labelList&, const pointField&)"
) << "Size of cell-to-coarse map " << agglom.size()
<< " differs from number of cells in mesh " << mesh_.nCells()
<< exit(FatalError);
}
// // For running sequential ...
// if (Pstream::nProcs() <= 1)
// {
// return scotchDecomp(decompositionDict_, mesh_)
// .decompose(agglom, agglomPoints, pointWeights);
// }
// Make Metis CSR (Compressed Storage Format) storage
// adjncy : contains neighbours (= edges in graph)
// xadj(celli) : start of information in adjncy for celli
List<int> adjncy;
List<int> xadj;
{
// Get cellCells on coarse mesh.
labelListList cellCells;
calcCellCells
(
mesh_,
agglom,
agglomPoints.size(),
cellCells
);
calcCSR(cellCells, adjncy, xadj);
}
// Decompose using weights
List<int> finalDecomp;
decompose(adjncy, xadj, pointWeights, finalDecomp);
// Rework back into decomposition for original mesh_
labelList fineDistribution(agglom.size());
forAll(fineDistribution, i)
{
fineDistribution[i] = finalDecomp[agglom[i]];
}
return fineDistribution;
}
Foam::labelList Foam::ptscotchDecomp::decompose
(
const labelListList& globalCellCells,
const pointField& cellCentres,
const scalarField& cWeights
)
{
if (cellCentres.size() != globalCellCells.size())
{
FatalErrorIn
(
"ptscotchDecomp::decompose(const pointField&, const labelListList&)"
) << "Inconsistent number of cells (" << globalCellCells.size()
<< ") and number of cell centres (" << cellCentres.size()
<< ")." << exit(FatalError);
}
// // For running sequential ...
// if (Pstream::nProcs() <= 1)
// {
// return scotchDecomp(decompositionDict_, mesh_)
// .decompose(globalCellCells, cellCentres, cWeights);
// }
// Make Metis CSR (Compressed Storage Format) storage
// adjncy : contains neighbours (= edges in graph)
// xadj(celli) : start of information in adjncy for celli
List<int> adjncy;
List<int> xadj;
calcCSR(globalCellCells, adjncy, xadj);
// Decompose using weights
List<int> finalDecomp;
decompose(adjncy, xadj, cWeights, finalDecomp);
// Copy back to labelList
labelList decomp(finalDecomp.size());
forAll(decomp, i)
{
decomp[i] = finalDecomp[i];
}
return decomp;
}
// ************************************************************************* //

View File

@ -23,17 +23,18 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::parMetisDecomp
Foam::ptscotchDecomp
Description
PTScotch domain decomposition
SourceFiles
parMetisDecomp.C
ptscotchDecomp.C
\*---------------------------------------------------------------------------*/
#ifndef parMetisDecomp_H
#define parMetisDecomp_H
#ifndef ptscotchDecomp_H
#define ptscotchDecomp_H
#include "decompositionMethod.H"
@ -41,10 +42,10 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class parMetisDecomp Declaration
Class ptscotchDecomp Declaration
\*---------------------------------------------------------------------------*/
class parMetisDecomp
class ptscotchDecomp
:
public decompositionMethod
{
@ -55,41 +56,32 @@ class parMetisDecomp
// Private Member Functions
//- Insert list in front of list.
template<class Type>
static void prepend(const UList<Type>&, List<Type>&);
//- Insert list at end of list.
template<class Type>
static void append(const UList<Type>&, List<Type>&);
//- Check and print error message
static void check(const int, const char*);
label decompose
(
Field<int>& xadj,
Field<int>& adjncy,
const pointField& cellCentres,
Field<int>& cellWeights,
Field<int>& faceWeights,
const List<int>& options,
List<int>& adjncy,
List<int>& xadj,
const scalarField& cWeights,
List<int>& finalDecomp
);
//- Disallow default bitwise copy construct and assignment
void operator=(const parMetisDecomp&);
parMetisDecomp(const parMetisDecomp&);
void operator=(const ptscotchDecomp&);
ptscotchDecomp(const ptscotchDecomp&);
public:
//- Runtime type information
TypeName("parMetis");
TypeName("ptscotch");
// Constructors
//- Construct given the decomposition dictionary and mesh
parMetisDecomp
ptscotchDecomp
(
const dictionary& decompositionDict,
const polyMesh& mesh
@ -98,23 +90,20 @@ public:
// Destructor
~parMetisDecomp()
virtual ~ptscotchDecomp()
{}
// Member Functions
//- parMetis handles Foam processor boundaries
virtual bool parallelAware() const
{
// ptscotch does not know about proc boundaries
return true;
}
//- Return for every coordinate the wanted processor number. Use the
// mesh connectivity (if needed)
// Weights get normalised so the minimum value is 1 before truncation
// to an integer so the weights should be multiples of the minimum
// value. The overall sum of weights might otherwise overflow.
virtual labelList decompose
(
const pointField& points,
@ -125,10 +114,9 @@ public:
// passed agglomeration map (from fine to coarse cells) and coarse cell
// location. Can be overridden by decomposers that provide this
// functionality natively.
// See note on weights above.
virtual labelList decompose
(
const labelList& cellToRegion,
const labelList& agglom,
const pointField& regionPoints,
const scalarField& regionWeights
);
@ -140,7 +128,6 @@ public:
// from 0 at processor0 and then incrementing all through the
// processors)
// - the connections are across coupled patches
// See note on weights above.
virtual labelList decompose
(
const labelListList& globalCellCells,
@ -148,13 +135,6 @@ public:
const scalarField& cWeights
);
//- Helper to convert mesh connectivity into distributed CSR
static void calcMetisDistributedCSR
(
const polyMesh&,
List<int>& adjncy,
List<int>& xadj
);
};
@ -164,12 +144,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "parMetisDecompTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -474,7 +474,7 @@ Foam::labelList Foam::scotchDecomp::decompose
{
FatalErrorIn
(
"parMetisDecomp::decompose(const labelList&, const pointField&)"
"scotchDecomp::decompose(const labelList&, const pointField&)"
) << "Size of cell-to-coarse map " << agglom.size()
<< " differs from number of cells in mesh " << mesh_.nCells()
<< exit(FatalError);
@ -556,144 +556,4 @@ Foam::labelList Foam::scotchDecomp::decompose
}
void Foam::scotchDecomp::calcCSR
(
const polyMesh& mesh,
List<int>& adjncy,
List<int>& xadj
)
{
// Make Metis CSR (Compressed Storage Format) storage
// adjncy : contains neighbours (= edges in graph)
// xadj(celli) : start of information in adjncy for celli
xadj.setSize(mesh.nCells()+1);
// Initialise the number of internal faces of the cells to twice the
// number of internal faces
label nInternalFaces = 2*mesh.nInternalFaces();
// Check the boundary for coupled patches and add to the number of
// internal faces
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
forAll(pbm, patchi)
{
if (isA<cyclicPolyPatch>(pbm[patchi]))
{
nInternalFaces += pbm[patchi].size();
}
}
// Create the adjncy array the size of the total number of internal and
// coupled faces
adjncy.setSize(nInternalFaces);
// Fill in xadj
// ~~~~~~~~~~~~
label freeAdj = 0;
for (label cellI = 0; cellI < mesh.nCells(); cellI++)
{
xadj[cellI] = freeAdj;
const labelList& cFaces = mesh.cells()[cellI];
forAll(cFaces, i)
{
label faceI = cFaces[i];
if
(
mesh.isInternalFace(faceI)
|| isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
)
{
freeAdj++;
}
}
}
xadj[mesh.nCells()] = freeAdj;
// Fill in adjncy
// ~~~~~~~~~~~~~~
labelList nFacesPerCell(mesh.nCells(), 0);
// Internal faces
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = mesh.faceOwner()[faceI];
label nei = mesh.faceNeighbour()[faceI];
adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
}
// Coupled faces. Only cyclics done.
forAll(pbm, patchi)
{
if (isA<cyclicPolyPatch>(pbm[patchi]))
{
const unallocLabelList& faceCells = pbm[patchi].faceCells();
label sizeby2 = faceCells.size()/2;
for (label facei=0; facei<sizeby2; facei++)
{
label own = faceCells[facei];
label nei = faceCells[facei + sizeby2];
adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
}
}
}
}
// From cell-cell connections to Metis format (like CompactListList)
void Foam::scotchDecomp::calcCSR
(
const labelListList& cellCells,
List<int>& adjncy,
List<int>& xadj
)
{
// Count number of internal faces
label nConnections = 0;
forAll(cellCells, coarseI)
{
nConnections += cellCells[coarseI].size();
}
// Create the adjncy array as twice the size of the total number of
// internal faces
adjncy.setSize(nConnections);
xadj.setSize(cellCells.size()+1);
// Fill in xadj
// ~~~~~~~~~~~~
label freeAdj = 0;
forAll(cellCells, coarseI)
{
xadj[coarseI] = freeAdj;
const labelList& cCells = cellCells[coarseI];
forAll(cCells, i)
{
adjncy[freeAdj++] = cCells[i];
}
}
xadj[cellCells.size()] = freeAdj;
}
// ************************************************************************* //

View File

@ -139,26 +139,6 @@ public:
const pointField& cc,
const scalarField& cWeights
);
//- Helper to convert local connectivity (supplied as owner,neighbour)
// into CSR (Metis,scotch) storage. Does cyclics but not processor
// patches
static void calcCSR
(
const polyMesh& mesh,
List<int>& adjncy,
List<int>& xadj
);
//- Helper to convert connectivity (supplied as cellcells) into
// CSR (Metis,scotch) storage
static void calcCSR
(
const labelListList& globalCellCells,
List<int>& adjncy,
List<int>& xadj
);
};