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>
116 class AlgParMA :
public Algorithm<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;
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 face,
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(
MESH_FACE));
207 num_ws = std::max(num_ws,adapter->getNumWeightsPerOf(entityAPFtoZ2(m->getDimension())));
208 apf::MeshTag* tag = m->createDoubleTag(
"parma_weight",num_ws);
210 setEntWeights(0,tag);
212 setEntWeights(1,tag);
214 setEntWeights(2,tag);
216 setEntWeights(m->getDimension(),tag);
223 void constructElements(
const gno_t* conn,
lno_t nelem,
const lno_t* offsets,
226 apf::ModelEntity* interior = m->findModelEntity(m->getDimension(), 0);
227 for (
lno_t i = 0; i < nelem; ++i) {
228 apf::Mesh::Type etype = topologyZ2toAPF(tops[i]);
230 for (
int j = offsets[i]; j < offsets[i+1]; ++j)
231 verts[j-offsets[i]] = globalToVert[conn[j]];
232 buildElement(m, interior, etype, verts);
235 int getMax(
const apf::GlobalToVert& globalToVert)
238 APF_CONST_ITERATE(apf::GlobalToVert, globalToVert, it)
239 max = std::max(max, it->first);
240 PCU_Max_Ints(&max, 1);
243 void constructResidence(apf::GlobalToVert& globalToVert)
245 int max = getMax(globalToVert);
247 int peers = PCU_Comm_Peers();
248 int quotient = total / peers;
249 int remainder = total % peers;
250 int mySize = quotient;
251 int self = PCU_Comm_Self();
252 if (
self == (peers - 1))
254 typedef std::vector< std::vector<int> > TmpParts;
255 TmpParts tmpParts(mySize);
259 APF_ITERATE(apf::GlobalToVert, globalToVert, it) {
261 int to = std::min(peers - 1, gid / quotient);
262 PCU_COMM_PACK(to, gid);
265 int myOffset =
self * quotient;
268 while (PCU_Comm_Receive()) {
270 PCU_COMM_UNPACK(gid);
271 int from = PCU_Comm_Sender();
272 tmpParts.at(gid - myOffset).push_back(from);
277 for (
int i = 0; i < mySize; ++i) {
278 std::vector<int>& parts = tmpParts[i];
279 for (
size_t j = 0; j < parts.size(); ++j) {
281 int gid = i + myOffset;
282 int nparts = parts.size();
283 PCU_COMM_PACK(to, gid);
284 PCU_COMM_PACK(to, nparts);
285 for (
size_t k = 0; k < parts.size(); ++k)
286 PCU_COMM_PACK(to, parts[k]);
293 while (PCU_Comm_Receive()) {
295 PCU_COMM_UNPACK(gid);
297 PCU_COMM_UNPACK(nparts);
298 apf::Parts residence;
299 for (
int i = 0; i < nparts; ++i) {
301 PCU_COMM_UNPACK(part);
302 residence.insert(part);
304 apf::MeshEntity* vert = globalToVert[gid];
305 m->setResidence(vert, residence);
312 void constructRemotes(apf::GlobalToVert& globalToVert)
314 int self = PCU_Comm_Self();
316 APF_ITERATE(apf::GlobalToVert, globalToVert, it) {
318 apf::MeshEntity* vert = it->second;
319 apf::Parts residence;
320 m->getResidence(vert, residence);
321 APF_ITERATE(apf::Parts, residence, rit)
323 PCU_COMM_PACK(*rit, gid);
324 PCU_COMM_PACK(*rit, vert);
328 while (PCU_Comm_Receive()) {
330 PCU_COMM_UNPACK(gid);
331 apf::MeshEntity* remote;
332 PCU_COMM_UNPACK(remote);
333 int from = PCU_Comm_Sender();
334 apf::MeshEntity* vert = globalToVert[gid];
335 m->addRemote(vert, from, remote);
346 AlgParMA(
const RCP<const Environment> &env__,
347 const RCP<
const Comm<int> > &problemComm__,
348 const RCP<
const IdentifierAdapter<user_t> > &adapter__)
350 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
353 AlgParMA(
const RCP<const Environment> &env__,
354 const RCP<
const Comm<int> > &problemComm__,
355 const RCP<
const VectorAdapter<user_t> > &adapter__)
357 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
360 AlgParMA(
const RCP<const Environment> &env__,
361 const RCP<
const Comm<int> > &problemComm__,
362 const RCP<
const GraphAdapter<user_t,userCoord_t> > &adapter__)
364 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
367 AlgParMA(
const RCP<const Environment> &env__,
368 const RCP<
const Comm<int> > &problemComm__,
369 const RCP<
const MatrixAdapter<user_t,userCoord_t> > &adapter__)
371 throw std::runtime_error(
"ParMA needs a MeshAdapter but you haven't given it one");
375 AlgParMA(
const RCP<const Environment> &env__,
376 const RCP<
const Comm<int> > &problemComm__,
377 const RCP<
const MeshAdapter<user_t> > &adapter__) :
378 env(env__), problemComm(problemComm__), adapter(adapter__)
380 setMPIComm(problemComm__);
386 if (!PCU_Comm_Initialized())
390 PCU_Switch_Comm(mpicomm);
397 else if (adapter->getLocalNumOf(
MESH_FACE)>0)
401 PCU_Max_Ints(&dim,1);
403 throw std::runtime_error(
"ParMA neeeds faces or region information");
406 if (dim!=adapter->getPrimaryEntityType())
407 throw std::runtime_error(
"ParMA only supports balancing primary type==mesh element");
411 gmi_model* g = gmi_load(
".null");
413 m = apf::makeEmptyMdsMesh(g,dim,
false);
418 adapter->getTopologyViewOf(primary_type,tops);
423 const gno_t* element_gids;
425 adapter->getIDsViewOf(primary_type,element_gids);
426 adapter->getPartsView(part_ids);
427 for (
size_t i =0;i<adapter->getLocalNumOf(primary_type);i++)
428 mapping_elm_gids_index[element_gids[i]] = i;
431 const gno_t* vertex_gids;
435 int c_dim = adapter->getDimension();
437 int* strides =
new int[c_dim];
438 for (
int i=0;i<c_dim;i++)
439 adapter->getCoordinatesViewOf(
MESH_VERTEX,vertex_coords[i],strides[i],i);
443 throw "APF needs adjacency information from elements to vertices";
444 const lno_t* offsets;
445 const gno_t* adjacent_vertex_gids;
446 adapter->getAdjsView(primary_type,
MESH_VERTEX,offsets,adjacent_vertex_gids);
449 apf::GlobalToVert vertex_mapping;
450 apf::ModelEntity* interior = m->findModelEntity(m->getDimension(), 0);
451 for (
size_t i=0;i<adapter->getLocalNumOf(
MESH_VERTEX);i++) {
452 apf::MeshEntity* vtx = m->createVert_(interior);
454 for (
int k=0;k<c_dim&&k<3;k++)
455 temp_coords[k] = vertex_coords[k][i*strides[k]];
457 for (
int k=c_dim;k<3;k++)
459 apf::Vector3 point(temp_coords[0],temp_coords[1],temp_coords[2]);
460 m->setPoint(vtx,0,point);
461 vertex_mapping[vertex_gids[i]] = vtx;
464 constructElements(adjacent_vertex_gids, adapter->getLocalNumOf(primary_type), offsets, tops, vertex_mapping);
465 constructResidence(vertex_mapping);
466 constructRemotes(vertex_mapping);
473 apf::FieldShape* s = apf::getConstant(dim);
474 gids = apf::createNumbering(m,
"global_ids",s,1);
475 origin_part_ids = apf::createNumbering(m,
"origin",s,1);
478 apf::MeshIterator* itr = m->begin(dim);
479 apf::MeshEntity* ent;
481 while ((ent = m->iterate(itr))) {
482 apf::number(gids,ent,0,0,element_gids[i]);
483 apf::number(origin_part_ids,ent,0,0,PCU_Comm_Self());
489 apf::alignMdsRemotes(m);
490 apf::deriveMdsModel(m);
495 delete [] vertex_coords;
498 void partition(
const RCP<PartitioningSolution<Adapter> > &solution);
503 template <
typename Adapter>
505 const RCP<PartitioningSolution<Adapter> > &solution
509 std::string alg_name =
"VtxElm";
510 double imbalance = 1.1;
513 int ghost_bridge=m->getDimension()-1;
516 const Teuchos::ParameterList &pl = env->getParameters();
518 const Teuchos::ParameterList &ppl = pl.sublist(
"parma_parameters");
519 for (ParameterList::ConstIterator iter = ppl.begin();
520 iter != ppl.end(); iter++) {
521 const std::string &zname = pl.name(iter);
522 if (zname ==
"parma_method") {
523 std::string zval = pl.entry(iter).getValue(&zval);
526 else if (zname ==
"step_size") {
527 double zval = pl.entry(iter).getValue(&zval);
530 else if (zname==
"ghost_layers" || zname==
"ghost_bridge") {
531 int zval = pl.entry(iter).getValue(&zval);
532 if (zname==
"ghost_layers")
539 catch (std::exception &e) {
543 const Teuchos::ParameterEntry *pe2 = pl.getEntryPtr(
"imbalance_tolerance");
545 imbalance = pe2->getValue<
double>(&imbalance);
549 bool weightVertex,weightEdge,weightFace,weightElement;
550 weightVertex=weightEdge=weightFace=weightElement=
false;
553 apf::Balancer* balancer;
554 const int verbose = 1;
555 if (alg_name==
"Vertex") {
556 balancer = Parma_MakeVtxBalancer(m, step, verbose);
559 else if (alg_name==
"Element") {
560 balancer = Parma_MakeElmBalancer(m, step, verbose);
563 else if (alg_name==
"VtxElm") {
564 balancer = Parma_MakeVtxElmBalancer(m,step,verbose);
565 weightVertex = weightElement=
true;
567 else if (alg_name==
"VtxEdgeElm") {
568 balancer = Parma_MakeVtxEdgeElmBalancer(m, step, verbose);
569 weightVertex=weightEdge=weightElement=
true;
571 else if (alg_name==
"Ghost") {
572 balancer = Parma_MakeGhostDiffuser(m, ghost_layers, ghost_bridge, step, verbose);
573 weightVertex=weightEdge=weightFace=
true;
574 if (3 == m->getDimension()) {
578 else if (alg_name==
"Shape") {
579 balancer = Parma_MakeShapeOptimizer(m,step,verbose);
582 else if (alg_name==
"Centroid") {
583 balancer = Parma_MakeCentroidDiffuser(m,step,verbose);
587 throw std::runtime_error(
"No such parma method defined");
591 apf::MeshTag*
weights = setWeights(weightVertex,weightEdge,weightFace,weightElement);
594 balancer->balance(
weights, imbalance);
598 int num_local = adapter->getLocalNumOf(adapter->getPrimaryEntityType());
599 ArrayRCP<part_t> partList(
new part_t[num_local], 0, num_local,
true);
603 apf::MeshEntity* ent;
604 apf::MeshIterator* itr = m->begin(m->getDimension());
606 while ((ent=m->iterate(itr))) {
607 if (m->isOwned(ent)) {
608 part_t target_part_id = apf::getNumber(origin_part_ids,ent,0,0);
609 gno_t element_gid = apf::getNumber(gids,ent,0,0);
610 PCU_COMM_PACK(target_part_id,element_gid);
619 while (PCU_Comm_Receive()) {
621 PCU_COMM_UNPACK(global_id);
622 lno_t local_id = mapping_elm_gids_index[global_id];
623 part_t new_part_id = PCU_Comm_Sender();
624 partList[local_id] = new_part_id;
627 solution->setParts(partList);
630 apf::destroyNumbering(gids);
631 apf::destroyNumbering(origin_part_ids);
632 apf::removeTagFromDimension(m,
weights, m->getDimension());
643 #endif // HAVE_ZOLTAN2_PARMA
#define Z2_FORWARD_EXCEPTIONS
Forward an exception back through call stack.
AlgParMA(const RCP< const Environment > &env, const RCP< const Comm< int > > &problemComm, const RCP< const BaseAdapter< user_t > > &adapter)
Defines the PartitioningSolution class.
Adapter::scalar_t scalar_t
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.