MueLu Version of the Day
Loading...
Searching...
No Matches
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
65namespace 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->getLocalNumElements()
287 << ", structured aggregation num nodes: " << geoData->getNumLocalFineNodes() << std::endl;
288 TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getLocalNumElements()
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_ */
#define SET_VALID_ENTRY(name)
Container class for aggregation information.
Algorithm for coarsening a graph with structured aggregation.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
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's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
void Build(Level &currentLevel) const
Build aggregates.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Namespace for MueLu classes and methods.
@ Errors
Errors.
@ Statistics1
Print more statistics.