MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_StructuredAggregationFactory_kokkos_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_KOKKOS_DEF_HPP
47#define MUELU_STRUCTUREDAGGREGATIONFACTORY_KOKKOS_DEF_HPP
48
49// Xpetra includes
50#include <Xpetra_Map.hpp>
51#include <Xpetra_CrsGraph.hpp>
52
53// MueLu generic includes
54#include "MueLu_Level.hpp"
55#include "MueLu_MasterList.hpp"
56#include "MueLu_Monitor.hpp"
57#include "MueLu_Utilities.hpp"
58
59// MueLu specific includes (kokkos version)
60#include "MueLu_LWGraph_kokkos.hpp"
61#include "MueLu_Aggregates_kokkos.hpp"
62#include "MueLu_IndexManager_kokkos.hpp"
63#include "MueLu_AggregationStructuredAlgorithm_kokkos.hpp"
64
66
67namespace MueLu {
68
69 template <class LocalOrdinal, class GlobalOrdinal, class Node>
72
73 template <class LocalOrdinal, class GlobalOrdinal, class Node>
76 RCP<ParameterList> validParamList = rcp(new ParameterList());
77
78#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
79 SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
80 SET_VALID_ENTRY("aggregation: allow user-specified singletons");
81 SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
82 SET_VALID_ENTRY("aggregation: phase3 avoid singletons");
83#undef SET_VALID_ENTRY
84
85 // general variables needed in StructuredAggregationFactory
86 validParamList->set<std::string> ("aggregation: output type", "Aggregates",
87 "Type of object holding the aggregation data: Aggregtes or CrsGraph");
88 validParamList->set<std::string> ("aggregation: coarsening rate", "{3}",
89 "Coarsening rate per spatial dimensions");
90 validParamList->set<int> ("aggregation: coarsening order", 0,
91 "The interpolation order used to construct grid transfer operators based off these aggregates.");
92 validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null,
93 "Graph of the matrix after amalgamation but without dropping.");
94 validParamList->set<RCP<const FactoryBase> >("DofsPerNode", Teuchos::null,
95 "Number of degrees of freedom per mesh node, provided by the coalsce drop factory.");
96 validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
97 "Number of spatial dimension provided by CoordinatesTransferFactory.");
98 validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
99 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
100
101 return validParamList;
102 } // GetValidParameterList()
103
104 template <class LocalOrdinal, class GlobalOrdinal, class Node>
106 DeclareInput(Level& currentLevel) const {
107 Input(currentLevel, "Graph");
108 Input(currentLevel, "DofsPerNode");
109
110 // Request the local number of nodes per dimensions
111 if(currentLevel.GetLevelID() == 0) {
112 if(currentLevel.IsAvailable("numDimensions", NoFactory::get())) {
113 currentLevel.DeclareInput("numDimensions", NoFactory::get(), this);
114 } else {
115 TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("numDimensions", NoFactory::get()),
117 "numDimensions was not provided by the user on level0!");
118 }
119 if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
120 currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
121 } else {
122 TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
124 "lNodesPerDim was not provided by the user on level0!");
125 }
126 } else {
127 Input(currentLevel, "lNodesPerDim");
128 Input(currentLevel, "numDimensions");
129 }
130 } // DeclareInput()
131
132 template <class LocalOrdinal, class GlobalOrdinal, class Node>
134 Build(Level &currentLevel) const {
135 FactoryMonitor m(*this, "Build", currentLevel);
136
137 RCP<Teuchos::FancyOStream> out;
138 if(const char* dbg = std::getenv("MUELU_STRUCTUREDAGGREGATION_DEBUG")) {
139 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
140 out->setShowAllFrontMatter(false).setShowProcRank(true);
141 } else {
142 out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
143 }
144
145 using device_type = typename LWGraph_kokkos::local_graph_type::device_type;
146 using execution_space = typename LWGraph_kokkos::local_graph_type::device_type::execution_space;
147 using memory_space = typename LWGraph_kokkos::local_graph_type::device_type::memory_space;
148
149 *out << "Entering structured aggregation" << std::endl;
150
151 ParameterList pL = GetParameterList();
152 bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
153
154 // General problem informations are gathered from data stored in the problem matix.
155 RCP<const LWGraph_kokkos> graph = Get<RCP<LWGraph_kokkos> >(currentLevel, "Graph");
156 RCP<const Map> fineMap = graph->GetDomainMap();
157 const int myRank = fineMap->getComm()->getRank();
158 const LO dofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
159
160 // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
161 // obtain a nodeMap.
162 const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
163 std::string outputType = pL.get<std::string>("aggregation: output type");
164 const bool outputAggregates = (outputType == "Aggregates" ? true : false);
165 Array<LO> lFineNodesPerDir(3);
166 int numDimensions;
167 if(currentLevel.GetLevelID() == 0) {
168 // On level 0, data is provided by applications and has no associated factory.
169 lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
170 numDimensions = currentLevel.Get<int>("numDimensions", NoFactory::get());
171 } else {
172 // On level > 0, data is provided directly by generating factories.
173 lFineNodesPerDir = Get<Array<LO> >(currentLevel, "lNodesPerDim");
174 numDimensions = Get<int>(currentLevel, "numDimensions");
175 }
176
177
178 // First make sure that input parameters are set logically based on dimension
179 for(int dim = 0; dim < 3; ++dim) {
180 if(dim >= numDimensions) {
181 lFineNodesPerDir[dim] = 1;
182 }
183 }
184
185 // Get the coarsening rate
186 std::string coarseningRate = pL.get<std::string>("aggregation: coarsening rate");
187 Teuchos::Array<LO> coarseRate;
188 try {
189 coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
190 } catch(const Teuchos::InvalidArrayStringRepresentation& e) {
191 GetOStream(Errors,-1) << " *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
192 << std::endl;
193 throw e;
194 }
195 TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
197 "\"aggregation: coarsening rate\" must have at least as many"
198 " components as the number of spatial dimensions in the problem.");
199
200 // Now that we have extracted info from the level, create the IndexManager
201 RCP<IndexManager_kokkos> geoData = rcp(new IndexManager_kokkos(numDimensions,
202 interpolationOrder, myRank,
203 lFineNodesPerDir,
204 coarseRate));
205
206 *out << "The index manager has now been built" << std::endl;
207 TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getLocalNumElements()
208 != static_cast<size_t>(geoData->getNumLocalFineNodes()),
210 "The local number of elements in the graph's map is not equal to "
211 "the number of nodes given by: lNodesPerDim!");
212
213 // Now we are ready for the big loop over the fine node that will assign each
214 // node on the fine grid to an aggregate and a processor.
215 RCP<AggregationStructuredAlgorithm_kokkos> myStructuredAlgorithm
217
218 if(interpolationOrder == 0 && outputAggregates){
219 RCP<Aggregates_kokkos> aggregates = rcp(new Aggregates_kokkos(graph->GetDomainMap()));
220 aggregates->setObjectLabel("ST");
221 aggregates->SetIndexManager(geoData);
222 aggregates->AggregatesCrossProcessors(false);
223 aggregates->SetNumAggregates(geoData->getNumCoarseNodes());
224
225 LO numNonAggregatedNodes = geoData->getNumLocalFineNodes();
226 Kokkos::View<unsigned*, device_type> aggStat("aggStat", numNonAggregatedNodes);
227 Kokkos::parallel_for("StructuredAggregation: initialize aggStat",
228 Kokkos::RangePolicy<execution_space>(0, numNonAggregatedNodes),
229 KOKKOS_LAMBDA(const LO nodeIdx) {aggStat(nodeIdx) = READY;});
230
231 myStructuredAlgorithm->BuildAggregates(pL, *graph, *aggregates, aggStat,
232 numNonAggregatedNodes);
233
234 *out << "numNonAggregatedNodes: " << numNonAggregatedNodes << std::endl;
235
236 TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError,
237 "MueLu::StructuredAggregationFactory::Build: Leftover nodes found! Error!");
238 aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
239 GetOStream(Statistics1) << aggregates->description() << std::endl;
240 Set(currentLevel, "Aggregates", aggregates);
241
242 } else {
243 // Create Coarse Data
244 RCP<CrsGraph> myGraph;
245 myStructuredAlgorithm->BuildGraph(*graph, geoData, dofsPerNode, myGraph);
246 Set(currentLevel, "prolongatorGraph", myGraph);
247 }
248
249 Set(currentLevel, "lCoarseNodesPerDim", geoData->getCoarseNodesPerDirArray());
250 Set(currentLevel, "indexManager", geoData);
251 Set(currentLevel, "structuredInterpolationOrder", interpolationOrder);
252 Set(currentLevel, "numDimensions", numDimensions);
253
254 } // Build()
255
256} //namespace MueLu
257
258#endif /* MUELU_STRUCTUREDAGGREGATIONFACTORY_KOKKOS_DEF_HPP */
#define SET_VALID_ENTRY(name)
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.
Container class for mesh layout and indices calculation.
Class that holds all level-specific information.
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.
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()
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.