Dune Core Modules (2.5.0)

parmetisgridpartitioner.hh
Go to the documentation of this file.
1#ifndef DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
2#define DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
3
8#include <algorithm>
9#include <vector>
10
13
14#include <dune/geometry/referenceelements.hh>
15
18
19#if HAVE_PARMETIS
20
21#include <parmetis.h>
22
23// only enable for ParMETIS because the implementation uses functions that
24// are not emulated by scotch
25#ifdef PARMETIS_MAJOR_VERSION
26
27namespace Dune
28{
29
34 template<class GridView>
35 struct ParMetisGridPartitioner {
36
37 // define index type as provided by ParMETIS
38#if PARMETIS_MAJOR_VERSION > 3
39 typedef idx_t idx_type;
40 typedef ::real_t real_type;
41#else
42 typedef int idx_type;
43 typedef float real_type;
44#endif // PARMETIS_MAJOR_VERSION > 3
45
46 typedef typename GridView::template Codim<0>::Iterator ElementIterator;
47 typedef typename GridView::template Codim<0>::template Partition<Interior_Partition>::Iterator InteriorElementIterator;
48 typedef typename GridView::IntersectionIterator IntersectionIterator;
49
50 enum {
51 dimension = GridView::dimension
52 };
53
54
65 static std::vector<unsigned> partition(const GridView& gv, const Dune::MPIHelper& mpihelper) {
66 const unsigned numElements = gv.size(0);
67
68 std::vector<unsigned> part(numElements);
69
70 // Setup parameters for ParMETIS
71 idx_type wgtflag = 0; // we don't use weights
72 idx_type numflag = 0; // we are using C-style arrays
73 idx_type ncon = 1; // number of balance constraints
74 idx_type ncommonnodes = 2; // number of nodes elements must have in common to be considered adjacent to each other
75 idx_type options[4] = {0, 0, 0, 0}; // use default values for random seed, output and coupling
76 idx_type edgecut; // will store number of edges cut by partition
77 idx_type nparts = mpihelper.size(); // number of parts equals number of processes
78 std::vector<real_type> tpwgts(ncon*nparts, 1./nparts); // load per subdomain and weight (same load on every process)
79 std::vector<real_type> ubvec(ncon, 1.05); // weight tolerance (same weight tolerance for every weight there is)
80
81 // The difference elmdist[i+1] - elmdist[i] is the number of nodes that are on process i
82 std::vector<idx_type> elmdist(nparts+1);
83 elmdist[0] = 0;
84 std::fill(elmdist.begin()+1, elmdist.end(), gv.size(0)); // all elements are on process zero
85
86 // Create and fill arrays "eptr", where eptr[i] is the number of vertices that belong to the i-th element, and
87 // "eind" contains the vertex-numbers of the i-the element in eind[eptr[i]] to eind[eptr[i+1]-1]
88 std::vector<idx_type> eptr, eind;
89 int numVertices = 0;
90 eptr.push_back(numVertices);
91
92 for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>(); eIt != gv.template end<0, Interior_Partition>(); ++eIt) {
93 const size_t curNumVertices = ReferenceElements<double, dimension>::general(eIt->type()).size(dimension);
94
95 numVertices += curNumVertices;
96 eptr.push_back(numVertices);
97
98 for (size_t k = 0; k < curNumVertices; ++k)
99 eind.push_back(gv.indexSet().subIndex(*eIt, k, dimension));
100 }
101
102 // Partition mesh using ParMETIS
103 if (0 == mpihelper.rank()) {
105
106#if PARMETIS_MAJOR_VERSION >= 4
107 const int OK =
108#endif
109 ParMETIS_V3_PartMeshKway(elmdist.data(), eptr.data(), eind.data(), NULL, &wgtflag, &numflag,
110 &ncon, &ncommonnodes, &nparts, tpwgts.data(), ubvec.data(),
111 options, &edgecut, reinterpret_cast<idx_type*>(part.data()), &comm);
112
113#if PARMETIS_MAJOR_VERSION >= 4
114 if (OK != METIS_OK)
115 DUNE_THROW(Dune::Exception, "ParMETIS returned an error code.");
116#endif
117 }
118
119 return part;
120 }
121
133 static std::vector<unsigned> repartition(const GridView& gv, const Dune::MPIHelper& mpihelper, real_type& itr = 1000) {
134
135 // Create global index map
136 GlobalIndexSet<GridView> globalIndex(gv,0);
137
138 int numElements = std::distance(gv.template begin<0, Interior_Partition>(),
139 gv.template end<0, Interior_Partition>());
140
141 std::vector<unsigned> interiorPart(numElements);
142
143 // Setup parameters for ParMETIS
144 idx_type wgtflag = 0; // we don't use weights
145 idx_type numflag = 0; // we are using C-style arrays
146 idx_type ncon = 1; // number of balance constraints
147 idx_type options[4] = {0, 0, 0, 0}; // use default values for random seed, output and coupling
148 idx_type edgecut; // will store number of edges cut by partition
149 idx_type nparts = mpihelper.size(); // number of parts equals number of processes
150 std::vector<real_type> tpwgts(ncon*nparts, 1./nparts); // load per subdomain and weight (same load on every process)
151 std::vector<real_type> ubvec(ncon, 1.05); // weight tolerance (same weight tolerance for every weight there is)
152
153 MPI_Comm comm = Dune::MPIHelper::getCommunicator();
154
155 // Make the number of interior elements of each processor available to all processors
156 std::vector<int> offset(gv.comm().size());
157 std::fill(offset.begin(), offset.end(), 0);
158
159 gv.comm().template allgather<int>(&numElements, 1, offset.data());
160
161 // The difference vtxdist[i+1] - vtxdist[i] is the number of elements that are on process i
162 std::vector<idx_type> vtxdist(gv.comm().size()+1);
163 vtxdist[0] = 0;
164
165 for (unsigned int i=1; i<vtxdist.size(); ++i)
166 vtxdist[i] = vtxdist[i-1] + offset[i-1];
167
168 // Set up element adjacency lists
169 std::vector<idx_type> xadj, adjncy;
170 xadj.push_back(0);
171
172 for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>(); eIt != gv.template end<0, Interior_Partition>(); ++eIt)
173 {
174 size_t numNeighbors = 0;
175
176 for (IntersectionIterator iIt = gv.template ibegin(*eIt); iIt != gv.template iend(*eIt); ++iIt) {
177 if (iIt->neighbor()) {
178 adjncy.push_back(globalIndex.index(iIt->outside()));
179
180 ++numNeighbors;
181 }
182 }
183
184 xadj.push_back(xadj.back() + numNeighbors);
185 }
186
187#if PARMETIS_MAJOR_VERSION >= 4
188 const int OK =
189#endif
190 ParMETIS_V3_AdaptiveRepart(vtxdist.data(), xadj.data(), adjncy.data(), NULL, NULL, NULL,
191 &wgtflag, &numflag, &ncon, &nparts, tpwgts.data(), ubvec.data(),
192 &itr, options, &edgecut, reinterpret_cast<idx_type*>(interiorPart.data()), &comm);
193
194#if PARMETIS_MAJOR_VERSION >= 4
195 if (OK != METIS_OK)
196 DUNE_THROW(Dune::Exception, "ParMETIS returned error code " << OK);
197#endif
198
199 // At this point, interiorPart contains a target rank for each interior element, and they are sorted
200 // by the order in which the grid view traverses them. Now we need to do two things:
201 // a) Add additional dummy entries for the ghost elements
202 // b) Use the element index for the actual ordering. Since there may be different types of elements,
203 // we cannot use the index set directly, but have to go through a Mapper.
204
205 typedef MultipleCodimMultipleGeomTypeMapper<GridView, MCMGElementLayout> ElementMapper;
206 ElementMapper elementMapper(gv);
207
208 std::vector<unsigned int> part(gv.size(0));
209 std::fill(part.begin(), part.end(), 0);
210 unsigned int c = 0;
211 for (InteriorElementIterator eIt = gv.template begin<0, Interior_Partition>();
212 eIt != gv.template end<0, Interior_Partition>();
213 ++eIt)
214 {
215 part[elementMapper.index(*eIt)] = interiorPart[c++];
216 }
217 return part;
218 }
219 };
220
221} // namespace Dune
222
223#else // PARMETIS_MAJOR_VERSION
224#warning "You seem to be using the ParMETIS emulation layer of scotch, which does not work with this file."
225#endif
226
227#else // HAVE_PARMETIS
228#warning "PARMETIS was not found, please check your configuration"
229#endif
230
231#endif // DUNE_GRID_UTILITY_PARMETISGRIDPARTITIONER_HH
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
A real mpi helper.
Definition: mpihelper.hh:160
int size() const
return number of processes
Definition: mpihelper.hh:231
int rank() const
return rank of process
Definition: mpihelper.hh:227
static MPICommunicator getCommunicator()
get the default communicator
Definition: mpihelper.hh:181
static MPICommunicator getLocalCommunicator()
get a local communicator
Definition: mpihelper.hh:192
A few common exception classes.
Provides a globally unique index for all entities of a distributed Dune grid.
typename FieldTraits< Type >::real_type real_t
Convenient access to FieldTraits<Type>::real_type.
Definition: typetraits.hh:472
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Traits::IntersectionIterator IntersectionIterator
type of the intersection iterator
Definition: gridview.hh:87
@ dimension
The dimension of the grid.
Definition: gridview.hh:128
Mapper for multiple codim and multiple geometry types.
Helpers for dealing with MPI.
Dune namespace.
Definition: alignment.hh:11
static const ReferenceElement< ctype, dim > & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:757
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)