46 #ifndef MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP_ 47 #define MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP_ 49 #include <Xpetra_Matrix.hpp> 50 #include <Xpetra_Map.hpp> 51 #include <Xpetra_Vector.hpp> 52 #include <Xpetra_MultiVectorFactory.hpp> 53 #include <Xpetra_VectorFactory.hpp> 58 #include "MueLu_InterfaceAggregationAlgorithm.hpp" 59 #include "MueLu_OnePtAggregationAlgorithm.hpp" 60 #include "MueLu_PreserveDirichletAggregationAlgorithm.hpp" 61 #include "MueLu_IsolatedNodeAggregationAlgorithm.hpp" 63 #include "MueLu_AggregationPhase1Algorithm.hpp" 64 #include "MueLu_AggregationPhase2aAlgorithm.hpp" 65 #include "MueLu_AggregationPhase2bAlgorithm.hpp" 66 #include "MueLu_AggregationPhase3Algorithm.hpp" 69 #include "MueLu_AggregationStructuredAlgorithm.hpp" 70 #include "MueLu_UncoupledIndexManager.hpp" 77 #include "MueLu_Aggregates.hpp" 80 #include "MueLu_Utilities.hpp" 81 #include "MueLu_AmalgamationInfo.hpp" 86 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
91 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 RCP<ParameterList> validParamList = rcp(
new ParameterList());
96 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
97 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 103 validParamList->getEntry(
"aggregation: ordering").setValidator(
104 rcp(
new validatorType(Teuchos::tuple<std::string>(
"natural",
"graph",
"random"),
"aggregation: ordering")));
111 SET_VALID_ENTRY(
"aggregation: error on nodes with no on-rank neighbors");
123 #undef SET_VALID_ENTRY 127 validParamList->set< RCP<const FactoryBase> >(
"Graph", null,
"Generating factory of the graph");
128 validParamList->set< RCP<const FactoryBase> >(
"DofsPerNode", null,
"Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
130 validParamList->set<std::string> (
"OnePt aggregate map name",
"",
131 "Name of input map for single node aggregates. (default='')");
132 validParamList->set<std::string> (
"OnePt aggregate map factory",
"",
133 "Generating factory of (DOF) map for single node aggregates.");
136 validParamList->set<std::string> (
"Interface aggregate map name",
"",
137 "Name of input map for interface aggregates. (default='')");
138 validParamList->set<std::string> (
"Interface aggregate map factory",
"",
139 "Generating factory of (DOF) map for interface aggregates.");
140 validParamList->set<RCP<const FactoryBase> > (
"interfacesDimensions", Teuchos::null,
141 "Describes the dimensions of all the interfaces on this rank.");
142 validParamList->set<RCP<const FactoryBase> > (
"nodeOnInterface", Teuchos::null,
143 "List the LIDs of the nodes on any interface.");
147 validParamList->set<RCP<const FactoryBase> >(
"numDimensions", Teuchos::null,
148 "Number of spatial dimension provided by CoordinatesTransferFactory.");
149 validParamList->set<RCP<const FactoryBase> >(
"lNodesPerDim", Teuchos::null,
150 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
154 validParamList->set<RCP<const FactoryBase> > (
"aggregationRegionType", Teuchos::null,
155 "Type of aggregation to use on the region (\"structured\" or \"uncoupled\")");
157 return validParamList;
160 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
163 Input(currentLevel,
"Graph");
165 ParameterList pL = GetParameterList();
178 "Aggregation region type was not provided by the user!");
185 "numDimensions was not provided by the user on level0!");
192 "lNodesPerDim was not provided by the user on level0!");
195 Input(currentLevel,
"aggregationRegionType");
196 Input(currentLevel,
"numDimensions");
197 Input(currentLevel,
"lNodesPerDim");
203 Input(currentLevel,
"DofsPerNode");
206 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true){
213 "interfacesDimensions was not provided by the user on level0!");
220 "nodeOnInterface was not provided by the user on level0!");
223 Input(currentLevel,
"interfacesDimensions");
224 Input(currentLevel,
"nodeOnInterface");
229 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
230 if (mapOnePtName.length() > 0) {
231 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
232 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
235 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
236 currentLevel.
DeclareInput(mapOnePtName, mapOnePtFact.get());
241 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
246 RCP<Teuchos::FancyOStream> out;
247 if(
const char* dbg = std::getenv(
"MUELU_HYBRIDAGGREGATION_DEBUG")) {
248 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
249 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
251 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
254 *out <<
"Entering hybrid aggregation" << std::endl;
256 ParameterList pL = GetParameterList();
257 bDefinitionPhase_ =
false;
259 if (pL.get<
int>(
"aggregation: max agg size") == -1)
260 pL.set(
"aggregation: max agg size", INT_MAX);
263 RCP<const FactoryBase> graphFact = GetFactory(
"Graph");
266 RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel,
"Graph");
267 RCP<const Map> fineMap = graph->GetDomainMap();
268 const int myRank = fineMap->getComm()->getRank();
269 const int numRanks = fineMap->getComm()->getSize();
271 out->setProcRankAndSize(graph->GetImportMap()->getComm()->getRank(),
272 graph->GetImportMap()->getComm()->getSize());
275 RCP<Aggregates> aggregates = rcp(
new Aggregates(*graph));
276 aggregates->setObjectLabel(
"HB");
279 const LO numRows = graph->GetNodeNumVertices();
280 std::vector<unsigned> aggStat(numRows,
READY);
283 std::string regionType;
284 if(currentLevel.GetLevelID() == 0) {
286 regionType = currentLevel.Get<std::string>(
"aggregationRegionType",
NoFactory::get());
289 regionType = Get< std::string >(currentLevel,
"aggregationRegionType");
292 int numDimensions = 0;
293 if(currentLevel.GetLevelID() == 0) {
295 numDimensions = currentLevel.Get<
int>(
"numDimensions",
NoFactory::get());
298 numDimensions = Get<int>(currentLevel,
"numDimensions");
302 std::string coarseningRate = pL.get<std::string>(
"aggregation: coarsening rate");
303 Teuchos::Array<LO> coarseRate;
305 coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
306 }
catch(
const Teuchos::InvalidArrayStringRepresentation& e) {
307 GetOStream(
Errors,-1) <<
" *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** " 311 TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
313 "\"aggregation: coarsening rate\" must have at least as many" 314 " components as the number of spatial dimensions in the problem.");
317 LO numNonAggregatedNodes = numRows;
318 if (regionType ==
"structured") {
324 const int interpolationOrder = pL.get<
int>(
"aggregation: coarsening order");
325 Array<LO> lFineNodesPerDir(3);
326 if(currentLevel.GetLevelID() == 0) {
328 lFineNodesPerDir = currentLevel.Get<Array<LO> >(
"lNodesPerDim",
NoFactory::get());
331 lFineNodesPerDir = Get<Array<LO> >(currentLevel,
"lNodesPerDim");
335 for(
int dim = numDimensions; dim < 3; ++dim) {
336 lFineNodesPerDir[dim] = 1;
340 RCP<MueLu::IndexManager<LO,GO,NO> > geoData;
351 TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
352 !=
static_cast<size_t>(geoData->getNumLocalFineNodes()),
354 "The local number of elements in the graph's map is not equal to " 355 "the number of nodes given by: lNodesPerDim!");
357 aggregates->SetIndexManager(geoData);
358 aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
360 Set(currentLevel,
"lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
364 if (regionType ==
"uncoupled"){
368 if (pL.get<
bool>(
"aggregation: allow user-specified singletons") ==
true) algos_.push_back(rcp(
new OnePtAggregationAlgorithm (graphFact)));
374 *out <<
" Build interface aggregates" << std::endl;
376 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true) {
377 BuildInterfaceAggregates(currentLevel, aggregates, aggStat, numNonAggregatedNodes,
381 *out <<
"Treat Dirichlet BC" << std::endl;
383 ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
384 if (dirichletBoundaryMap != Teuchos::null)
385 for (LO i = 0; i < numRows; i++)
386 if (dirichletBoundaryMap[i] ==
true)
390 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
391 RCP<Map> OnePtMap = Teuchos::null;
392 if (mapOnePtName.length()) {
393 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
394 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
395 OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName,
NoFactory::get());
397 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
398 OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
402 LO nDofsPerNode = Get<LO>(currentLevel,
"DofsPerNode");
403 GO indexBase = graph->GetDomainMap()->getIndexBase();
404 if (OnePtMap != Teuchos::null) {
405 for (LO i = 0; i < numRows; i++) {
407 GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
408 for (LO kr = 0; kr < nDofsPerNode; kr++)
409 if (OnePtMap->isNodeGlobalElement(grid + kr))
415 Array<LO> lCoarseNodesPerDir(3,-1);
416 Set(currentLevel,
"lCoarseNodesPerDim", lCoarseNodesPerDir);
419 aggregates->AggregatesCrossProcessors(
false);
421 *out <<
"Run all the algorithms on the local rank" << std::endl;
422 for (
size_t a = 0; a < algos_.size(); a++) {
423 std::string phase = algos_[a]->description();
425 *out << regionType <<
" | Executing phase " << a << std::endl;
427 int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
428 algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
429 algos_[a]->SetProcRankVerbose(oldRank);
430 *out << regionType <<
" | Done Executing phase " << a << std::endl;
433 *out <<
"Compute statistics on aggregates" << std::endl;
434 aggregates->ComputeAggregateSizes(
true);
436 Set(currentLevel,
"Aggregates", aggregates);
437 Set(currentLevel,
"numDimensions", numDimensions);
438 Set(currentLevel,
"aggregationRegionTypeCoarse", regionType);
440 GetOStream(
Statistics1) << aggregates->description() << std::endl;
441 *out <<
"HybridAggregation done!" << std::endl;
444 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
447 std::vector<unsigned>& aggStat, LO& numNonAggregatedNodes,
448 Array<LO> coarseRate)
const {
449 FactoryMonitor m(*
this,
"BuildInterfaceAggregates", currentLevel);
451 RCP<Teuchos::FancyOStream> out;
452 if(
const char* dbg = std::getenv(
"MUELU_HYBRIDAGGREGATION_DEBUG")) {
453 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
454 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
456 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
460 if(coarseRate.size() == 1) {coarseRate.resize(3, coarseRate[0]);}
461 ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
462 ArrayRCP<LO> procWinner = aggregates->GetProcWinner() ->getDataNonConst(0);
463 Array<LO> interfacesDimensions = Get<Array<LO> >(currentLevel,
"interfacesDimensions");
464 Array<LO> nodesOnInterfaces = Get<Array<LO> >(currentLevel,
"nodeOnInterface");
465 const int numInterfaces = interfacesDimensions.size() / 3;
466 const int myRank = aggregates->GetMap()->getComm()->getRank();
469 Array<LO> coarseInterfacesDimensions(interfacesDimensions.size());
470 Array<LO> nodesOnCoarseInterfaces;
472 LO endRate, totalNumCoarseNodes = 0, numCoarseNodes;
473 for(
int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
475 for(
int dim = 0; dim < 3; ++dim) {
476 endRate = (interfacesDimensions[3*interfaceIdx + dim] - 1) % coarseRate[dim];
477 if(interfacesDimensions[3*interfaceIdx + dim] == 1) {
478 coarseInterfacesDimensions[3*interfaceIdx + dim] = 1;
480 coarseInterfacesDimensions[3*interfaceIdx + dim]
481 = (interfacesDimensions[3*interfaceIdx+dim]-1) / coarseRate[dim] + 2;
482 if(endRate==0){ coarseInterfacesDimensions[3*interfaceIdx + dim]--;}
484 numCoarseNodes *= coarseInterfacesDimensions[3*interfaceIdx + dim];
486 totalNumCoarseNodes += numCoarseNodes;
488 nodesOnCoarseInterfaces.resize(totalNumCoarseNodes, -1);
491 Array<LO> endRate(3);
492 LO interfaceOffset = 0, aggregateCount = 0, coarseNodeCount = 0;
493 for(
int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
494 ArrayView<LO> fineNodesPerDim = interfacesDimensions(3*interfaceIdx, 3);
495 ArrayView<LO> coarseNodesPerDim = coarseInterfacesDimensions(3*interfaceIdx, 3);
496 LO numInterfaceNodes = 1, numCoarseNodes = 1;
497 for(
int dim = 0; dim < 3; ++dim) {
498 numInterfaceNodes *= fineNodesPerDim[dim];
499 numCoarseNodes *= coarseNodesPerDim[dim];
500 endRate[dim] = (fineNodesPerDim[dim]-1) % coarseRate[dim];
502 ArrayView<LO> interfaceNodes = nodesOnInterfaces(interfaceOffset, numInterfaceNodes);
504 interfaceOffset += numInterfaceNodes;
506 LO rem, rate, fineNodeIdx;
507 Array<LO> nodeIJK(3), coarseIJK(3), rootIJK(3);
510 for(LO coarseNodeIdx = 0; coarseNodeIdx < numCoarseNodes; ++coarseNodeIdx) {
511 coarseIJK[2] = coarseNodeIdx / (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
512 rem = coarseNodeIdx % (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
513 coarseIJK[1] = rem / coarseNodesPerDim[0];
514 coarseIJK[0] = rem % coarseNodesPerDim[0];
516 for(LO dim = 0; dim < 3; ++dim) {
517 if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
518 nodeIJK[dim] = fineNodesPerDim[dim] - 1;
520 nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
523 fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
525 if(aggStat[interfaceNodes[fineNodeIdx]] ==
READY) {
526 vertex2AggId[interfaceNodes[fineNodeIdx]] = aggregateCount;
527 procWinner[interfaceNodes[fineNodeIdx]] = myRank;
528 aggStat[interfaceNodes[fineNodeIdx]] =
AGGREGATED;
530 --numNonAggregatedNodes;
532 nodesOnCoarseInterfaces[coarseNodeCount] = vertex2AggId[interfaceNodes[fineNodeIdx]];
539 for(LO nodeIdx = 0; nodeIdx < numInterfaceNodes; ++nodeIdx) {
542 if(aggStat[interfaceNodes[nodeIdx]] ==
AGGREGATED) {
continue;}
544 nodeIJK[2] = nodeIdx / (fineNodesPerDim[0]*fineNodesPerDim[1]);
545 rem = nodeIdx % (fineNodesPerDim[0]*fineNodesPerDim[1]);
546 nodeIJK[1] = rem / fineNodesPerDim[0];
547 nodeIJK[0] = rem % fineNodesPerDim[0];
549 for(
int dim = 0; dim < 3; ++dim) {
550 coarseIJK[dim] = nodeIJK[dim] / coarseRate[dim];
551 rem = nodeIJK[dim] % coarseRate[dim];
552 if(nodeIJK[dim] < fineNodesPerDim[dim] - endRate[dim]) {
553 rate = coarseRate[dim];
557 if(rem > (rate / 2)) {++coarseIJK[dim];}
560 for(LO dim = 0; dim < 3; ++dim) {
561 if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
562 nodeIJK[dim] = fineNodesPerDim[dim] - 1;
564 nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
567 fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
569 vertex2AggId[interfaceNodes[nodeIdx]] = vertex2AggId[interfaceNodes[fineNodeIdx]];
570 procWinner[interfaceNodes[nodeIdx]] = myRank;
571 aggStat[interfaceNodes[nodeIdx]] =
AGGREGATED;
572 --numNonAggregatedNodes;
577 aggregates->SetNumAggregates(aggregateCount);
580 Set(currentLevel,
"coarseInterfacesDimensions", coarseInterfacesDimensions);
581 Set(currentLevel,
"nodeOnCoarseInterface", nodesOnCoarseInterfaces);
Algorithm for coarsening a graph with uncoupled aggregation. keep special marked nodes as singleton n...
Container class for aggregation information.
Timer to be used in factories. Similar to Monitor but with additional timers.
Algorithm for coarsening a graph with structured aggregation.
Namespace for MueLu classes and methods.
static const NoFactory * get()
int GetLevelID() const
Return level number.
Algorithm for coarsening a graph with uncoupled aggregation. creates aggregates along an interface us...
Builds one-to-one aggregates for all Dirichlet boundary nodes. For some applications this might be ne...
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void DeclareInput(Level ¤tLevel) const
Input.
void Build(Level ¤tLevel) const
Build aggregates.
Among unaggregated points, see if we can make a reasonable size aggregate out of it.IdeaAmong unaggregated points, see if we can make a reasonable size aggregate out of it. We do this by looking at neighbors and seeing how many are unaggregated and on my processor. Loosely, base the number of new aggregates created on the percentage of unaggregated nodes.
void BuildInterfaceAggregates(Level ¤tLevel, RCP< Aggregates > aggregates, std::vector< unsigned > &aggStat, LO &numNonAggregatedNodes, Array< LO > coarseRate) const
Specifically build aggregates along interfaces.
Add leftovers to existing aggregatesIdeaIn phase 2b non-aggregated nodes are added to existing aggreg...
Algorithm for coarsening a graph with uncoupled aggregation.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Exception throws to report errors in the internal logical of the program.
HybridAggregationFactory()
Constructor.
Handle leftover nodes. Try to avoid singleton nodesIdeaIn phase 3 we try to stick unaggregated nodes ...
#define SET_VALID_ENTRY(name)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()