MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_StructuredAggregationFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_STRUCTUREDAGGREGATIONFACTORY_DEF_HPP_
47 #define MUELU_STRUCTUREDAGGREGATIONFACTORY_DEF_HPP_
48 
49 #include <Xpetra_Map.hpp>
50 #include <Xpetra_CrsGraph.hpp>
51 
52 #include "MueLu_AggregationStructuredAlgorithm.hpp"
53 #include "MueLu_Level.hpp"
54 #include "MueLu_GraphBase.hpp"
55 #include "MueLu_Aggregates.hpp"
56 #include "MueLu_MasterList.hpp"
57 #include "MueLu_Monitor.hpp"
58 #include "MueLu_Utilities.hpp"
59 #include "MueLu_UncoupledIndexManager.hpp"
60 #include "MueLu_LocalLexicographicIndexManager.hpp"
61 #include "MueLu_GlobalLexicographicIndexManager.hpp"
62 
64 
65 namespace MueLu {
66 
67  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69  StructuredAggregationFactory() : bDefinitionPhase_(true)
70  { }
71 
72  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75  RCP<ParameterList> validParamList = rcp(new ParameterList());
76 
77 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
78  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
79  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
80  SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
81  SET_VALID_ENTRY("aggregation: phase3 avoid singletons");
82 
83  // general variables needed in StructuredAggregationFactory
84  SET_VALID_ENTRY("aggregation: mesh layout");
85  SET_VALID_ENTRY("aggregation: mode");
86  SET_VALID_ENTRY("aggregation: output type");
87  SET_VALID_ENTRY("aggregation: coarsening rate");
88  SET_VALID_ENTRY("aggregation: coarsening order");
89 #undef SET_VALID_ENTRY
90  validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null,
91  "Graph of the matrix after amalgamation but without dropping.");
92  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
93  "Number of spatial dimension provided by CoordinatesTransferFactory.");
94  validParamList->set<RCP<const FactoryBase> >("gNodesPerDim", Teuchos::null,
95  "Global number of nodes per spatial dimension provided by CoordinatesTransferFactory.");
96  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
97  "Local number of nodes per spatial dimension provided by CoordinatesTransferFactory.");
98  validParamList->set<RCP<const FactoryBase> >("DofsPerNode", Teuchos::null,
99  "Generating factory for variable \'DofsPerNode\', usually the same as the \'Graph\' factory");
100 
101  return validParamList;
102  } // GetValidParameterList()
103 
104  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
106  DeclareInput(Level& currentLevel) const {
107  Input(currentLevel, "Graph");
108  Input(currentLevel, "DofsPerNode");
109 
110  ParameterList pL = GetParameterList();
111  std::string coupling = pL.get<std::string>("aggregation: mode");
112  const bool coupled = (coupling == "coupled" ? true : false);
113  if(coupled) {
114  // Request the global number of nodes per dimensions
115  if(currentLevel.GetLevelID() == 0) {
116  if(currentLevel.IsAvailable("gNodesPerDim", NoFactory::get())) {
117  currentLevel.DeclareInput("gNodesPerDim", NoFactory::get(), this);
118  } else {
119  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("gNodesPerDim", NoFactory::get()),
121  "gNodesPerDim was not provided by the user on level0!");
122  }
123  } else {
124  Input(currentLevel, "gNodesPerDim");
125  }
126  }
127 
128  // Request the local number of nodes per dimensions
129  if(currentLevel.GetLevelID() == 0) {
130  if(currentLevel.IsAvailable("numDimensions", NoFactory::get())) {
131  currentLevel.DeclareInput("numDimensions", NoFactory::get(), this);
132  } else {
133  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("numDimensions", NoFactory::get()),
135  "numDimensions was not provided by the user on level0!");
136  }
137  if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
138  currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
139  } else {
140  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
142  "lNodesPerDim was not provided by the user on level0!");
143  }
144  } else {
145  Input(currentLevel, "numDimensions");
146  Input(currentLevel, "lNodesPerDim");
147  }
148  } // DeclareInput()
149 
150  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
152  Build(Level &currentLevel) const {
153  FactoryMonitor m(*this, "Build", currentLevel);
154 
156  if(const char* dbg = std::getenv("MUELU_STRUCTUREDAGGREGATION_DEBUG")) {
157  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
158  out->setShowAllFrontMatter(false).setShowProcRank(true);
159  } else {
160  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
161  }
162 
163  *out << "Entering structured aggregation" << std::endl;
164 
165  ParameterList pL = GetParameterList();
166  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
167 
168  // General problem informations are gathered from data stored in the problem matix.
169  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
170  RCP<const Map> fineMap = graph->GetDomainMap();
171  const int myRank = fineMap->getComm()->getRank();
172  const int numRanks = fineMap->getComm()->getSize();
173  const GO minGlobalIndex = fineMap->getMinGlobalIndex();
174  const LO dofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
175 
176  // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
177  // obtain a nodeMap.
178  const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
179  std::string meshLayout = pL.get<std::string>("aggregation: mesh layout");
180  std::string coupling = pL.get<std::string>("aggregation: mode");
181  const bool coupled = (coupling == "coupled" ? true : false);
182  std::string outputType = pL.get<std::string>("aggregation: output type");
183  const bool outputAggregates = (outputType == "Aggregates" ? true : false);
184  int numDimensions;
185  Array<GO> gFineNodesPerDir(3);
186  Array<LO> lFineNodesPerDir(3);
187  if(currentLevel.GetLevelID() == 0) {
188  // On level 0, data is provided by applications and has no associated factory.
189  numDimensions = currentLevel.Get<int>("numDimensions", NoFactory::get());
190  lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
191  if(coupled) {
192  gFineNodesPerDir = currentLevel.Get<Array<GO> >("gNodesPerDim", NoFactory::get());
193  }
194  } else {
195  // On level > 0, data is provided directly by generating factories.
196  numDimensions = Get<int>(currentLevel, "numDimensions");
197  lFineNodesPerDir = Get<Array<LO> >(currentLevel, "lNodesPerDim");
198  if(coupled) {
199  gFineNodesPerDir = Get<Array<GO> >(currentLevel, "gNodesPerDim");
200  }
201  }
202 
203 
204  // First make sure that input parameters are set logically based on dimension
205  for(int dim = 0; dim < 3; ++dim) {
206  if(dim >= numDimensions) {
207  gFineNodesPerDir[dim] = 1;
208  lFineNodesPerDir[dim] = 1;
209  }
210  }
211 
212  // Get the coarsening rate
213  std::string coarseningRate = pL.get<std::string>("aggregation: coarsening rate");
214  Teuchos::Array<LO> coarseRate;
215  try {
216  coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
218  GetOStream(Errors,-1) << " *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
219  << std::endl;
220  throw e;
221  }
222  TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
224  "\"aggregation: coarsening rate\" must have at least as many"
225  " components as the number of spatial dimensions in the problem.");
226 
227  // Now that we have extracted info from the level, create the IndexManager
228  RCP<IndexManager > geoData;
229  if(!coupled) {
230  geoData = rcp(new MueLu::UncoupledIndexManager<LO,GO,NO>(fineMap->getComm(),
231  coupled,
232  numDimensions,
233  interpolationOrder,
234  myRank,
235  numRanks,
236  gFineNodesPerDir,
237  lFineNodesPerDir,
238  coarseRate));
239  } else if(meshLayout == "Local Lexicographic") {
240  Array<GO> meshData;
241  if(currentLevel.GetLevelID() == 0) {
242  // On level 0, data is provided by applications and has no associated factory.
243  meshData = currentLevel.Get<Array<GO> >("aggregation: mesh data", NoFactory::get());
245  "The meshData array is empty, somehow the input for structured"
246  " aggregation are not captured correctly.");
247  } else {
248  // On level > 0, data is provided directly by generating factories.
249  meshData = Get<Array<GO> >(currentLevel, "aggregation: mesh data");
250  }
251  // Note, LBV Feb 5th 2018:
252  // I think that it might make sense to pass ghostInterface rather than interpolationOrder.
253  // For that I need to make sure that ghostInterface can be computed with minimal mesh
254  // knowledge outside of the IndexManager...
255  geoData = rcp(new MueLu::LocalLexicographicIndexManager<LO,GO,NO>(fineMap->getComm(),
256  coupled,
257  numDimensions,
258  interpolationOrder,
259  myRank,
260  numRanks,
261  gFineNodesPerDir,
262  lFineNodesPerDir,
263  coarseRate,
264  meshData));
265  } else if(meshLayout == "Global Lexicographic") {
266  // Note, LBV Feb 5th 2018:
267  // I think that it might make sense to pass ghostInterface rather than interpolationOrder.
268  // For that I need to make sure that ghostInterface can be computed with minimal mesh
269  // knowledge outside of the IndexManager...
270  geoData = rcp(new MueLu::GlobalLexicographicIndexManager<LO,GO,NO>(fineMap->getComm(),
271  coupled,
272  numDimensions,
273  interpolationOrder,
274  gFineNodesPerDir,
275  lFineNodesPerDir,
276  coarseRate,
277  minGlobalIndex));
278  }
279 
280 
281  *out << "The index manager has now been built" << std::endl;
282  *out << "graph num nodes: " << fineMap->getNodeNumElements()
283  << ", structured aggregation num nodes: " << geoData->getNumLocalFineNodes() << std::endl;
284  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
285  != static_cast<size_t>(geoData->getNumLocalFineNodes()),
287  "The local number of elements in the graph's map is not equal to "
288  "the number of nodes given by: lNodesPerDim!");
289  if(coupled) {
290  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getGlobalNumElements()
291  != static_cast<size_t>(geoData->getNumGlobalFineNodes()),
293  "The global number of elements in the graph's map is not equal to "
294  "the number of nodes given by: gNodesPerDim!");
295  }
296 
297  *out << "Compute coarse mesh data" << std::endl;
298  std::vector<std::vector<GO> > coarseMeshData = geoData->getCoarseMeshData();
299 
300  // Now we are ready for the big loop over the fine node that will assign each
301  // node on the fine grid to an aggregate and a processor.
302  RCP<const FactoryBase> graphFact = GetFactory("Graph");
303  RCP<const Map> coarseCoordinatesFineMap, coarseCoordinatesMap;
305  myStructuredAlgorithm = rcp(new AggregationStructuredAlgorithm(graphFact));
306 
307  if(interpolationOrder == 0 && outputAggregates){
308  // Create aggregates for prolongation
309  RCP<Aggregates> aggregates = rcp(new Aggregates(graph->GetDomainMap()));
310  aggregates->setObjectLabel("ST");
311  aggregates->SetIndexManager(geoData);
312  aggregates->AggregatesCrossProcessors(coupled);
313  aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
314  std::vector<unsigned> aggStat(geoData->getNumLocalFineNodes(), READY);
315  LO numNonAggregatedNodes = geoData->getNumLocalFineNodes();
316 
317  myStructuredAlgorithm->BuildAggregates(pL, *graph, *aggregates, aggStat,
318  numNonAggregatedNodes);
319 
321  "MueLu::StructuredAggregationFactory::Build: Leftover nodes found! Error!");
322  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
323  GetOStream(Statistics1) << aggregates->description() << std::endl;
324  Set(currentLevel, "Aggregates", aggregates);
325 
326  } else {
327  // Create the graph of the prolongator
328  RCP<CrsGraph> myGraph;
329  myStructuredAlgorithm->BuildGraph(*graph, geoData, dofsPerNode, myGraph,
330  coarseCoordinatesFineMap, coarseCoordinatesMap);
331  Set(currentLevel, "prolongatorGraph", myGraph);
332  }
333 
334  if(coupled) {
335  Set(currentLevel, "gCoarseNodesPerDim", geoData->getGlobalCoarseNodesPerDir());
336  }
337  Set(currentLevel, "lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
338  Set(currentLevel, "coarseCoordinatesFineMap", coarseCoordinatesFineMap);
339  Set(currentLevel, "coarseCoordinatesMap", coarseCoordinatesMap);
340  Set(currentLevel, "interpolationOrder", interpolationOrder);
341  Set(currentLevel, "numDimensions", numDimensions);
342 
343  } // Build()
344 } //namespace MueLu
345 
346 
347 #endif /* MUELU_STRUCTUREDAGGREGATIONFACTORY_DEF_HPP_ */
Container class for aggregation information.
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
bool empty() const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Algorithm for coarsening a graph with structured aggregation.
Print more statistics.
virtual const RCP< const Map > GetDomainMap() const =0
void SetIndexManager(RCP< IndexManager > &geoData)
Get the index manager used by structured aggregation algorithms.
Teuchos::ArrayRCP< LO > ComputeAggregateSizes(bool forceRecompute=false) const
Compute sizes of aggregates.
static const NoFactory * get()
void DeclareInput(Level &currentLevel) const
Input.
Array< GO > getGlobalCoarseNodesPerDir() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
virtual void setObjectLabel(const std::string &objectLabel)
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
virtual std::vector< std::vector< GO > > getCoarseMeshData() const =0
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Exception throws to report errors in the internal logical of the program.
std::string description() const
Return a simple one-line description of this object.
void AggregatesCrossProcessors(const bool &flag)
Record whether aggregates include DOFs from other processes.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
void Build(Level &currentLevel) const
Build aggregates.
#define SET_VALID_ENTRY(name)
Array< LO > getLocalCoarseNodesPerDir() const
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.