Nektar++
MeshPartition.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: MeshPartition.cpp
4//
5// For more information, please see: http://www.nektar.info/
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description:
32//
33//
34////////////////////////////////////////////////////////////////////////////////
35
36#ifndef TIXML_USE_STL
37#define TIXML_USE_STL
38#endif
39
42
43#include <iomanip>
44#include <iostream>
45#include <map>
46#include <vector>
47
48#include <tinyxml.h>
49
54
56
57#include <boost/algorithm/string.hpp>
58#include <boost/format.hpp>
59#include <boost/graph/adjacency_iterator.hpp>
60#include <boost/graph/adjacency_list.hpp>
61#include <boost/graph/detail/edge.hpp>
62
64{
65
67{
68 static MeshPartitionFactory instance;
69 return instance;
70}
71
73 LibUtilities::CommSharedPtr comm, int meshDim,
74 std::map<int, MeshEntity> element,
75 CompositeDescriptor compMap)
76 : m_session(session), m_comm(comm), m_dim(meshDim), m_numFields(0),
77 m_elements(element), m_compMap(compMap), m_fieldNameToId(),
78 m_weightingRequired(false), m_weightBnd(false), m_weightDofs(false),
79 m_parallel(false)
80{
81 // leave the meshpartition method of reading expansions and conditions
84
85 for (auto elIt = m_elements.cbegin(); elIt != m_elements.cend();)
86 {
87 if (elIt->second.ghost)
88 {
89 m_ghostElmts[elIt->first] = elIt->second;
90 elIt = m_elements.erase(elIt);
91 }
92 else
93 {
94 ++elIt;
95 }
96 }
97}
98
99void MeshPartition::PartitionMesh(int nParts, bool shared, bool overlapping,
100 [[maybe_unused]] int nLocal)
101{
102 ASSERTL0(m_parallel || m_elements.size() >= nParts,
103 "Too few elements for this many processes.");
104 m_shared = shared;
105
106 ASSERTL0(!m_parallel || shared,
107 "Parallel partitioning requires shared filesystem.");
108
110 {
112 }
113 CreateGraph();
114
115 PartitionGraph(nParts, overlapping);
116}
117
119{
120 // Find the Expansions tag
121 TiXmlElement *expansionTypes = m_session->GetElement("Nektar/Expansions");
122
123 // Use fine-level expansion for mesh partition (Parallel-in-Time)
125
126 // Find the Expansion type
127 TiXmlElement *expansion = expansionTypes->FirstChildElement();
128
129 std::string expType = expansion->Value();
130
131 // Expansiontypes will contain plenty of data, where relevant at this stage
132 // are composite ID(s) that this expansion type describes, nummodes and a
133 // list of fields that this expansion relates to. If this does not exist the
134 // variable is only set to "DefaultVar".
135 if (expType == "E")
136 {
137 while (expansion)
138 {
139 std::vector<unsigned int> composite;
140 std::vector<unsigned int> nummodes;
141 std::vector<std::string> fieldName;
142
143 const char *nModesStr = expansion->Attribute("NUMMODES");
144 ASSERTL0(nModesStr,
145 "NUMMODES was not defined in EXPANSION section of input");
146 std::string numModesStr = nModesStr;
147 bool valid =
148 ParseUtils::GenerateVector(numModesStr.c_str(), nummodes);
149 ASSERTL0(valid, "Unable to correctly parse the number of modes.");
150
151 if (nummodes.size() == 1)
152 {
153 for (int i = 1; i < m_dim; i++)
154 {
155 nummodes.push_back(nummodes[0]);
156 }
157 }
158 ASSERTL0(nummodes.size() == m_dim,
159 "Number of modes should match mesh dimension");
160
161 const char *fStr = expansion->Attribute("FIELDS");
162 if (fStr)
163 {
164 std::string fieldStr = fStr;
165 bool valid =
166 ParseUtils::GenerateVector(fieldStr.c_str(), fieldName);
167 ASSERTL0(valid, "Unable to correctly parse the field string in "
168 "ExpansionTypes.");
169
170 for (int i = 0; i < fieldName.size(); ++i)
171 {
172 if (m_fieldNameToId.count(fieldName[i]) == 0)
173 {
174 int k = m_fieldNameToId.size();
175 m_fieldNameToId[fieldName[i]] = k;
176 m_numFields++;
177 }
178 }
179 }
180 else
181 {
182 fieldName.push_back("DefaultVar");
183 int k = m_fieldNameToId.size();
184
185 if (m_fieldNameToId.count("DefaultVar") == 0)
186 {
187 ASSERTL0(
188 k == 0,
189 "Omitting field variables and explicitly listing "
190 "them in different ExpansionTypes is wrong practise");
191
192 m_fieldNameToId["DefaultVar"] = k;
193 m_numFields++;
194 }
195 }
196
197 std::string compositeStr = expansion->Attribute("COMPOSITE");
198 ASSERTL0(compositeStr.length() > 3,
199 "COMPOSITE must be specified in expansion definition");
200 int beg = compositeStr.find_first_of("[");
201 int end = compositeStr.find_first_of("]");
202 std::string compositeListStr =
203 compositeStr.substr(beg + 1, end - beg - 1);
204 bool parseGood = ParseUtils::GenerateSeqVector(
205 compositeListStr.c_str(), composite);
206 ASSERTL0(parseGood && !composite.empty(),
207 (std::string("Unable to read composite index range: ") +
208 compositeListStr)
209 .c_str());
210
211 // construct mapping (elmt id, field name) -> nummodes
212 for (int i = 0; i < composite.size(); ++i)
213 {
214 auto &shapeType = m_compMap[composite[i]].first;
215 auto &elmtIds = m_compMap[composite[i]].second;
216
217 for (int j = 0; j < fieldName.size(); j++)
218 {
219 for (auto &elid : elmtIds)
220 {
221 m_expansions[elid][fieldName[j]] = nummodes;
222 m_shape[elid] = shapeType;
223 }
224 }
225 }
226
227 expansion = expansion->NextSiblingElement("E");
228 }
229 }
230 else if (expType == "F")
231 {
232 ASSERTL0(expansion->Attribute("FILE"),
233 "Attribute FILE expected for type F expansion");
234 std::string filenameStr = expansion->Attribute("FILE");
235 ASSERTL0(!filenameStr.empty(),
236 "A filename must be specified for the FILE "
237 "attribute of expansion");
238
239 // Create fieldIO object to load file
240 // need a serial communicator to avoid problems with
241 // shared file system
244 std::string iofmt =
245 LibUtilities::FieldIO::GetFileType(filenameStr, comm);
248 iofmt, comm, m_session->GetSharedFilesystem());
249
250 // Load field definitions from file
251 std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs;
252 f->Import(filenameStr, fielddefs);
253
254 // Parse field definitions
255 for (int i = 0; i < fielddefs.size(); ++i)
256 {
257 // Name of fields
258 for (int j = 0; j < fielddefs[i]->m_fields.size(); ++j)
259 {
260 std::string fieldName = fielddefs[i]->m_fields[j];
261 if (m_fieldNameToId.count(fieldName) == 0)
262 {
263 int k = m_fieldNameToId.size();
264 m_fieldNameToId[fieldName] = k;
265 m_numFields++;
266 }
267 }
268 // Number of modes and shape for each element
269 int numHomoDir = fielddefs[i]->m_numHomogeneousDir;
270 int cnt = 0;
271 for (int j = 0; j < fielddefs[i]->m_elementIDs.size(); ++j)
272 {
273 int elid = fielddefs[i]->m_elementIDs[j];
274 std::vector<unsigned int> nummodes;
275 for (int k = 0; k < m_dim; k++)
276 {
277 nummodes.push_back(fielddefs[i]->m_numModes[cnt++]);
278 }
279 if (fielddefs[i]->m_uniOrder)
280 {
281 cnt = 0;
282 }
283 else
284 {
285 cnt += numHomoDir;
286 }
287 for (int k = 0; k < fielddefs[i]->m_fields.size(); k++)
288 {
289 std::string fieldName = fielddefs[i]->m_fields[k];
290 m_expansions[elid][fieldName] = nummodes;
291 }
292 m_shape[elid] = fielddefs[i]->m_shapeType;
293 }
294 }
295 }
296 else
297 {
298 ASSERTL0(false,
299 "Expansion type not defined or not supported at the moment");
300 }
301}
302
303void MeshPartition::PrintPartInfo(std::ostream &out)
304{
305 int nElmt = boost::num_vertices(m_graph);
306 int nPart = m_localPartition.size();
307
308 out << "# Partition information:" << std::endl;
309 out << "# No. elements : " << nElmt << std::endl;
310 out << "# No. partitions: " << nPart << std::endl;
311 out << "# ID nElmt nLocDof nBndDof" << std::endl;
312
313 BoostVertexIterator vertit, vertit_end;
314 std::vector<int> partElmtCount(nPart, 0);
315 std::vector<int> partLocCount(nPart, 0);
316 std::vector<int> partBndCount(nPart, 0);
317
318 std::map<int, int> elmtSizes;
319 std::map<int, int> elmtBndSizes;
320
321 for (std::map<int, NummodesPerField>::iterator expIt = m_expansions.begin();
322 expIt != m_expansions.end(); ++expIt)
323 {
324 int elid = expIt->first;
325 NummodesPerField npf = expIt->second;
326
327 for (NummodesPerField::iterator it = npf.begin(); it != npf.end(); ++it)
328 {
329 ASSERTL0(it->second.size() == m_dim,
330 " Number of directional"
331 " modes in expansion spec for element id = " +
332 std::to_string(elid) + " and field " + it->first +
333 " does not correspond to mesh dimension");
334
335 int na = it->second[0];
336 int nb = 0;
337 int nc = 0;
338 if (m_dim >= 2)
339 {
340 nb = it->second[1];
341 }
342 if (m_dim == 3)
343 {
344 nc = it->second[2];
345 }
346
347 elmtSizes[elid] =
348 CalculateElementWeight(m_shape[elid], false, na, nb, nc);
349 elmtBndSizes[elid] =
350 CalculateElementWeight(m_shape[elid], true, na, nb, nc);
351 }
352 }
353
354 for (boost::tie(vertit, vertit_end) = boost::vertices(m_graph);
355 vertit != vertit_end; ++vertit)
356 {
357 int partId = m_graph[*vertit].partition;
358 partElmtCount[partId]++;
359 partLocCount[partId] += elmtSizes[m_graph[*vertit].id];
360 partBndCount[partId] += elmtBndSizes[m_graph[*vertit].id];
361 }
362
363 for (int i = 0; i < nPart; ++i)
364 {
365 out << i << " " << partElmtCount[i] << " " << partLocCount[i] << " "
366 << partBndCount[i] << std::endl;
367 }
368}
369
371{
372 if (!m_session->DefinesElement("Nektar/Conditions/SolverInfo"))
373 {
374 // No SolverInfo = no change of default action to weight
375 // mesh graph.
376 return;
377 }
378
379 TiXmlElement *solverInfoElement =
380 m_session->GetElement("Nektar/Conditions/SolverInfo");
381
382 // Use fine-level solver info for mesh partition (Parallel-in-Time)
384
385 TiXmlElement *solverInfo = solverInfoElement->FirstChildElement("I");
386 ASSERTL0(solverInfo, "Cannot read SolverInfo tags");
387
388 while (solverInfo)
389 {
390 // read the property name
391 ASSERTL0(solverInfo->Attribute("PROPERTY"),
392 "Missing PROPERTY attribute in solver info "
393 "section. ");
394 std::string solverProperty = solverInfo->Attribute("PROPERTY");
395 ASSERTL0(!solverProperty.empty(),
396 "Solver info properties must have a non-empty "
397 "name. ");
398 // make sure that solver property is capitalised
399 std::string solverPropertyUpper = boost::to_upper_copy(solverProperty);
400
401 // read the value
402 ASSERTL0(solverInfo->Attribute("VALUE"),
403 "Missing VALUE attribute in solver info section. ");
404 std::string solverValue = solverInfo->Attribute("VALUE");
405 ASSERTL0(!solverValue.empty(),
406 "Solver info properties must have a non-empty value");
407 // make sure that property value is capitalised
408 std::string propertyValueUpper = boost::to_upper_copy(solverValue);
409
410 if (solverPropertyUpper == "WEIGHTPARTITIONS")
411 {
412 if (propertyValueUpper == "DOF")
413 {
414 m_weightingRequired = true;
415 m_weightDofs = true;
416 }
417 else if (propertyValueUpper == "BOUNDARY")
418 {
419 m_weightingRequired = true;
420 m_weightBnd = true;
421 }
422 else if (propertyValueUpper == "BOTH")
423 {
424 m_weightingRequired = true;
425 m_weightDofs = true;
426 m_weightBnd = true;
427 }
428 return;
429 }
430 solverInfo = solverInfo->NextSiblingElement("I");
431 }
432}
433
434/*
435 * Calculate element weights based on
436 * - element type (Q,T,H,P,R,A)
437 * - nummodes in expansion which this element belongs to via composite.
438 *
439 * For each element we prepare two vertex weightings, one associated
440 * with the number of matrix elements associated with it (to balance
441 * matrix multiplication work) and another associated
442 * with all work which scales linearly with the number of its
443 * coefficients: communication, vector updates etc.
444 *
445 * \todo Refactor this code to explicitly represent performance model
446 * and flexibly generate graph vertex weights depending on perf data.
447 */
449{
450 std::vector<unsigned int> weight(m_numFields, 1);
451 std::map<int, MeshEntity>::iterator eIt;
452 for (eIt = m_elements.begin(); eIt != m_elements.end(); ++eIt)
453 {
454 m_vertWeights[eIt->second.origId] = weight;
455 m_vertBndWeights[eIt->second.origId] = weight;
456 m_edgeWeights[eIt->second.origId] = weight;
457 }
458
459 for (std::map<int, NummodesPerField>::iterator expIt = m_expansions.begin();
460 expIt != m_expansions.end(); ++expIt)
461 {
462 int elid = expIt->first;
463 NummodesPerField npf = expIt->second;
464
465 for (NummodesPerField::iterator it = npf.begin(); it != npf.end(); ++it)
466 {
467 ASSERTL0(it->second.size() == m_dim,
468 " Number of directional"
469 " modes in expansion spec for element id = " +
470 std::to_string(elid) + " and field " + it->first +
471 " does not correspond to mesh dimension");
472
473 int na = it->second[0];
474 int nb = 0;
475 int nc = 0;
476 if (m_dim >= 2)
477 {
478 nb = it->second[1];
479 }
480 if (m_dim == 3)
481 {
482 nc = it->second[2];
483 }
484
485 // Assume for parallel partitioning that this is just missing from
486 // our partition.
487 if (m_vertWeights.find(elid) == m_vertWeights.end())
488 {
489 continue;
490 }
491
492 m_vertWeights[elid][m_fieldNameToId[it->first]] =
493 CalculateElementWeight(m_shape[elid], false, na, nb, nc);
494 m_vertBndWeights[elid][m_fieldNameToId[it->first]] =
495 CalculateElementWeight(m_shape[elid], true, na, nb, nc);
496 m_edgeWeights[elid][m_fieldNameToId[it->first]] =
497 CalculateEdgeWeight(m_shape[elid], na, nb, nc);
498 }
499 } // for i
500}
501
503{
504 // Maps edge/face to first mesh element id.
505 // On locating second mesh element id, graph edge is created instead.
506 std::unordered_map<int, std::vector<int>> vGraphEdges;
507 int vcnt = 0;
508
509 for (auto &elmt : m_elements)
510 {
511 auto vert = boost::add_vertex(m_graph);
512 m_graph[vert].id = elmt.first;
513 m_graph[vert].partition = 0;
514
516 {
517 m_graph[vert].weight = m_vertWeights[elmt.second.origId];
518 m_graph[vert].bndWeight = m_vertBndWeights[elmt.second.origId];
519 m_graph[vert].edgeWeight = m_edgeWeights[elmt.second.origId];
520 }
521
522 // Process element entries and add graph edges
523 for (auto &eId : elmt.second.list)
524 {
525 // Look to see if we've examined this edge/face before
526 // If so, we've got both graph vertices so add edge
527 auto edgeIt = vGraphEdges.find(eId);
528 if (edgeIt != vGraphEdges.end())
529 {
530 for (auto &iId : edgeIt->second)
531 {
532 BoostEdge e = boost::add_edge(vcnt, iId, m_graph).first;
533 m_graph[e].id = vcnt;
534 }
535 vGraphEdges[eId].push_back(vcnt);
536 }
537 else
538 {
539 std::vector<int> Id;
540 Id.push_back(vcnt);
541 vGraphEdges[eId] = Id;
542 }
543 }
544
545 // Increment counter for graph vertex id.
546 ++vcnt;
547 }
548
549 // Now process ghost elements.
550 for (auto &ghost : m_ghostElmts)
551 {
552 auto vert = boost::add_vertex(m_graph);
553 m_graph[vert].id = ghost.first;
554 m_graph[vert].partition = -1;
555
556 for (auto &facet : ghost.second.list)
557 {
558 auto edgeIt = vGraphEdges.find(facet);
559 if (edgeIt != vGraphEdges.end())
560 {
561 for (auto &iId : edgeIt->second)
562 {
563 BoostEdge e = boost::add_edge(vcnt, iId, m_graph).first;
564 m_graph[e].id = vcnt;
565 }
566 vGraphEdges[facet].push_back(vcnt);
567 }
568 }
569
570 // Increment counter for graph vertex id.
571 ++vcnt;
572 }
573}
574
575/**
576 * @brief Partition the graph.
577 *
578 * This routine partitions the graph @p pGraph into @p nParts, producing
579 * subgraphs that are populated in @p pLocalPartition. If the @p
580 * overlapping option is set (which is used for post-processing
581 * purposes), the resulting partitions are extended to cover
582 * neighbouring elements by additional vertex on the dual graph, which
583 * produces overlapping partitions (i.e. the intersection of two
584 * connected partitions is non-empty).
585 *
586 * @param nParts Number of partitions.
587 * @param pLocalPartition Vector of sub-graphs representing each
588 * @param overlapping True if resulting partitions should overlap.
589 */
590void MeshPartition::PartitionGraph(int nParts, bool overlapping)
591{
592 int i;
593 int nGraphVerts = boost::num_vertices(m_graph);
594 int nGhost = m_ghostElmts.size();
595 int nLocal = nGraphVerts - nGhost;
596
597 int ncon = 1;
599 {
600 ncon = 2;
601 }
602
603 // Convert boost graph into CSR format
604 BoostVertexIterator vertit, vertit_end;
605 BoostAdjacencyIterator adjvertit, adjvertit_end;
606 Array<OneD, int> part(nGraphVerts, 0);
607
608 if (m_comm->GetRowComm()->TreatAsRankZero() || m_parallel)
609 {
610 int acnt = 0;
611 int vcnt = 0;
612 int nWeight = ncon * nLocal;
613
614 Array<OneD, int> xadj(nLocal + 1);
615 std::vector<int> adjncy_tmp, adjwgt_tmp;
616 Array<OneD, int> vwgt(nWeight, 1);
617 Array<OneD, int> vsize(nLocal, 1);
618
619 // Initialise starting point of adjacency array.
620 xadj[0] = 0;
621
622 for (boost::tie(vertit, vertit_end) = boost::vertices(m_graph);
623 vertit != vertit_end && vcnt < nLocal; ++vertit)
624 {
625 for (boost::tie(adjvertit, adjvertit_end) =
626 boost::adjacent_vertices(*vertit, m_graph);
627 adjvertit != adjvertit_end; ++adjvertit, ++acnt)
628 {
629 adjncy_tmp.push_back(m_graph[*adjvertit].id);
631 {
632 adjwgt_tmp.push_back(m_graph[*vertit].edgeWeight[0]);
633 }
634 else
635 {
636 adjwgt_tmp.push_back(1);
637 }
638 }
639
640 xadj[++vcnt] = acnt;
641
643 {
644 int ccnt = 0;
645 if (m_weightDofs)
646 {
647 vwgt[ncon * (vcnt - 1) + ccnt] = m_graph[*vertit].weight[0];
648 ccnt++;
649 }
650 if (m_weightBnd)
651 {
652 vwgt[ncon * (vcnt - 1) + ccnt] =
653 m_graph[*vertit].bndWeight[0];
654 }
655 }
656 }
657
658 Array<OneD, int> adjncy(adjncy_tmp.size(), &adjncy_tmp[0]);
659 Array<OneD, int> adjwgt(adjwgt_tmp.size(), &adjwgt_tmp[0]);
660
661 // Call partitioner to partition graph
662 int vol = 0;
663
664 try
665 {
666 //////////////////////////////////////////////////////
667 // On a cartesian communicator do mesh partition just on the first
668 // column
669 // so there is no doubt the partitions are all the same in all the
670 // columns
671 if (m_comm->GetColumnComm()->GetRank() == 0)
672 {
673 // Attempt partitioning.
674 v_PartitionGraphImpl(nLocal, ncon, xadj, adjncy, vwgt, vsize,
675 adjwgt, nParts, vol, part);
676
677 // Check the partitioner produced a valid partition and fix if
678 // not.
679 if (!m_parallel)
680 {
681 CheckPartitions(nParts, part);
682 }
683
684 if (!m_shared)
685 {
686 // distribute among columns
687 for (i = 1; i < m_comm->GetColumnComm()->GetSize(); ++i)
688 {
689 m_comm->GetColumnComm()->Send(i, part);
690 }
691 }
692 }
693 else
694 {
695 m_comm->GetColumnComm()->Recv(0, part);
696 }
697 }
698 catch (...)
699 {
700 NEKERROR(ErrorUtil::efatal, "Error in calling graph partitioner.");
701 }
702 }
703 else if (!m_parallel)
704 {
705 m_comm->GetRowComm()->Recv(0, part);
706 }
707
708 // Create storage for this (and possibly other) process's partitions.
709 i = 0;
710
711 if (!m_parallel)
712 {
713 // Populate subgraph(s) for each rank.
714 for (boost::tie(vertit, vertit_end) = boost::vertices(m_graph);
715 vertit != vertit_end; ++vertit, ++i)
716 {
717 m_localPartition[part[i]].push_back(m_graph[*vertit].id);
718 }
719 }
720 else
721 {
722 // Figure out how many vertices we're going to get from each processor.
723 int nproc = m_comm->GetSpaceComm()->GetSize(),
724 rank = m_comm->GetSpaceComm()->GetRank();
725 std::vector<int> numToSend(nproc, 0), numToRecv(nproc);
726 std::map<int, std::vector<int>> procMap;
727
728 for (boost::tie(vertit, vertit_end) = boost::vertices(m_graph);
729 vertit != vertit_end && i < nLocal; ++vertit, ++i)
730 {
731 int toProc = part[i];
732 numToSend[toProc]++;
733 procMap[toProc].push_back(m_graph[*vertit].id);
734 }
735
736 m_comm->GetSpaceComm()->AlltoAll(numToSend, numToRecv);
737
738 // Build offsets for all-to-all communication
739 std::vector<int> sendOffsetMap(nproc), recvOffsetMap(nproc);
740
741 sendOffsetMap[0] = 0;
742 recvOffsetMap[0] = 0;
743 for (int i = 1; i < nproc; ++i)
744 {
745 sendOffsetMap[i] = sendOffsetMap[i - 1] + numToSend[i - 1];
746 recvOffsetMap[i] = recvOffsetMap[i - 1] + numToRecv[i - 1];
747 }
748
749 // Build data to send
750 int totalSend = Vmath::Vsum(nproc, &numToSend[0], 1);
751 int totalRecv = Vmath::Vsum(nproc, &numToRecv[0], 1);
752
753 std::vector<int> sendData(totalSend), recvData(totalRecv);
754
755 int cnt = 0;
756 for (auto &verts : procMap)
757 {
758 for (auto &vert : verts.second)
759 {
760 sendData[cnt++] = vert;
761 }
762 }
763
764 // Send ID map to processors
765 m_comm->GetSpaceComm()->AlltoAllv(sendData, numToSend, sendOffsetMap,
766 recvData, numToRecv, recvOffsetMap);
767
768 // Finally, populate m_localPartition for this processor. Could contain
769 // duplicates so erase those first.
770 std::unordered_set<int> uniqueIDs;
771 for (auto &id : recvData)
772 {
773 uniqueIDs.insert(id);
774 }
775
776 m_localPartition[rank].insert(m_localPartition[rank].begin(),
777 uniqueIDs.begin(), uniqueIDs.end());
778 }
779
780 // If the overlapping option is set (for post-processing purposes),
781 // add vertices that correspond to the neighbouring elements.
782 if (overlapping)
783 {
784 ASSERTL0(!m_parallel, "Overlapping partitioning not supported in "
785 "parallel execution");
786
787 for (boost::tie(vertit, vertit_end) = boost::vertices(m_graph);
788 vertit != vertit_end; ++vertit)
789 {
790 for (boost::tie(adjvertit, adjvertit_end) =
791 boost::adjacent_vertices(*vertit, m_graph);
792 adjvertit != adjvertit_end; ++adjvertit)
793 {
794 if (part[*adjvertit] != part[*vertit])
795 {
796 m_localPartition[part[*vertit]].push_back(
797 m_graph[*adjvertit].id);
798 }
799 }
800 }
801 }
802}
803
805{
806 unsigned int i = 0;
807 unsigned int cnt = 0;
808 bool valid = true;
809
810 // Check that every process has at least one element assigned
811 for (i = 0; i < nParts; ++i)
812 {
813 cnt = std::count(pPart.begin(), pPart.end(), i);
814 if (cnt == 0)
815 {
816 valid = false;
817 }
818 }
819
820 // If the graph partitioner produced an invalid partition, repartition
821 // naively. Elements are assigned to processes in a round-robin fashion.
822 // It is assumed that graph partitioner failure only occurs when the number
823 // of elements is approx. the number of processes, so this approach should
824 // not be too inefficient communication-wise.
825 if (!valid)
826 {
827 for (i = 0; i < pPart.size(); ++i)
828 {
829 pPart[i] = i % nParts;
830 }
831 }
832}
833
834void MeshPartition::GetElementIDs(const int procid,
835 std::vector<unsigned int> &elmtid)
836{
837 BoostVertexIterator vertit, vertit_end;
838
839 auto it = m_localPartition.find(procid);
840
841 ASSERTL0(it != m_localPartition.end(), "Unable to find local partition");
842
843 elmtid = m_localPartition[procid];
844}
845
847 bool bndWeight, int na, int nb,
848 int nc)
849{
850 int weight = 0;
851
852 switch (elmtType)
853 {
855 weight = bndWeight
857 na, nb, nc)
859 na, nb, nc);
860 break;
862 weight =
863 bndWeight
865 na, nb, nc)
867 na, nb, nc);
868 break;
870 weight = bndWeight
872 na, nb, nc)
874 na, nb, nc);
875 break;
877 weight = bndWeight
879 na, nb, nc)
881 na, nb, nc);
882 break;
884 weight =
885 bndWeight
887 nb)
889 nb);
890 break;
892 weight =
893 bndWeight
895 nb)
897 break;
899 weight =
900 bndWeight
903 break;
905 weight = 1;
906 break;
907 default:
908 break;
909 }
910
911 return weight;
912}
913
914/**
915 * Calculate the number of modes needed for communication when
916 * in partition boundary, to be used as weighting for edges.
917 * Since we do not know exactly which face this refers to, assume
918 * the max order and quad face (for prisms) as arbitrary choices
919 */
921 int nb, int nc)
922{
923 int weight = 0;
924 int n = std::max(na, std::max(nb, nc));
925 switch (elmtType)
926 {
929 break;
932 break;
935 break;
938 break;
941 weight = n;
942 break;
944 weight = 1;
945 break;
946 default:
947 break;
948 }
949
950 return weight;
951}
952} // namespace Nektar::SpatialDomains
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#define SPATIAL_DOMAINS_EXPORT
static const std::string GetFileType(const std::string &filename, CommSharedPtr comm)
Determine file type of given input file.
Definition: FieldIO.cpp:94
Provides a generic Factory class.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static void GetXMLElementTimeLevel(TiXmlElement *&element, const size_t timeLevel, const bool enableCheck=true)
Get XML elment time level (Parallel-in-Time)
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
Definition: ParseUtils.cpp:130
static bool GenerateSeqVector(const std::string &str, std::vector< unsigned int > &out)
Takes a comma-separated compressed string and converts it to entries in a vector.
Definition: ParseUtils.cpp:104
boost::graph_traits< BoostGraph >::vertex_iterator BoostVertexIterator
virtual void v_PartitionGraphImpl(int &nVerts, int &nVertConds, Nektar::Array< Nektar::OneD, int > &xadj, Nektar::Array< Nektar::OneD, int > &adjcy, Nektar::Array< Nektar::OneD, int > &vertWgt, Nektar::Array< Nektar::OneD, int > &vertSize, Nektar::Array< Nektar::OneD, int > &edgeWgt, int &nparts, int &volume, Nektar::Array< Nektar::OneD, int > &part)=0
void PartitionGraph(int nParts, bool overlapping=false)
Partition the graph.
std::map< int, NummodesPerField > m_expansions
boost::graph_traits< BoostGraph >::adjacency_iterator BoostAdjacencyIterator
std::map< int, MultiWeight > m_edgeWeights
std::map< std::string, int > m_fieldNameToId
std::map< int, MeshEntity > m_ghostElmts
LibUtilities::SessionReaderSharedPtr m_session
std::map< int, MultiWeight > m_vertBndWeights
std::map< int, MultiWeight > m_vertWeights
void PrintPartInfo(std::ostream &out)
void GetElementIDs(const int procid, std::vector< unsigned int > &tmp)
boost::graph_traits< BoostGraph >::edge_descriptor BoostEdge
std::map< int, std::vector< unsigned int > > m_localPartition
int CalculateElementWeight(LibUtilities::ShapeType elmtType, bool bndWeight, int na, int nb, int nc)
std::map< std::string, NumModes > NummodesPerField
int CalculateEdgeWeight(LibUtilities::ShapeType elmtType, int na, int nb, int nc)
void PartitionMesh(int nParts, bool shared=false, bool overlapping=false, int nLocal=0)
std::map< int, MeshEntity > m_elements
MeshPartition(const LibUtilities::SessionReaderSharedPtr session, LibUtilities::CommSharedPtr comm, int meshDim, std::map< int, MeshEntity > element, CompositeDescriptor compMap)
LibUtilities::CommSharedPtr m_comm
void CheckPartitions(int nParts, Array< OneD, int > &pPart)
std::map< int, LibUtilities::ShapeType > m_shape
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:153
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:161
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:279
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:290
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:232
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:259
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:133
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:143
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:210
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:187
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:109
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:121
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:322
std::shared_ptr< SessionReader > SessionReaderSharedPtr
FieldIOFactory & GetFieldIOFactory()
Returns the FieldIO factory.
Definition: FieldIO.cpp:69
CommFactory & GetCommFactory()
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
std::map< int, std::pair< LibUtilities::ShapeType, std::vector< int > > > CompositeDescriptor
Definition: MeshGraph.h:61
MeshPartitionFactory & GetMeshPartitionFactory()
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608