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;
286 regionType = currentLevel.
Get<std::string>(
"aggregationRegionType",
NoFactory::get());
289 regionType = Get< std::string >(currentLevel,
"aggregationRegionType");
292 int numDimensions = 0;
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);
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->getLocalNumElements()
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") {
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;
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);