Dune Core Modules (2.4.1)

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