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");
122 #undef SET_VALID_ENTRY 126 validParamList->set< RCP<const FactoryBase> >(
"Graph", null,
"Generating factory of the graph");
127 validParamList->set< RCP<const FactoryBase> >(
"DofsPerNode", null,
"Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
129 validParamList->set<std::string> (
"OnePt aggregate map name",
"",
130 "Name of input map for single node aggregates. (default='')");
131 validParamList->set<std::string> (
"OnePt aggregate map factory",
"",
132 "Generating factory of (DOF) map for single node aggregates.");
135 validParamList->set<std::string> (
"Interface aggregate map name",
"",
136 "Name of input map for interface aggregates. (default='')");
137 validParamList->set<std::string> (
"Interface aggregate map factory",
"",
138 "Generating factory of (DOF) map for interface aggregates.");
139 validParamList->set<RCP<const FactoryBase> > (
"interfacesDimensions", Teuchos::null,
140 "Describes the dimensions of all the interfaces on this rank.");
141 validParamList->set<RCP<const FactoryBase> > (
"nodeOnInterface", Teuchos::null,
142 "List the LIDs of the nodes on any interface.");
146 validParamList->set<RCP<const FactoryBase> >(
"numDimensions", Teuchos::null,
147 "Number of spatial dimension provided by CoordinatesTransferFactory.");
148 validParamList->set<RCP<const FactoryBase> >(
"lNodesPerDim", Teuchos::null,
149 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
153 validParamList->set<RCP<const FactoryBase> > (
"aggregationRegionType", Teuchos::null,
154 "Type of aggregation to use on the region (\"structured\" or \"uncoupled\")");
156 return validParamList;
159 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
162 Input(currentLevel,
"Graph");
164 ParameterList pL = GetParameterList();
177 "Aggregation region type was not provided by the user!");
184 "numDimensions was not provided by the user on level0!");
191 "lNodesPerDim was not provided by the user on level0!");
194 Input(currentLevel,
"aggregationRegionType");
195 Input(currentLevel,
"numDimensions");
196 Input(currentLevel,
"lNodesPerDim");
202 Input(currentLevel,
"DofsPerNode");
205 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true){
212 "interfacesDimensions was not provided by the user on level0!");
219 "nodeOnInterface was not provided by the user on level0!");
222 Input(currentLevel,
"interfacesDimensions");
223 Input(currentLevel,
"nodeOnInterface");
228 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
229 if (mapOnePtName.length() > 0) {
230 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
231 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
234 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
235 currentLevel.
DeclareInput(mapOnePtName, mapOnePtFact.get());
240 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
245 RCP<Teuchos::FancyOStream> out;
246 if(
const char* dbg = std::getenv(
"MUELU_HYBRIDAGGREGATION_DEBUG")) {
247 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
248 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
250 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
253 *out <<
"Entering hybrid aggregation" << std::endl;
255 ParameterList pL = GetParameterList();
256 bDefinitionPhase_ =
false;
258 if (pL.get<
int>(
"aggregation: max agg size") == -1)
259 pL.set(
"aggregation: max agg size", INT_MAX);
262 RCP<const FactoryBase> graphFact = GetFactory(
"Graph");
265 RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel,
"Graph");
266 RCP<const Map> fineMap = graph->GetDomainMap();
267 const int myRank = fineMap->getComm()->getRank();
268 const int numRanks = fineMap->getComm()->getSize();
270 out->setProcRankAndSize(graph->GetImportMap()->getComm()->getRank(),
271 graph->GetImportMap()->getComm()->getSize());
274 RCP<Aggregates> aggregates = rcp(
new Aggregates(*graph));
275 aggregates->setObjectLabel(
"HB");
278 const LO numRows = graph->GetNodeNumVertices();
279 std::vector<unsigned> aggStat(numRows,
READY);
282 std::string regionType;
283 if(currentLevel.GetLevelID() == 0) {
285 regionType = currentLevel.Get<std::string>(
"aggregationRegionType",
NoFactory::get());
288 regionType = Get< std::string >(currentLevel,
"aggregationRegionType");
291 int numDimensions = 0;
292 if(currentLevel.GetLevelID() == 0) {
294 numDimensions = currentLevel.Get<
int>(
"numDimensions",
NoFactory::get());
297 numDimensions = Get<int>(currentLevel,
"numDimensions");
301 std::string coarseningRate = pL.get<std::string>(
"aggregation: coarsening rate");
302 Teuchos::Array<LO> coarseRate;
304 coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
305 }
catch(
const Teuchos::InvalidArrayStringRepresentation& e) {
306 GetOStream(
Errors,-1) <<
" *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** " 310 TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
312 "\"aggregation: coarsening rate\" must have at least as many" 313 " components as the number of spatial dimensions in the problem.");
316 LO numNonAggregatedNodes = numRows;
317 if (regionType ==
"structured") {
323 const int interpolationOrder = pL.get<
int>(
"aggregation: coarsening order");
324 Array<LO> lFineNodesPerDir(3);
325 if(currentLevel.GetLevelID() == 0) {
327 lFineNodesPerDir = currentLevel.Get<Array<LO> >(
"lNodesPerDim",
NoFactory::get());
330 lFineNodesPerDir = Get<Array<LO> >(currentLevel,
"lNodesPerDim");
334 for(
int dim = numDimensions; dim < 3; ++dim) {
335 lFineNodesPerDir[dim] = 1;
339 RCP<MueLu::IndexManager<LO,GO,NO> > geoData;
350 TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
351 !=
static_cast<size_t>(geoData->getNumLocalFineNodes()),
353 "The local number of elements in the graph's map is not equal to " 354 "the number of nodes given by: lNodesPerDim!");
356 aggregates->SetIndexManager(geoData);
357 aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
359 Set(currentLevel,
"lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
363 if (regionType ==
"uncoupled"){
367 if (pL.get<
bool>(
"aggregation: allow user-specified singletons") ==
true) algos_.push_back(rcp(
new OnePtAggregationAlgorithm (graphFact)));
373 *out <<
" Build interface aggregates" << std::endl;
375 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true) {
376 BuildInterfaceAggregates(currentLevel, aggregates, aggStat, numNonAggregatedNodes,
380 *out <<
"Treat Dirichlet BC" << std::endl;
382 ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
383 if (dirichletBoundaryMap != Teuchos::null)
384 for (LO i = 0; i < numRows; i++)
385 if (dirichletBoundaryMap[i] ==
true)
389 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
390 RCP<Map> OnePtMap = Teuchos::null;
391 if (mapOnePtName.length()) {
392 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
393 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
394 OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName,
NoFactory::get());
396 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
397 OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
401 LO nDofsPerNode = Get<LO>(currentLevel,
"DofsPerNode");
402 GO indexBase = graph->GetDomainMap()->getIndexBase();
403 if (OnePtMap != Teuchos::null) {
404 for (LO i = 0; i < numRows; i++) {
406 GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
407 for (LO kr = 0; kr < nDofsPerNode; kr++)
408 if (OnePtMap->isNodeGlobalElement(grid + kr))
414 Array<LO> lCoarseNodesPerDir(3,-1);
415 Set(currentLevel,
"lCoarseNodesPerDim", lCoarseNodesPerDir);
418 aggregates->AggregatesCrossProcessors(
false);
420 *out <<
"Run all the algorithms on the local rank" << std::endl;
421 for (
size_t a = 0; a < algos_.size(); a++) {
422 std::string phase = algos_[a]->description();
424 *out << regionType <<
" | Executing phase " << a << std::endl;
426 int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
427 algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
428 algos_[a]->SetProcRankVerbose(oldRank);
429 *out << regionType <<
" | Done Executing phase " << a << std::endl;
432 *out <<
"Compute statistics on aggregates" << std::endl;
433 aggregates->ComputeAggregateSizes(
true);
435 Set(currentLevel,
"Aggregates", aggregates);
436 Set(currentLevel,
"numDimensions", numDimensions);
437 Set(currentLevel,
"aggregationRegionTypeCoarse", regionType);
439 GetOStream(
Statistics1) << aggregates->description() << std::endl;
440 *out <<
"HybridAggregation done!" << std::endl;
443 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
446 std::vector<unsigned>& aggStat, LO& numNonAggregatedNodes,
447 Array<LO> coarseRate)
const {
448 FactoryMonitor m(*
this,
"BuildInterfaceAggregates", currentLevel);
450 RCP<Teuchos::FancyOStream> out;
451 if(
const char* dbg = std::getenv(
"MUELU_HYBRIDAGGREGATION_DEBUG")) {
452 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
453 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
455 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
459 if(coarseRate.size() == 1) {coarseRate.resize(3, coarseRate[0]);}
460 ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
461 ArrayRCP<LO> procWinner = aggregates->GetProcWinner() ->getDataNonConst(0);
462 Array<LO> interfacesDimensions = Get<Array<LO> >(currentLevel,
"interfacesDimensions");
463 Array<LO> nodesOnInterfaces = Get<Array<LO> >(currentLevel,
"nodeOnInterface");
464 const int numInterfaces = interfacesDimensions.size() / 3;
465 const int myRank = aggregates->GetMap()->getComm()->getRank();
468 Array<LO> coarseInterfacesDimensions(interfacesDimensions.size());
469 Array<LO> nodesOnCoarseInterfaces;
471 LO endRate, totalNumCoarseNodes = 0, numCoarseNodes;
472 for(
int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
474 for(
int dim = 0; dim < 3; ++dim) {
475 endRate = (interfacesDimensions[3*interfaceIdx + dim] - 1) % coarseRate[dim];
476 if(interfacesDimensions[3*interfaceIdx + dim] == 1) {
477 coarseInterfacesDimensions[3*interfaceIdx + dim] = 1;
479 coarseInterfacesDimensions[3*interfaceIdx + dim]
480 = (interfacesDimensions[3*interfaceIdx+dim]-1) / coarseRate[dim] + 2;
481 if(endRate==0){ coarseInterfacesDimensions[3*interfaceIdx + dim]--;}
483 numCoarseNodes *= coarseInterfacesDimensions[3*interfaceIdx + dim];
485 totalNumCoarseNodes += numCoarseNodes;
487 nodesOnCoarseInterfaces.resize(totalNumCoarseNodes, -1);
490 Array<LO> endRate(3);
491 LO interfaceOffset = 0, aggregateCount = 0, coarseNodeCount = 0;
492 for(
int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
493 ArrayView<LO> fineNodesPerDim = interfacesDimensions(3*interfaceIdx, 3);
494 ArrayView<LO> coarseNodesPerDim = coarseInterfacesDimensions(3*interfaceIdx, 3);
495 LO numInterfaceNodes = 1, numCoarseNodes = 1;
496 for(
int dim = 0; dim < 3; ++dim) {
497 numInterfaceNodes *= fineNodesPerDim[dim];
498 numCoarseNodes *= coarseNodesPerDim[dim];
499 endRate[dim] = fineNodesPerDim[dim] % coarseRate[dim];
501 ArrayView<LO> interfaceNodes = nodesOnInterfaces(interfaceOffset, numInterfaceNodes);
503 interfaceOffset += numInterfaceNodes;
505 LO rem, rate, fineNodeIdx;
506 Array<LO> nodeIJK(3), coarseIJK(3), rootIJK(3);
509 for(LO coarseNodeIdx = 0; coarseNodeIdx < numCoarseNodes; ++coarseNodeIdx) {
510 coarseIJK[2] = coarseNodeIdx / (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
511 rem = coarseNodeIdx % (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
512 coarseIJK[1] = rem / coarseNodesPerDim[0];
513 coarseIJK[0] = rem % coarseNodesPerDim[0];
515 for(LO dim = 0; dim < 3; ++dim) {
516 if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
517 nodeIJK[dim] = fineNodesPerDim[dim] - 1;
519 nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
522 fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
524 if(aggStat[interfaceNodes[fineNodeIdx]] ==
READY) {
525 vertex2AggId[interfaceNodes[fineNodeIdx]] = aggregateCount;
526 procWinner[interfaceNodes[fineNodeIdx]] = myRank;
527 aggStat[interfaceNodes[fineNodeIdx]] =
AGGREGATED;
529 --numNonAggregatedNodes;
531 nodesOnCoarseInterfaces[coarseNodeCount] = vertex2AggId[interfaceNodes[fineNodeIdx]];
538 for(LO nodeIdx = 0; nodeIdx < numInterfaceNodes; ++nodeIdx) {
541 if(aggStat[interfaceNodes[nodeIdx]] ==
AGGREGATED) {
continue;}
543 nodeIJK[2] = nodeIdx / (fineNodesPerDim[0]*fineNodesPerDim[1]);
544 rem = nodeIdx % (fineNodesPerDim[0]*fineNodesPerDim[1]);
545 nodeIJK[1] = rem / fineNodesPerDim[0];
546 nodeIJK[0] = rem % fineNodesPerDim[0];
548 for(
int dim = 0; dim < 3; ++dim) {
549 coarseIJK[dim] = nodeIJK[dim] / coarseRate[dim];
550 rem = nodeIJK[dim] % coarseRate[dim];
551 if(nodeIJK[dim] < fineNodesPerDim[dim] - endRate[dim]) {
552 rate = coarseRate[dim];
556 if(rem > (rate / 2)) {++coarseIJK[dim];}
557 if(coarseNodesPerDim[dim] - coarseIJK[dim] > fineNodesPerDim[dim]-nodeIJK[dim]){
562 for(LO dim = 0; dim < 3; ++dim) {
563 if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
564 nodeIJK[dim] = fineNodesPerDim[dim] - 1;
566 nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
569 fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];
571 vertex2AggId[interfaceNodes[nodeIdx]] = vertex2AggId[interfaceNodes[fineNodeIdx]];
572 procWinner[interfaceNodes[nodeIdx]] = myRank;
573 aggStat[interfaceNodes[nodeIdx]] =
AGGREGATED;
574 --numNonAggregatedNodes;
579 aggregates->SetNumAggregates(aggregateCount);
582 Set(currentLevel,
"coarseInterfacesDimensions", coarseInterfacesDimensions);
583 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()