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")));
112 SET_VALID_ENTRY(
"aggregation: error on nodes with no on-rank neighbors");
114 #undef SET_VALID_ENTRY 118 validParamList->set< RCP<const FactoryBase> >(
"Graph", null,
"Generating factory of the graph");
119 validParamList->set< RCP<const FactoryBase> >(
"DofsPerNode", null,
"Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
121 validParamList->set<std::string> (
"OnePt aggregate map name",
"",
122 "Name of input map for single node aggregates. (default='')");
123 validParamList->set<std::string> (
"OnePt aggregate map factory",
"",
124 "Generating factory of (DOF) map for single node aggregates.");
127 validParamList->set< std::string > (
"Interface aggregate map name",
"",
128 "Name of input map for interface aggregates. (default='')");
129 validParamList->set< std::string > (
"Interface aggregate map factory",
"",
130 "Generating factory of (DOF) map for interface aggregates.");
131 validParamList->set<RCP<const FactoryBase> > (
"nodeOnInterface", Teuchos::null,
132 "Array specifying whether or not a node is on the interface (1 or 0).");
136 validParamList->set<std::string > (
"aggregation: mesh layout",
"Global Lexicographic",
137 "Type of mesh ordering");
138 validParamList->set<std::string > (
"aggregation: coupling",
"uncoupled",
139 "aggregation coupling mode: coupled or uncoupled");
140 validParamList->set<
int> (
"aggregation: number of spatial dimensions", 3,
141 "The number of spatial dimensions in the problem");
142 validParamList->set<
int> (
"aggregation: coarsening order", 0,
143 "The interpolation order used to construct grid transfer operators based off these aggregates.");
144 validParamList->set<std::string> (
"aggregation: coarsening rate",
"{3}",
145 "Coarsening rate per spatial dimensions");
146 validParamList->set<RCP<const FactoryBase> >(
"aggregation: mesh data", Teuchos::null,
147 "Mesh ordering associated data");
149 validParamList->set<RCP<const FactoryBase> >(
"gNodesPerDim", Teuchos::null,
150 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
151 validParamList->set<RCP<const FactoryBase> >(
"lNodesPerDim", Teuchos::null,
152 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
156 validParamList->set<RCP<const FactoryBase> > (
"aggregationRegionType", Teuchos::null,
157 "Type of aggregation to use on the region (\"structured\" or \"uncoupled\")");
159 return validParamList;
162 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
165 Input(currentLevel,
"Graph");
167 ParameterList pL = GetParameterList();
176 "Aggregation region type was not provided by the user!");
179 Input(currentLevel,
"aggregationRegionType");
184 Input(currentLevel,
"DofsPerNode");
187 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true){
194 "nodeOnInterface was not provided by the user on level0!");
197 Input(currentLevel,
"nodeOnInterface");
202 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
203 if (mapOnePtName.length() > 0) {
204 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
205 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
208 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
209 currentLevel.
DeclareInput(mapOnePtName, mapOnePtFact.get());
216 std::string coupling = pL.get<std::string>(
"aggregation: coupling");
217 const bool coupled = (coupling ==
"coupled" ? true :
false);
226 "gNodesPerDim was not provided by the user on level0!");
229 Input(currentLevel,
"gNodesPerDim");
240 "lNodesPerDim was not provided by the user on level0!");
243 Input(currentLevel,
"lNodesPerDim");
249 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
254 RCP<Teuchos::FancyOStream> out;
255 if(
const char* dbg = std::getenv(
"MUELU_HYBRIDAGGREGATION_DEBUG")) {
256 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
258 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
260 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
262 *out <<
"Entering hybrid aggregation" << std::endl;
264 ParameterList pL = GetParameterList();
265 bDefinitionPhase_ =
false;
267 if (pL.get<
int>(
"aggregation: max agg size") == -1)
268 pL.set(
"aggregation: max agg size", INT_MAX);
271 RCP<const FactoryBase> graphFact = GetFactory(
"Graph");
274 RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel,
"Graph");
275 RCP<const Map> fineMap = graph->GetDomainMap();
276 const int myRank = fineMap->getComm()->getRank();
277 const int numRanks = fineMap->getComm()->getSize();
278 const GO minGlobalIndex = fineMap->getMinGlobalIndex();
281 RCP<Aggregates> aggregates = rcp(
new Aggregates(*graph));
282 aggregates->setObjectLabel(
"HB");
285 const LO numRows = graph->GetNodeNumVertices();
286 std::vector<unsigned> aggStat(numRows,
READY);
289 std::string regionType;
290 if(currentLevel.GetLevelID() == 0) {
292 regionType = currentLevel.Get< std::string >(
"aggregationRegionType",
NoFactory::get());
295 regionType = Get< std::string >(currentLevel,
"aggregationRegionType");
297 *out<<
"p="<< myRank <<
" | "<<regionType<<
" | regionType determined" << std::endl;
300 if (regionType ==
"structured") {
306 const int numDimensions = pL.get<
int>(
"aggregation: number of spatial dimensions");
307 const int interpolationOrder = pL.get<
int>(
"aggregation: coarsening order");
308 std::string meshLayout = pL.get<std::string>(
"aggregation: mesh layout");
309 std::string coupling = pL.get<std::string>(
"aggregation: coupling");
310 const bool coupled =
false;
311 Array<GO> gFineNodesPerDir(3);
312 Array<LO> lFineNodesPerDir(3);
313 if(currentLevel.GetLevelID() == 0) {
316 gFineNodesPerDir = currentLevel.Get<Array<GO> >(
"gNodesPerDim",
NoFactory::get());
318 lFineNodesPerDir = currentLevel.Get<Array<LO> >(
"lNodesPerDim",
NoFactory::get());
322 gFineNodesPerDir = Get<Array<GO> >(currentLevel,
"gNodesPerDim");
324 lFineNodesPerDir = Get<Array<LO> >(currentLevel,
"lNodesPerDim");
328 for(
int dim = numDimensions; dim < 3; ++dim) {
329 lFineNodesPerDir[dim] = 1;
333 std::string coarseningRate = pL.get<std::string>(
"aggregation: coarsening rate");
334 Teuchos::Array<LO> coarseRate;
336 coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
337 }
catch(
const Teuchos::InvalidArrayStringRepresentation e) {
338 GetOStream(
Errors,-1) <<
" *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** " 342 TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
344 "\"aggregation: coarsening rate\" must have at least as many" 345 " components as the number of spatial dimensions in the problem.");
348 RCP<MueLu::IndexManager<LO,GO,NO> > geoData;
361 "Coupled aggregation is not yet implemented in hybrid aggregation");
364 *out <<
"The index manager has now been built" << std::endl;
365 TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
366 !=
static_cast<size_t>(geoData->getNumLocalFineNodes()),
368 "The local number of elements in the graph's map is not equal to " 369 "the number of nodes given by: lNodesPerDim!");
371 TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getGlobalNumElements()
372 !=
static_cast<size_t>(geoData->getNumGlobalFineNodes()),
374 "The global number of elements in the graph's map is not equal to " 375 "the number of nodes given by: gNodesPerDim!");
378 *out <<
"Compute coarse mesh data" << std::endl;
379 std::vector<std::vector<GO> > coarseMeshData = geoData->getCoarseMeshData();
381 aggregates->SetIndexManager(geoData);
382 aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
384 Set(currentLevel,
"gCoarseNodesPerDim", geoData->getGlobalCoarseNodesPerDir());
385 Set(currentLevel,
"lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
389 if (regionType ==
"uncoupled"){
393 if (pL.get<
bool>(
"aggregation: allow user-specified singletons") ==
true) algos_.push_back(rcp(
new OnePtAggregationAlgorithm (graphFact)));
401 std::string mapInterfaceName = pL.get<std::string>(
"Interface aggregate map name");
402 RCP<Map> InterfaceMap = Teuchos::null;
404 if (pL.get<
bool>(
"aggregation: use interface aggregation") ==
true){
405 Teuchos::Array<LO> nodeOnInterface = Get<Array<LO>>(currentLevel,
"nodeOnInterface");
406 for (LO i = 0; i < numRows; i++) {
407 if (nodeOnInterface[i])
413 ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
414 if (dirichletBoundaryMap != Teuchos::null)
415 for (LO i = 0; i < numRows; i++)
416 if (dirichletBoundaryMap[i] ==
true)
420 std::string mapOnePtName = pL.get<std::string>(
"OnePt aggregate map name");
421 RCP<Map> OnePtMap = Teuchos::null;
422 if (mapOnePtName.length()) {
423 std::string mapOnePtFactName = pL.get<std::string>(
"OnePt aggregate map factory");
424 if (mapOnePtFactName ==
"" || mapOnePtFactName ==
"NoFactory") {
425 OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName,
NoFactory::get());
427 RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
428 OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
431 LO nDofsPerNode = Get<LO>(currentLevel,
"DofsPerNode");
432 GO indexBase = graph->GetDomainMap()->getIndexBase();
433 if (OnePtMap != Teuchos::null) {
434 for (LO i = 0; i < numRows; i++) {
436 GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
437 for (LO kr = 0; kr < nDofsPerNode; kr++)
438 if (OnePtMap->isNodeGlobalElement(grid + kr))
444 Array<LO> lCoarseNodesPerDir(3,-1);
445 Set(currentLevel,
"lCoarseNodesPerDim", lCoarseNodesPerDir);
448 aggregates->AggregatesCrossProcessors(
false);
450 LO numNonAggregatedNodes = numRows;
451 for (
size_t a = 0; a < algos_.size(); a++) {
452 std::string phase = algos_[a]->description();
454 *out <<
"p=" << myRank <<
" | "<<regionType<<
" | Executing phase " << a << std::endl;
456 int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
457 algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
458 algos_[a]->SetProcRankVerbose(oldRank);
459 *out <<
"p=" << myRank <<
" | "<<regionType<<
" | Done Executing phase " << a << std::endl;
462 aggregates->ComputeAggregateSizes(
true);
464 Set(currentLevel,
"Aggregates", aggregates);
466 Set(currentLevel,
"aggregationRegionTypeCoarse", regionType);
468 GetOStream(
Statistics1) << aggregates->description() << std::endl;
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.
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()