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 #undef SET_VALID_ENTRY
82 
83  // general variables needed in StructuredAggregationFactory
84  validParamList->set<std::string> ("aggregation: mesh layout","Global Lexicographic",
85  "Type of mesh ordering");
86  validParamList->set<std::string> ("aggregation: coupling","coupled",
87  "aggregation coupling mode: coupled or uncoupled");
88  validParamList->set<std::string> ("aggregation: output type", "Aggregates",
89  "Type of object holding the aggregation data: Aggregtes or CrsGraph");
90  validParamList->set<std::string> ("aggregation: coarsening rate", "{3}",
91  "Coarsening rate per spatial dimensions");
92  validParamList->set<int> ("aggregation: number of spatial dimensions", 3,
93  "The number of spatial dimensions in the problem");
94  validParamList->set<int> ("aggregation: coarsening order", 0,
95  "The interpolation order used to construct grid transfer operators based off these aggregates.");
96 
97  validParamList->set<RCP<const FactoryBase> >("aggregation: mesh data", Teuchos::null,
98  "Mesh ordering associated data");
99 
100  validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null,
101  "Graph of the matrix after amalgamation but without dropping.");
102  validParamList->set<RCP<const FactoryBase> >("gNodesPerDim", Teuchos::null,
103  "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
104  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
105  "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
106 
107  return validParamList;
108  } // GetValidParameterList
109 
110  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112  DeclareInput(Level& currentLevel) const {
113  Input(currentLevel, "Graph");
114 
115  ParameterList pL = GetParameterList();
116  std::string coupling = pL.get<std::string>("aggregation: coupling");
117  const bool coupled = (coupling == "coupled" ? true : false);
118  if(coupled) {
119  // Request the global number of nodes per dimensions
120  if(currentLevel.GetLevelID() == 0) {
121  if(currentLevel.IsAvailable("gNodesPerDim", NoFactory::get())) {
122  currentLevel.DeclareInput("gNodesPerDim", NoFactory::get(), this);
123  } else {
124  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("gNodesPerDim", NoFactory::get()),
126  "gNodesPerDim was not provided by the user on level0!");
127  }
128  } else {
129  Input(currentLevel, "gNodesPerDim");
130  }
131  }
132 
133  // Request the local number of nodes per dimensions
134  if(currentLevel.GetLevelID() == 0) {
135  if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
136  currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
137  } else {
138  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
140  "lNodesPerDim was not provided by the user on level0!");
141  }
142  } else {
143  Input(currentLevel, "lNodesPerDim");
144  }
145  } // DeclareInput
146 
147  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
149  Build(Level &currentLevel) const {
150  FactoryMonitor m(*this, "Build", currentLevel);
151 
153  if(const char* dbg = std::getenv("MUELU_STRUCTUREDAGGREGATION_DEBUG")) {
154  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
155  out->setShowAllFrontMatter(false).setShowProcRank(true);
156  } else {
157  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
158  }
159 
160  *out << "Entering structured aggregation" << std::endl;
161 
162  ParameterList pL = GetParameterList();
163  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
164 
165  // General problem informations are gathered from data stored in the problem matix.
166  RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
167  RCP<const Map> fineMap = graph->GetDomainMap();
168  const int myRank = fineMap->getComm()->getRank();
169  const int numRanks = fineMap->getComm()->getSize();
170  const GO minGlobalIndex = fineMap->getMinGlobalIndex();
171 
172  // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
173  // obtain a nodeMap.
174  const int numDimensions = pL.get<int>("aggregation: number of spatial dimensions");
175  const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
176  std::string meshLayout = pL.get<std::string>("aggregation: mesh layout");
177  std::string coupling = pL.get<std::string>("aggregation: coupling");
178  const bool coupled = (coupling == "coupled" ? true : false);
179  std::string outputType = pL.get<std::string>("aggregation: output type");
180  const bool outputAggregates = (outputType == "Aggregates" ? true : false);
181  Array<GO> gFineNodesPerDir(3);
182  Array<LO> lFineNodesPerDir(3);
183  if(currentLevel.GetLevelID() == 0) {
184  // On level 0, data is provided by applications and has no associated factory.
185  if(coupled) {
186  gFineNodesPerDir = currentLevel.Get<Array<GO> >("gNodesPerDim", NoFactory::get());
187  }
188  lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
189  } else {
190  // On level > 0, data is provided directly by generating factories.
191  if(coupled) {
192  gFineNodesPerDir = Get<Array<GO> >(currentLevel, "gNodesPerDim");
193  }
194  lFineNodesPerDir = Get<Array<LO> >(currentLevel, "lNodesPerDim");
195  }
196 
197 
198  // First make sure that input parameters are set logically based on dimension
199  for(int dim = 0; dim < 3; ++dim) {
200  if(dim >= numDimensions) {
201  gFineNodesPerDir[dim] = 1;
202  lFineNodesPerDir[dim] = 1;
203  }
204  }
205 
206  // Get the coarsening rate
207  std::string coarseningRate = pL.get<std::string>("aggregation: coarsening rate");
208  Teuchos::Array<LO> coarseRate;
209  try {
210  coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
212  GetOStream(Errors,-1) << " *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
213  << std::endl;
214  throw e;
215  }
216  TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
218  "\"aggregation: coarsening rate\" must have at least as many"
219  " components as the number of spatial dimensions in the problem.");
220 
221  // Now that we have extracted info from the level, create the IndexManager
222  RCP<IndexManager > geoData;
223  if(!coupled) {
224  geoData = rcp(new MueLu::UncoupledIndexManager<LO,GO,NO>(fineMap->getComm(),
225  coupled,
226  numDimensions,
227  interpolationOrder,
228  myRank,
229  numRanks,
230  gFineNodesPerDir,
231  lFineNodesPerDir,
232  coarseRate));
233  } else if(meshLayout == "Local Lexicographic") {
234  Array<GO> meshData;
235  if(currentLevel.GetLevelID() == 0) {
236  // On level 0, data is provided by applications and has no associated factory.
237  meshData = currentLevel.Get<Array<GO> >("aggregation: mesh data", NoFactory::get());
239  "The meshData array is empty, somehow the input for structured"
240  " aggregation are not captured correctly.");
241  } else {
242  // On level > 0, data is provided directly by generating factories.
243  meshData = Get<Array<GO> >(currentLevel, "aggregation: mesh data");
244  }
245  // Note, LBV Feb 5th 2018:
246  // I think that it might make sense to pass ghostInterface rather than interpolationOrder.
247  // For that I need to make sure that ghostInterface can be computed with minimal mesh
248  // knowledge outside of the IndexManager...
249  geoData = rcp(new MueLu::LocalLexicographicIndexManager<LO,GO,NO>(fineMap->getComm(),
250  coupled,
251  numDimensions,
252  interpolationOrder,
253  myRank,
254  numRanks,
255  gFineNodesPerDir,
256  lFineNodesPerDir,
257  coarseRate,
258  meshData));
259  } else if(meshLayout == "Global Lexicographic") {
260  // Note, LBV Feb 5th 2018:
261  // I think that it might make sense to pass ghostInterface rather than interpolationOrder.
262  // For that I need to make sure that ghostInterface can be computed with minimal mesh
263  // knowledge outside of the IndexManager...
264  geoData = rcp(new MueLu::GlobalLexicographicIndexManager<LO,GO,NO>(fineMap->getComm(),
265  coupled,
266  numDimensions,
267  interpolationOrder,
268  gFineNodesPerDir,
269  lFineNodesPerDir,
270  coarseRate,
271  minGlobalIndex));
272  }
273 
274 
275  *out << "The index manager has now been built" << std::endl;
276  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
277  != static_cast<size_t>(geoData->getNumLocalFineNodes()),
279  "The local number of elements in the graph's map is not equal to "
280  "the number of nodes given by: lNodesPerDim!");
281  if(coupled) {
282  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getGlobalNumElements()
283  != static_cast<size_t>(geoData->getNumGlobalFineNodes()),
285  "The global number of elements in the graph's map is not equal to "
286  "the number of nodes given by: gNodesPerDim!");
287  }
288 
289  *out << "Compute coarse mesh data" << std::endl;
290  std::vector<std::vector<GO> > coarseMeshData = geoData->getCoarseMeshData();
291 
292  // Now we are ready for the big loop over the fine node that will assign each
293  // node on the fine grid to an aggregate and a processor.
294  RCP<const FactoryBase> graphFact = GetFactory("Graph");
295  RCP<const Map> coarseCoordinatesFineMap, coarseCoordinatesMap;
297  myStructuredAlgorithm = rcp(new AggregationStructuredAlgorithm(graphFact));
298 
299  if(interpolationOrder == 0 && outputAggregates){
300  RCP<Aggregates> aggregates = rcp(new Aggregates(graph->GetDomainMap()));
301  aggregates->setObjectLabel("ST");
302  aggregates->SetIndexManager(geoData);
303  aggregates->AggregatesCrossProcessors(coupled);
304  aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());
305  std::vector<unsigned> aggStat(geoData->getNumLocalFineNodes(), READY);
306  LO numNonAggregatedNodes = geoData->getNumLocalFineNodes();
307 
308  myStructuredAlgorithm->BuildAggregates(pL, *graph, *aggregates, aggStat,
309  numNonAggregatedNodes);
310 
312  "MueLu::StructuredAggregationFactory::Build: Leftover nodes found! Error!");
313  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
314  GetOStream(Statistics1) << aggregates->description() << std::endl;
315  Set(currentLevel, "Aggregates", aggregates);
316 
317  } else {
318  // Create Coarse Data
319  RCP<CrsGraph> myGraph;
320  myStructuredAlgorithm->BuildGraph(*graph, geoData, myGraph, coarseCoordinatesFineMap,
321  coarseCoordinatesMap);
322  Set(currentLevel, "prolongatorGraph", myGraph);
323  }
324 
325  if(coupled) {
326  Set(currentLevel, "gCoarseNodesPerDim", geoData->getGlobalCoarseNodesPerDir());
327  }
328  Set(currentLevel, "lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());
329  Set(currentLevel, "coarseCoordinatesFineMap", coarseCoordinatesFineMap);
330  Set(currentLevel, "coarseCoordinatesMap", coarseCoordinatesMap);
331  Set(currentLevel, "interpolationOrder", interpolationOrder);
332  Set(currentLevel, "numDimensions", numDimensions);
333 
334  } // Build
335 } //namespace MueLu
336 
337 
338 #endif /* MUELU_STRUCTUREDAGGREGATIONFACTORY_DEF_HPP_ */
std::string description() const
Return a simple one-line description of this object.
void Build(Level &currentLevel) const
Build aggregates.
Container class for aggregation information.
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
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
Namespace for MueLu classes and methods.
void SetIndexManager(RCP< IndexManager > &geoData)
Get the index manager used by structured aggregation algorithms.
static const NoFactory * get()
Array< LO > getLocalCoarseNodesPerDir() const
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
virtual void setObjectLabel(const std::string &objectLabel)
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
bool empty() const
virtual std::vector< std::vector< GO > > getCoarseMeshData() const =0
Array< GO > getGlobalCoarseNodesPerDir() const
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 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()
#define SET_VALID_ENTRY(name)
Teuchos::ArrayRCP< LO > ComputeAggregateSizes(bool forceRecompute=false) const
Compute sizes of aggregates.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.