MueLu  Version of the Day
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  validParamList->set<const bool>("aggregation: single coarse point", false,
101  "Allows the aggreagtion process to reduce spacial dimensions to a single layer");
102 
103  return validParamList;
104  } // GetValidParameterList()
105 
106  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
108  DeclareInput(Level& currentLevel) const {
109  Input(currentLevel, "Graph");
110  Input(currentLevel, "DofsPerNode");
111 
112  ParameterList pL = GetParameterList();
113  std::string coupling = pL.get<std::string>("aggregation: mode");
114  const bool coupled = (coupling == "coupled" ? true : false);
115  if(coupled) {
116  // Request the global number of nodes per dimensions
117  if(currentLevel.GetLevelID() == 0) {
118  if(currentLevel.IsAvailable("gNodesPerDim", NoFactory::get())) {
119  currentLevel.DeclareInput("gNodesPerDim", NoFactory::get(), this);
120  } else {
121  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("gNodesPerDim", NoFactory::get()),
123  "gNodesPerDim was not provided by the user on level0!");
124  }
125  } else {
126  Input(currentLevel, "gNodesPerDim");
127  }
128  }
129 
130  // Request the local number of nodes per dimensions
131  if(currentLevel.GetLevelID() == 0) {
132  if(currentLevel.IsAvailable("numDimensions", NoFactory::get())) {
133  currentLevel.DeclareInput("numDimensions", NoFactory::get(), this);
134  } else {
135  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("numDimensions", NoFactory::get()),
137  "numDimensions was not provided by the user on level0!");
138  }
139  if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
140  currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
141  } else {
142  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
144  "lNodesPerDim was not provided by the user on level0!");
145  }
146  } else {
147  Input(currentLevel, "numDimensions");
148  Input(currentLevel, "lNodesPerDim");
149  }
150  } // DeclareInput()
151 
152  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
154  Build(Level &currentLevel) const {
155  FactoryMonitor m(*this, "Build", currentLevel);
156 
157  RCP<Teuchos::FancyOStream> out;
158  if(const char* dbg = std::getenv("MUELU_STRUCTUREDAGGREGATION_DEBUG")) {
159  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
160  out->setShowAllFrontMatter(false).setShowProcRank(true);
161  } else {
162  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
163  }
164 
165  *out << "Entering structured aggregation" << std::endl;
166 
167  ParameterList pL = GetParameterList();
168  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
169 
170  // General problem informations are gathered from data stored in the problem matix.
171  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
172  RCP<const Map> fineMap = graph->GetDomainMap();
173  const int myRank = fineMap->getComm()->getRank();
174  const int numRanks = fineMap->getComm()->getSize();
175  const GO minGlobalIndex = fineMap->getMinGlobalIndex();
176  const LO dofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
177 
178  // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
179  // obtain a nodeMap.
180  const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
181  std::string meshLayout = pL.get<std::string>("aggregation: mesh layout");
182  std::string coupling = pL.get<std::string>("aggregation: mode");
183  const bool coupled = (coupling == "coupled" ? true : false);
184  std::string outputType = pL.get<std::string>("aggregation: output type");
185  const bool outputAggregates = (outputType == "Aggregates" ? true : false);
186  const bool singleCoarsePoint = pL.get<bool>("aggregation: single coarse point");
187  int numDimensions;
188  Array<GO> gFineNodesPerDir(3);
189  Array<LO> lFineNodesPerDir(3);
190  if(currentLevel.GetLevelID() == 0) {
191  // On level 0, data is provided by applications and has no associated factory.
192  numDimensions = currentLevel.Get<int>("numDimensions", NoFactory::get());
193  lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
194  if(coupled) {
195  gFineNodesPerDir = currentLevel.Get<Array<GO> >("gNodesPerDim", NoFactory::get());
196  }
197  } else {
198  // On level > 0, data is provided directly by generating factories.
199  numDimensions = Get<int>(currentLevel, "numDimensions");
200  lFineNodesPerDir = Get<Array<LO> >(currentLevel, "lNodesPerDim");
201  if(coupled) {
202  gFineNodesPerDir = Get<Array<GO> >(currentLevel, "gNodesPerDim");
203  }
204  }
205 
206 
207  // First make sure that input parameters are set logically based on dimension
208  for(int dim = 0; dim < 3; ++dim) {
209  if(dim >= numDimensions) {
210  gFineNodesPerDir[dim] = 1;
211  lFineNodesPerDir[dim] = 1;
212  }
213  }
214 
215  // Get the coarsening rate
216  std::string coarseningRate = pL.get<std::string>("aggregation: coarsening rate");
217  Teuchos::Array<LO> coarseRate;
218  try {
219  coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
220  } catch(const Teuchos::InvalidArrayStringRepresentation& e) {
221  GetOStream(Errors,-1) << " *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
222  << std::endl;
223  throw e;
224  }
225  TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
227  "\"aggregation: coarsening rate\" must have at least as many"
228  " components as the number of spatial dimensions in the problem.");
229 
230  // Now that we have extracted info from the level, create the IndexManager
231  RCP<IndexManager > geoData;
232  if(!coupled) {
233  geoData = rcp(new MueLu::UncoupledIndexManager<LO,GO,NO>(fineMap->getComm(),
234  coupled,
235  numDimensions,
236  interpolationOrder,
237  myRank,
238  numRanks,
239  gFineNodesPerDir,
240  lFineNodesPerDir,
241  coarseRate,
242  singleCoarsePoint));
243  } else if(meshLayout == "Local Lexicographic") {
244  Array<GO> meshData;
245  if(currentLevel.GetLevelID() == 0) {
246  // On level 0, data is provided by applications and has no associated factory.
247  meshData = currentLevel.Get<Array<GO> >("aggregation: mesh data", NoFactory::get());
248  TEUCHOS_TEST_FOR_EXCEPTION(meshData.empty() == true, Exceptions::RuntimeError,
249  "The meshData array is empty, somehow the input for structured"
250  " aggregation are not captured correctly.");
251  } else {
252  // On level > 0, data is provided directly by generating factories.
253  meshData = Get<Array<GO> >(currentLevel, "aggregation: mesh data");
254  }
255  // Note, LBV Feb 5th 2018:
256  // I think that it might make sense to pass ghostInterface rather than interpolationOrder.
257  // For that I need to make sure that ghostInterface can be computed with minimal mesh
258  // knowledge outside of the IndexManager...
259  geoData = rcp(new MueLu::LocalLexicographicIndexManager<LO,GO,NO>(fineMap->getComm(),
260  coupled,
261  numDimensions,
262  interpolationOrder,
263  myRank,
264  numRanks,
265  gFineNodesPerDir,
266  lFineNodesPerDir,
267  coarseRate,
268  meshData));
269  } else if(meshLayout == "Global Lexicographic") {
270  // Note, LBV Feb 5th 2018:
271  // I think that it might make sense to pass ghostInterface rather than interpolationOrder.
272  // For that I need to make sure that ghostInterface can be computed with minimal mesh
273  // knowledge outside of the IndexManager...
274  geoData = rcp(new MueLu::GlobalLexicographicIndexManager<LO,GO,NO>(fineMap->getComm(),
275  coupled,
276  numDimensions,
277  interpolationOrder,
278  gFineNodesPerDir,
279  lFineNodesPerDir,
280  coarseRate,
281  minGlobalIndex));
282  }
283 
284 
285  *out << "The index manager has now been built" << std::endl;
286  *out << "graph num nodes: " << fineMap->getNodeNumElements()
287  << ", structured aggregation num nodes: " << geoData->getNumLocalFineNodes() << std::endl;
288  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
289  != static_cast<size_t>(geoData->getNumLocalFineNodes()),
291  "The local number of elements in the graph's map is not equal to "
292  "the number of nodes given by: lNodesPerDim!");
293  if(coupled) {
294  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getGlobalNumElements()
295  != static_cast<size_t>(geoData->getNumGlobalFineNodes()),
297  "The global number of elements in the graph's map is not equal to "
298  "the number of nodes given by: gNodesPerDim!");
299  }
300 
301  *out << "Compute coarse mesh data" << std::endl;
302  std::vector<std::vector<GO> > coarseMeshData = geoData->getCoarseMeshData();
303 
304  // Now we are ready for the big loop over the fine node that will assign each
305  // node on the fine grid to an aggregate and a processor.
306  RCP<const FactoryBase> graphFact = GetFactory("Graph");
307  RCP<const Map> coarseCoordinatesFineMap, coarseCoordinatesMap;
308  RCP<MueLu::AggregationStructuredAlgorithm<LocalOrdinal, GlobalOrdinal, Node> >
309  myStructuredAlgorithm = rcp(new AggregationStructuredAlgorithm(graphFact));
310 
311  if(interpolationOrder == 0 && outputAggregates){
312  // Create aggregates for prolongation
313  *out << "Compute Aggregates" << std::endl;
314  RCP<Aggregates> aggregates = rcp(new Aggregates(graph->GetDomainMap()));
315  aggregates->setObjectLabel("ST");
316  aggregates->SetIndexManager(geoData);
317  aggregates->AggregatesCrossProcessors(coupled);
318  aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
319  std::vector<unsigned> aggStat(geoData->getNumLocalFineNodes(), READY);
320  LO numNonAggregatedNodes = geoData->getNumLocalFineNodes();
321 
322  myStructuredAlgorithm->BuildAggregates(pL, *graph, *aggregates, aggStat,
323  numNonAggregatedNodes);
324 
325  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError,
326  "MueLu::StructuredAggregationFactory::Build: Leftover nodes found! Error!");
327  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
328  GetOStream(Statistics1) << aggregates->description() << std::endl;
329  Set(currentLevel, "Aggregates", aggregates);
330 
331  } else {
332  // Create the graph of the prolongator
333  *out << "Compute CrsGraph" << std::endl;
334  RCP<CrsGraph> myGraph;
335  myStructuredAlgorithm->BuildGraph(*graph, geoData, dofsPerNode, myGraph,
336  coarseCoordinatesFineMap, coarseCoordinatesMap);
337  Set(currentLevel, "prolongatorGraph", myGraph);
338  }
339 
340  if(coupled) {
341  Set(currentLevel, "gCoarseNodesPerDim", geoData->getGlobalCoarseNodesPerDir());
342  }
343  Set(currentLevel, "lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
344  Set(currentLevel, "coarseCoordinatesFineMap", coarseCoordinatesFineMap);
345  Set(currentLevel, "coarseCoordinatesMap", coarseCoordinatesMap);
346  Set(currentLevel, "structuredInterpolationOrder", interpolationOrder);
347  Set(currentLevel, "numDimensions", numDimensions);
348 
349  } // Build()
350 } //namespace MueLu
351 
352 
353 #endif /* MUELU_STRUCTUREDAGGREGATIONFACTORY_DEF_HPP_ */
void Build(Level &currentLevel) const
Build aggregates.
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.
Print more statistics.
Namespace for MueLu classes and methods.
static const NoFactory * get()
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report errors in the internal logical of the program.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
#define SET_VALID_ENTRY(name)