MueLu  Version of the Day
MueLu_CoarseMapFactory_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 /*
47  * MueLu_CoarseMapFactory_def.hpp
48  *
49  * Created on: Oct 12, 2012
50  * Author: wiesner
51  */
52 
53 #ifndef MUELU_COARSEMAPFACTORY_DEF_HPP_
54 #define MUELU_COARSEMAPFACTORY_DEF_HPP_
55 
56 #include <Teuchos_Array.hpp>
57 
58 #include <Xpetra_MultiVector.hpp>
59 #include <Xpetra_StridedMapFactory.hpp>
60 
62 #include "MueLu_Level.hpp"
63 #include "MueLu_Aggregates.hpp"
64 #include "MueLu_Monitor.hpp"
65 
66 namespace MueLu {
67 
68  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  {
71  RCP<ParameterList> validParamList = rcp(new ParameterList());
72 
73  validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Generating factory for aggregates.");
74  validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory for null space.");
75 
76  validParamList->set< std::string >("Striding info", "{}", "Striding information");
77  validParamList->set< LocalOrdinal >("Strided block id", -1, "Strided block id");
78 
79  // Domain GID offsets: list of offset values for domain gids (coarse gids) of tentative prolongator (default = 0).
80  // For each multigrid level we can define a different offset value, ie for the prolongator between
81  // level 0 and level 1 we use the offset entry domainGidOffsets_[0], for the prolongator between
82  // level 1 and level 2 we use domainGidOffsets_[1]...
83  // If the vector domainGidOffsets_ is too short and does not contain a sufficient number of gid offset
84  // values for all levels, a gid offset of 0 is used for all the remaining levels!
85  // The GIDs for the domain dofs of Ptent start with domainGidOffset, are contiguous and distributed
86  // equally over the procs (unless some reordering is done).
87  validParamList->set< std::string > ("Domain GID offsets", "{0}", "vector with offsets for GIDs for each level. If no offset GID value is given for the level we use 0 as default.");
88 
89  return validParamList;
90  }
91 
92  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
94  {
95  Input(currentLevel, "Aggregates");
96  Input(currentLevel, "Nullspace");
97  }
98 
99  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
101  {
102  // store striding map in internal variable
103  stridingInfo_ = stridingInfo;
104 
105  // try to remove string "Striding info" from parameter list to make sure,
106  // that stridingInfo_ is not replaced in the Build call.
107  std::string strStridingInfo; strStridingInfo.clear();
108  SetParameter("Striding info", ParameterEntry(strStridingInfo));
109  }
110 
111  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
113  {
114  FactoryMonitor m(*this, "Build", currentLevel);
115 
116  GlobalOrdinal domainGIDOffset = GetDomainGIDOffset(currentLevel);
117  BuildCoarseMap(currentLevel, domainGIDOffset);
118  }
119 
120  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
122  Level& currentLevel, const GlobalOrdinal domainGIDOffset) const
123  {
124  RCP<Aggregates> aggregates = Get< RCP<Aggregates> >(currentLevel, "Aggregates");
125  RCP<MultiVector> nullspace = Get< RCP<MultiVector> >(currentLevel, "Nullspace");
126 
127  GlobalOrdinal numAggs = aggregates->GetNumAggregates();
128  const size_t NSDim = nullspace->getNumVectors();
129  RCP<const Teuchos::Comm<int> > comm = aggregates->GetMap()->getComm();
130  const ParameterList & pL = GetParameterList();
131 
132  LocalOrdinal stridedBlockId = pL.get<LocalOrdinal>("Strided block id");
133 
134  // read in stridingInfo from parameter list and fill the internal member variable
135  // read the data only if the parameter "Striding info" exists and is non-empty
136  if(pL.isParameter("Striding info")) {
137  std::string strStridingInfo = pL.get<std::string>("Striding info");
138  if(strStridingInfo.empty() == false) {
139  Teuchos::Array<size_t> arrayVal = Teuchos::fromStringToArray<size_t>(strStridingInfo);
140  stridingInfo_ = Teuchos::createVector(arrayVal);
141  }
142  }
143 
144  CheckForConsistentStridingInformation(stridedBlockId, NSDim);
145 
146  GetOStream(Statistics2) << "domainGIDOffset: " << domainGIDOffset << " block size: " << getFixedBlockSize() << " stridedBlockId: " << stridedBlockId << std::endl;
147 
148  // number of coarse level dofs (fixed by number of aggregates and blocksize data)
149  GlobalOrdinal nCoarseDofs = numAggs * getFixedBlockSize();
150  GlobalOrdinal indexBase = aggregates->GetMap()->getIndexBase();
151 
152  RCP<const Map> coarseMap = StridedMapFactory::Build(aggregates->GetMap()->lib(),
153  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
154  nCoarseDofs,
155  indexBase,
156  stridingInfo_,
157  comm,
158  stridedBlockId,
159  domainGIDOffset);
160 
161  Set(currentLevel, "CoarseMap", coarseMap);
162  }
163 
164  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
166  Level& currentLevel) const
167  {
168  GlobalOrdinal domainGidOffset = Teuchos::ScalarTraits<GlobalOrdinal>::zero();
169 
170  std::vector<GlobalOrdinal> domainGidOffsets;
171  domainGidOffsets.clear();
172  const ParameterList & pL = GetParameterList();
173  if(pL.isParameter("Domain GID offsets")) {
174  std::string strDomainGIDs = pL.get<std::string>("Domain GID offsets");
175  if(strDomainGIDs.empty() == false) {
176  Teuchos::Array<GlobalOrdinal> arrayVal = Teuchos::fromStringToArray<GlobalOrdinal>(strDomainGIDs);
177  domainGidOffsets = Teuchos::createVector(arrayVal);
178  if(currentLevel.GetLevelID() < Teuchos::as<int>(domainGidOffsets.size()) ) {
179  domainGidOffset = domainGidOffsets[currentLevel.GetLevelID()];
180  }
181  }
182  }
183 
184  return domainGidOffset;
185  }
186 
187  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
189  const LocalOrdinal stridedBlockId, const size_t nullspaceDimension) const
190  {
191  // check for consistency of striding information with NSDim and nCoarseDofs
192  if (stridedBlockId == -1) {
193  // this means we have no real strided map but only a block map with constant blockSize "nullspaceDimension"
194  TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo_.size() > 1, Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): stridingInfo_.size() but must be one");
195  stridingInfo_.clear();
196  stridingInfo_.push_back(nullspaceDimension);
197  TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo_.size() != 1, Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): stridingInfo_.size() but must be one");
198 
199  } else {
200  // stridedBlockId > -1, set by user
201  TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockId > Teuchos::as<LO>(stridingInfo_.size() - 1), Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): it is stridingInfo_.size() <= stridedBlockId_. error.");
202  size_t stridedBlockSize = stridingInfo_[stridedBlockId];
203  TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockSize != nullspaceDimension , Exceptions::RuntimeError, "MueLu::CoarseMapFactory::Build(): dimension of strided block != nullspaceDimension. error.");
204  }
205  }
206 
207 } //namespace MueLu
208 
209 #endif /* MUELU_COARSEMAPFACTORY_DEF_HPP_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
virtual GlobalOrdinal GetDomainGIDOffset(Level &currentLevel) const
Extract domain GID offset from user data.
Timer to be used in factories. Similar to Monitor but with additional timers.
virtual void BuildCoarseMap(Level &currentLevel, const GlobalOrdinal domainGIDOffset) const
Build the coarse map using the domain GID offset.
Namespace for MueLu classes and methods.
RCP< const ParameterList > GetValidParameterList() const override
Return a const parameter list of valid parameters that setParameterList() will accept.
Print even more statistics.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Build(Level &currentLevel) const override
Build an object with this factory.
virtual void setStridingData(std::vector< size_t > stridingInfo)
setStridingData set striding vector for the domain DOF map (= coarse map), e.g. (2,1) for 2 velocity dofs and 1 pressure dof
Exception throws to report errors in the internal logical of the program.
virtual void CheckForConsistentStridingInformation(LocalOrdinal stridedBlockId, const size_t nullspaceDimension) const
void DeclareInput(Level &currentLevel) const override
Specifies the data that this class needs, and the factories that generate that data.