45 #ifndef _ZOLTAN2_ALGPARMA_HPP_ 46 #define _ZOLTAN2_ALGPARMA_HPP_ 73 #ifndef HAVE_ZOLTAN2_PARMA 79 template <
typename Adapter>
83 typedef typename Adapter::user_t
user_t;
86 const RCP<
const Comm<int> > &problemComm,
89 throw std::runtime_error(
90 "BUILD ERROR: ParMA requested but not compiled into Zoltan2.\n" 91 "Please set CMake flag Zoltan2_ENABLE_ParMA:BOOL=ON.");
98 #ifdef HAVE_ZOLTAN2_PARMA 102 #include <gmi_null.h> 104 #include <apfMesh2.h> 105 #include <apfNumbering.h> 108 #include <apfConvert.h> 109 #include <apfShape.h> 115 template <
typename Adapter>
121 typedef typename Adapter::lno_t
lno_t;
122 typedef typename Adapter::gno_t
gno_t;
123 typedef typename Adapter::scalar_t
scalar_t;
124 typedef typename Adapter::part_t
part_t;
125 typedef typename Adapter::user_t
user_t;
126 typedef typename Adapter::userCoord_t userCoord_t;
128 const RCP<const Environment> env;
129 const RCP<const Comm<int> > problemComm;
130 const RCP<const MeshAdapter<user_t> > adapter;
133 apf::Numbering* gids;
134 apf::Numbering* origin_part_ids;
135 std::map<gno_t, lno_t> mapping_elm_gids_index;
140 void setMPIComm(
const RCP<
const Comm<int> > &problemComm__) {
141 # ifdef HAVE_ZOLTAN2_MPI 142 mpicomm = Teuchos::getRawMpiComm(*problemComm__);
144 mpicomm = MPI_COMM_WORLD;
154 return apf::Mesh::VERTEX;
156 return apf::Mesh::EDGE;
160 return apf::Mesh::QUAD;
162 return apf::Mesh::TET;
164 return apf::Mesh::HEX;
165 else if (ttype==
PRISM)
170 throw std::runtime_error(
"APF does not support this topology type");
176 void setEntWeights(
int dim, apf::MeshTag* tag) {
178 for (
int i=0;i<m->getTagSize(tag);i++) {
179 apf::MeshIterator* itr = m->begin(dim);
180 apf::MeshEntity* ent;
181 const scalar_t* ws=NULL;
183 if (i<adapter->getNumWeightsPerOf(etype))
184 adapter->getWeightsViewOf(etype,ws,stride,i);
186 while ((ent= m->iterate(itr))) {
189 w =
static_cast<double>(ws[j]);
190 m->setDoubleTag(ent,tag,&w);
198 apf::MeshTag* setWeights(
bool vtx,
bool edge,
bool elm) {
201 num_ws = std::max(num_ws,adapter->getNumWeightsPerOf(
MESH_VERTEX));
203 num_ws = std::max(num_ws,adapter->getNumWeightsPerOf(
MESH_EDGE));
205 num_ws = std::max(num_ws,adapter->getNumWeightsPerOf(entityAPFtoZ2(m->getDimension())));
206 apf::MeshTag* tag = m->createDoubleTag(
"parma_weight",num_ws);
208 setEntWeights(0,tag);
210 setEntWeights(1,tag);
212 setEntWeights(m->getDimension(),tag);
219 void constructElements(
const gno_t* conn, lno_t nelem,
const lno_t* offsets,
222 apf::ModelEntity* interior = m->findModelEntity(m->getDimension(), 0);
223 for (lno_t i = 0; i < nelem; ++i) {
224 apf::Mesh::Type etype = topologyZ2toAPF(tops[i]);
226 for (
int j = offsets[i]; j < offsets[i+1]; ++j)
227 verts[j-offsets[i]] = globalToVert[conn[j]];
228 buildElement(m, interior, etype, verts);
231 int getMax(
const apf::GlobalToVert& globalToVert)
234 APF_CONST_ITERATE(apf::GlobalToVert, globalToVert, it)
235 max = std::max(max, it->first);
236 PCU_Max_Ints(&max, 1);
239 void constructResidence(apf::GlobalToVert& globalToVert)
241 int max = getMax(globalToVert);
243 int peers = PCU_Comm_Peers();
244 int quotient = total / peers;
245 int remainder = total % peers;
246 int mySize = quotient;
247 int self = PCU_Comm_Self();
248 if (
self == (peers - 1))
250 typedef std::vector< std::vector<int> > TmpParts;
251 TmpParts tmpParts(mySize);
255 APF_ITERATE(apf::GlobalToVert, globalToVert, it) {
257 int to = std::min(peers - 1, gid / quotient);
258 PCU_COMM_PACK(to, gid);
261 int myOffset =
self * quotient;
264 while (PCU_Comm_Receive()) {
266 PCU_COMM_UNPACK(gid);
267 int from = PCU_Comm_Sender();
268 tmpParts.at(gid - myOffset).push_back(from);
273 for (
int i = 0; i < mySize; ++i) {
274 std::vector<int>& parts = tmpParts[i];
275 for (
size_t j = 0; j < parts.size(); ++j) {
277 int gid = i + myOffset;
278 int nparts = parts.size();
279 PCU_COMM_PACK(to, gid);
280 PCU_COMM_PACK(to, nparts);
281 for (
size_t k = 0; k < parts.size(); ++k)
282 PCU_COMM_PACK(to, parts[k]);
289 while (PCU_Comm_Receive()) {
291 PCU_COMM_UNPACK(gid);
293 PCU_COMM_UNPACK(nparts);
294 apf::Parts residence;
295 for (
int i = 0; i < nparts; ++i) {
297 PCU_COMM_UNPACK(part);
298 residence.insert(part);
300 apf::MeshEntity* vert = globalToVert[gid];
301 m->setResidence(vert, residence);
308 void constructRemotes(apf::GlobalToVert& globalToVert)
310 int self = PCU_Comm_Self();
312 APF_ITERATE(apf::GlobalToVert, globalToVert, it) {
314 apf::MeshEntity* vert = it->second;
315 apf::Parts residence;
316 m->getResidence(vert, residence);
317 APF_ITERATE(apf::Parts, residence, rit)
319 PCU_COMM_PACK(*rit, gid);
320 PCU_COMM_PACK(*rit, vert);
324 while (PCU_Comm_Receive()) {
326 PCU_COMM_UNPACK(gid);
327 apf::MeshEntity* remote;
328 PCU_COMM_UNPACK(remote);
329 int from = PCU_Comm_Sender();
330 apf::MeshEntity* vert = globalToVert[gid];
331 m->addRemote(vert, from, remote);
342 AlgParMA(
const RCP<const Environment> &env__,
343 const RCP<
const Comm<int> > &problemComm__,
346 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
349 AlgParMA(
const RCP<const Environment> &env__,
350 const RCP<
const Comm<int> > &problemComm__,
353 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
356 AlgParMA(
const RCP<const Environment> &env__,
357 const RCP<
const Comm<int> > &problemComm__,
360 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
363 AlgParMA(
const RCP<const Environment> &env__,
364 const RCP<
const Comm<int> > &problemComm__,
367 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
371 AlgParMA(
const RCP<const Environment> &env__,
372 const RCP<
const Comm<int> > &problemComm__,
374 env(env__), problemComm(problemComm__), adapter(adapter__)
376 setMPIComm(problemComm__);
382 if (!PCU_Comm_Initialized())
386 PCU_Switch_Comm(mpicomm);
393 else if (adapter->getLocalNumOf(
MESH_FACE)>0)
397 PCU_Max_Ints(&dim,1);
399 throw std::runtime_error(
"ParMA neeeds faces or region information");
402 if (dim!=adapter->getPrimaryEntityType())
403 throw std::runtime_error(
"ParMA only supports balancing primary type==mesh element");
407 gmi_model* g = gmi_load(
".null");
409 m = apf::makeEmptyMdsMesh(g,dim,
false);
414 adapter->getTopologyViewOf(primary_type,tops);
419 const gno_t* element_gids;
420 const part_t* part_ids;
421 adapter->getIDsViewOf(primary_type,element_gids);
422 adapter->getPartsView(part_ids);
423 for (
size_t i =0;i<adapter->getLocalNumOf(primary_type);i++)
424 mapping_elm_gids_index[element_gids[i]] = i;
427 const gno_t* vertex_gids;
431 int c_dim = adapter->getDimension();
432 const scalar_t ** vertex_coords =
new const scalar_t*[c_dim];
433 int* strides =
new int[c_dim];
434 for (
int i=0;i<c_dim;i++)
435 adapter->getCoordinatesViewOf(
MESH_VERTEX,vertex_coords[i],strides[i],i);
439 throw "APF needs adjacency information from elements to vertices";
440 const lno_t* offsets;
441 const gno_t* adjacent_vertex_gids;
442 adapter->getAdjsView(primary_type,
MESH_VERTEX,offsets,adjacent_vertex_gids);
445 apf::GlobalToVert vertex_mapping;
446 apf::ModelEntity* interior = m->findModelEntity(m->getDimension(), 0);
447 for (
size_t i=0;i<adapter->getLocalNumOf(
MESH_VERTEX);i++) {
448 apf::MeshEntity* vtx = m->createVert_(interior);
449 scalar_t temp_coords[3];
450 for (
int k=0;k<c_dim&&k<3;k++)
451 temp_coords[k] = vertex_coords[k][i*strides[k]];
453 for (
int k=c_dim;k<3;k++)
455 apf::Vector3 point(temp_coords[0],temp_coords[1],temp_coords[2]);
456 m->setPoint(vtx,0,point);
457 vertex_mapping[vertex_gids[i]] = vtx;
460 constructElements(adjacent_vertex_gids, adapter->getLocalNumOf(primary_type), offsets, tops, vertex_mapping);
461 constructResidence(vertex_mapping);
462 constructRemotes(vertex_mapping);
469 apf::FieldShape* s = apf::getConstant(dim);
470 gids = apf::createNumbering(m,
"global_ids",s,1);
471 origin_part_ids = apf::createNumbering(m,
"origin",s,1);
474 apf::MeshIterator* itr = m->begin(dim);
475 apf::MeshEntity* ent;
477 while ((ent = m->iterate(itr))) {
478 apf::number(gids,ent,0,0,element_gids[i]);
479 apf::number(origin_part_ids,ent,0,0,PCU_Comm_Self());
485 apf::alignMdsRemotes(m);
486 apf::deriveMdsModel(m);
491 delete [] vertex_coords;
499 template <
typename Adapter>
505 std::string alg_name =
"VtxElm";
506 double imbalance = 1.1;
509 int ghost_bridge=m->getDimension()-1;
512 const Teuchos::ParameterList &pl = env->getParameters();
514 const Teuchos::ParameterList &ppl = pl.sublist(
"parma_parameters");
515 for (ParameterList::ConstIterator iter = ppl.begin();
516 iter != ppl.end(); iter++) {
517 const std::string &zname = pl.name(iter);
518 if (zname ==
"parma_method") {
519 std::string zval = pl.entry(iter).getValue(&zval);
522 else if (zname ==
"step_size") {
523 double zval = pl.entry(iter).getValue(&zval);
526 else if (zname==
"ghost_layers" || zname==
"ghost_bridge") {
527 int zval = pl.entry(iter).getValue(&zval);
528 if (zname==
"ghost_layers")
535 catch (std::exception &e) {
539 const Teuchos::ParameterEntry *pe2 = pl.getEntryPtr(
"imbalance_tolerance");
541 imbalance = pe2->getValue<
double>(&imbalance);
545 bool weightVertex,weightEdge,weightElement;
546 weightVertex=weightEdge=weightElement=
false;
549 apf::Balancer* balancer;
550 const int verbose = 1;
551 if (alg_name==
"Vertex") {
552 balancer = Parma_MakeVtxBalancer(m, step, verbose);
555 else if (alg_name==
"Element") {
556 balancer = Parma_MakeElmBalancer(m, step, verbose);
559 else if (alg_name==
"VtxElm") {
560 balancer = Parma_MakeVtxElmBalancer(m,step,verbose);
561 weightVertex = weightElement=
true;
563 else if (alg_name==
"VtxEdgeElm") {
564 balancer = Parma_MakeVtxEdgeElmBalancer(m, step, verbose);
565 weightVertex=weightEdge=weightElement=
true;
567 else if (alg_name==
"Ghost") {
568 balancer = Parma_MakeGhostDiffuser(m, ghost_layers, ghost_bridge, step, verbose);
571 else if (alg_name==
"Shape") {
572 balancer = Parma_MakeShapeOptimizer(m,step,verbose);
575 else if (alg_name==
"Centroid") {
576 balancer = Parma_MakeCentroidDiffuser(m,step,verbose);
580 throw std::runtime_error(
"No such parma method defined");
584 apf::MeshTag*
weights = setWeights(weightVertex,weightEdge,weightElement);
587 balancer->balance(weights, imbalance);
591 int num_local = adapter->getLocalNumOf(adapter->getPrimaryEntityType());
592 ArrayRCP<part_t> partList(
new part_t[num_local], 0, num_local,
true);
596 apf::MeshEntity* ent;
597 apf::MeshIterator* itr = m->begin(m->getDimension());
599 while ((ent=m->iterate(itr))) {
600 if (m->isOwned(ent)) {
601 part_t target_part_id = apf::getNumber(origin_part_ids,ent,0,0);
602 gno_t element_gid = apf::getNumber(gids,ent,0,0);
603 PCU_COMM_PACK(target_part_id,element_gid);
612 while (PCU_Comm_Receive()) {
614 PCU_COMM_UNPACK(global_id);
615 lno_t local_id = mapping_elm_gids_index[global_id];
616 part_t new_part_id = PCU_Comm_Sender();
617 partList[local_id] = new_part_id;
620 solution->setParts(partList);
623 apf::destroyNumbering(gids);
624 apf::destroyNumbering(origin_part_ids);
625 apf::removeTagFromDimension(m, weights, m->getDimension());
626 m->destroyTag(weights);
636 #endif // HAVE_ZOLTAN2_PARMA
IdentifierAdapter defines the interface for identifiers.
#define Z2_FORWARD_EXCEPTIONS
Forward an exception back through call stack.
MatrixAdapter defines the adapter interface for matrices.
GraphAdapter defines the interface for graph-based user data.
AlgParMA(const RCP< const Environment > &env, const RCP< const Comm< int > > &problemComm, const RCP< const BaseAdapter< user_t > > &adapter)
MeshAdapter defines the interface for mesh input.
Defines the PartitioningSolution class.
Adapter::scalar_t scalar_t
A PartitioningSolution is a solution to a partitioning problem.
VectorAdapter defines the interface for vector input.
Algorithm defines the base class for all algorithms.
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals, polyhedrons, tetrahedrons, hexhedrons, prisms, or pyramids.
virtual void partition(const RCP< PartitioningSolution< Adapter > > &solution)
Partitioning method.
Traits class to handle conversions between gno_t/lno_t and TPL data types (e.g., ParMETIS's idx_t...
BaseAdapter defines methods required by all Adapters.
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
A gathering of useful namespace methods.