MueLu Version of the Day
Loading...
Searching...
No Matches
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
66namespace 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
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();
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
MueLu::DefaultGlobalOrdinal GlobalOrdinal
virtual void BuildCoarseMap(Level &currentLevel, const GlobalOrdinal domainGIDOffset) const
Build the coarse map using the domain GID offset.
virtual void CheckForConsistentStridingInformation(LocalOrdinal stridedBlockId, const size_t nullspaceDimension) const
RCP< const ParameterList > GetValidParameterList() const override
Return a const parameter list of valid parameters that setParameterList() will accept.
virtual GlobalOrdinal GetDomainGIDOffset(Level &currentLevel) const
Extract domain GID offset from user data.
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,...
void DeclareInput(Level &currentLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
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
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.