MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_NullspaceFactory_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_NULLSPACEFACTORY_DEF_HPP
47#define MUELU_NULLSPACEFACTORY_DEF_HPP
48
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_MultiVectorFactory.hpp>
51#include <Xpetra_BlockedMultiVector.hpp>
52#include <Xpetra_VectorFactory.hpp>
53
55#include "MueLu_Level.hpp"
56#include "MueLu_Monitor.hpp"
57#include "MueLu_MasterList.hpp"
58
59namespace MueLu {
60
61 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63 RCP<ParameterList> validParamList = rcp(new ParameterList());
64
65#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
66 SET_VALID_ENTRY("nullspace: calculate rotations");
67#undef SET_VALID_ENTRY
68 validParamList->set< std::string >("Fine level nullspace", "Nullspace", "Variable name which is used to store null space multi vector on the finest level (default=\"Nullspace\"). For block matrices also \"Nullspace1\" to \"Nullspace9\" are accepted to describe the null space vectors for the (i,i) block (i=1..9).");
69
70 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the fine level matrix (only needed if default null space is generated)");
71 validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the fine level null space");
72 validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory of the coordinates");
73
74 // TODO not very elegant.
75 // 1/20/2016: we could add a sublist (e.g. "Nullspaces" which is excluded from parameter validation)
76 validParamList->set< RCP<const FactoryBase> >("Nullspace1", Teuchos::null, "Generating factory of the fine level null space associated with the (1,1) block in your n x n block matrix.");
77 validParamList->set< RCP<const FactoryBase> >("Nullspace2", Teuchos::null, "Generating factory of the fine level null space associated with the (2,2) block in your n x n block matrix.");
78 validParamList->set< RCP<const FactoryBase> >("Nullspace3", Teuchos::null, "Generating factory of the fine level null space associated with the (3,3) block in your n x n block matrix.");
79 validParamList->set< RCP<const FactoryBase> >("Nullspace4", Teuchos::null, "Generating factory of the fine level null space associated with the (4,4) block in your n x n block matrix.");
80 validParamList->set< RCP<const FactoryBase> >("Nullspace5", Teuchos::null, "Generating factory of the fine level null space associated with the (5,5) block in your n x n block matrix.");
81 validParamList->set< RCP<const FactoryBase> >("Nullspace6", Teuchos::null, "Generating factory of the fine level null space associated with the (6,6) block in your n x n block matrix.");
82 validParamList->set< RCP<const FactoryBase> >("Nullspace7", Teuchos::null, "Generating factory of the fine level null space associated with the (7,7) block in your n x n block matrix.");
83 validParamList->set< RCP<const FactoryBase> >("Nullspace8", Teuchos::null, "Generating factory of the fine level null space associated with the (8,8) block in your n x n block matrix.");
84 validParamList->set< RCP<const FactoryBase> >("Nullspace9", Teuchos::null, "Generating factory of the fine level null space associated with the (9,9) block in your n x n block matrix.");
85
86 return validParamList;
87 }
88
89 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
91
92 const ParameterList & pL = GetParameterList();
93 std::string nspName = pL.get<std::string>("Fine level nullspace");
94
95 // only request "A" in DeclareInput if
96 // 1) there is not nspName (e.g. "Nullspace") is available in Level, AND
97 // 2) it is the finest level (i.e. LevelID == 0)
98 if (currentLevel.IsAvailable(nspName, NoFactory::get()) == false && currentLevel.GetLevelID() == 0)
99 Input(currentLevel, "A");
100
101 if( currentLevel.GetLevelID() == 0 &&
102 currentLevel.IsAvailable("Coordinates", NoFactory::get()) && // we have coordinates (provided by user app)
103 pL.get<bool>("nullspace: calculate rotations") ) { // and we want to calculate rotation modes
104 calculateRotations_ = true;
105 Input(currentLevel, "Coordinates");
106 }
107
108 if (currentLevel.GetLevelID() != 0) {
109 // validate nullspaceFact_
110 // 1) The factory for "Nullspace" (or nspName) must not be Teuchos::null, since the default factory
111 // for "Nullspace" is a NullspaceFactory
112 // 2) The factory for "Nullspace" (or nspName) must be a TentativePFactory or any other TwoLevelFactoryBase derived object
113 // which generates the variable "Nullspace" as output
114 TEUCHOS_TEST_FOR_EXCEPTION(GetFactory(nspName)==Teuchos::null, Exceptions::RuntimeError, "MueLu::NullspaceFactory::DeclareInput(): You must declare an existing factory which produces the variable \"Nullspace\" in the NullspaceFactory (e.g. a TentativePFactory).");
115 currentLevel.DeclareInput("Nullspace", GetFactory(nspName).get(), this); /* ! "Nullspace" and nspName mismatch possible here */
116 }
117 }
118
119 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121 FactoryMonitor m(*this, "Nullspace factory", currentLevel);
122
123 RCP<MultiVector> nullspace;
124
125 //TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.GetLevelID() != 0, Exceptions::RuntimeError, "MueLu::NullspaceFactory::Build(): NullspaceFactory can be used for finest level (LevelID == 0) only.");
126 const ParameterList & pL = GetParameterList();
127 std::string nspName = pL.get<std::string>("Fine level nullspace");
128
129 // get coordinates and compute mean of coordinates. (or centroid).
130
131 typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType coordinate_type;
132 typedef Xpetra::MultiVector<coordinate_type,LO,GO,NO> RealValuedMultiVector;
133 RCP<RealValuedMultiVector> Coords;
134 ArrayRCP<const coordinate_type> xvals, yvals, zvals;
135 Scalar cx, cy, cz;
136
137 cx = 0.; cy = 0.; cz = 0.;
138 if(calculateRotations_) {
139 Coords = Get< RCP<RealValuedMultiVector> >(currentLevel, "Coordinates");
140
141 cx = Coords->getVector(0)->meanValue();
142 if (Coords->getNumVectors() > 1)
143 cy = Coords->getVector(1)->meanValue();
144 if (Coords->getNumVectors() > 2)
145 cz = Coords->getVector(2)->meanValue();
146
147 xvals = Coords->getData(0);
148 if (Coords->getNumVectors() > 1)
149 yvals = Coords->getData(1);
150 if (Coords->getNumVectors() > 2)
151 zvals = Coords->getData(2);
152 }
153
154 if (currentLevel.GetLevelID() == 0) {
155
156 if (currentLevel.IsAvailable(nspName, NoFactory::get())) {
157 // When a fine nullspace have already been defined by user using Set("Nullspace", ...) or
158 // Set("Nullspace1", ...), we use it.
159 nullspace = currentLevel.Get< RCP<MultiVector> >(nspName, NoFactory::get());
160 GetOStream(Runtime1) << "Use user-given nullspace " << nspName << ": nullspace dimension=" << nullspace->getNumVectors() << " nullspace length=" << nullspace->getGlobalLength() << std::endl;
161 } else {
162 // "Nullspace" (nspName) is not available
163 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
164
165 // determine numPDEs
166 LocalOrdinal numPDEs = 1;
167 if(A->IsView("stridedMaps")==true) {
168 Xpetra::viewLabel_t oldView = A->SwitchToView("stridedMaps"); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
169 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap()) == Teuchos::null, Exceptions::BadCast, "MueLu::CoalesceFactory::Build: cast to strided row map failed.");
170 numPDEs = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap())->getFixedBlockSize();
171 oldView = A->SwitchToView(oldView);
172 }
173
174 LO nullspaceDim = numPDEs;
175
176 if(calculateRotations_) {
177 if (Coords->getNumVectors() > 1) nullspaceDim++;
178 if (Coords->getNumVectors() > 2) nullspaceDim += 2;
179 GetOStream(Runtime1) << "Generating nullspace with rotations: dimension = " << nullspaceDim << std::endl;
180 }
181 else GetOStream(Runtime1) << "Generating canonical nullspace: dimension = " << numPDEs << std::endl;
182
183 nullspace = MultiVectorFactory::Build(A->getDomainMap(), nullspaceDim);
184
185 RCP<BlockedMultiVector> bnsp = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(nullspace);
186 if(bnsp.is_null() == true) {
187 for (int i=0; i<numPDEs; ++i) {
188 ArrayRCP<Scalar> nsValues = nullspace->getDataNonConst(i);
189 int numBlocks = nsValues.size() / numPDEs;
190 for (int j=0; j< numBlocks; ++j) {
191 nsValues[j*numPDEs + i] = 1.0;
192 }
193 }
194 if (( (int) nullspaceDim > numPDEs ) && ((int) numPDEs > 1)) {
195 /* xy rotation */
196 ArrayRCP<Scalar> nsValues = nullspace->getDataNonConst(numPDEs);
197 int numBlocks = nsValues.size() / numPDEs;
198 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks != (int) xvals.size(), Exceptions::RuntimeError, "MueLu::NullspaceFactory::Build(): number of coordinates does not match ndofs/numPDEs.");
199 for (int j=0; j< numBlocks; ++j) {
200 nsValues[j*numPDEs + 0] = -(yvals[j]-cy);
201 nsValues[j*numPDEs + 1] = (xvals[j]-cx);
202 }
203 }
204 if (( (int) nullspaceDim == numPDEs+3 ) && ((int) numPDEs > 2)) {
205 /* xz rotation */
206 ArrayRCP<Scalar> nsValues = nullspace->getDataNonConst(numPDEs+1);
207 int numBlocks = nsValues.size() / numPDEs;
208 for (int j=0; j< numBlocks; ++j) {
209 nsValues[j*numPDEs + 1] = -(zvals[j]-cz);
210 nsValues[j*numPDEs + 2] = (yvals[j]-cy);
211 }
212 /* yz rotation */
213 nsValues = nullspace->getDataNonConst(numPDEs+2);
214 numBlocks = nsValues.size() / numPDEs;
215 for (int j=0; j< numBlocks; ++j) {
216 nsValues[j*numPDEs + 0] = -(zvals[j]-cz);
217 nsValues[j*numPDEs + 2] = (xvals[j]-cx);
218 }
219 }
220 /*
221 // Scale columns to match what Galeri does. Not sure how important this is as QR factorization
222 // when creating ptent should take care of scaling issues ... but leaving it just in case.
223 if ( (int) nullspaceDim > numPDEs ) {
224 Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> norms2(nullspaceDim);
225 nullspace->norm2(norms2);
226 Teuchos::Array<Scalar> norms2scalar(nullspaceDim);
227 for (int i = 0; i < nullspaceDim; i++)
228 norms2scalar[i] = norms2[0] / norms2[i];
229 nullspace->scale(norms2scalar);
230 }
231 */
232
233 } else {
234 fillNullspaceVector(bnsp,numPDEs, xvals, yvals, zvals, nullspaceDim, cx, cy, cz);
235 }
236 } // end if "Nullspace" not available
237 } else {
238 // on coarser levels always use "Nullspace" as variable name, since it is expected by
239 // tentative P factory to be "Nullspace"
240
241 nullspace = currentLevel.Get< RCP<MultiVector> >("Nullspace", GetFactory(nspName).get()); /* ! "Nullspace" and nspName mismatch possible here */
242
243 }
244
245 // provide "Nullspace" variable on current level (used by TentativePFactory)
246 Set(currentLevel, "Nullspace", nullspace);
247
248 } // Build
249
250 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
251 void NullspaceFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::fillNullspaceVector(const RCP<BlockedMultiVector>& nsp, LocalOrdinal numPDEs, ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> xvals, ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> yvals, ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> zvals, LocalOrdinal nullspaceDim, Scalar cx, Scalar cy, Scalar cz) const {
252 RCP< const BlockedMap> bmap = nsp->getBlockedMap();
253
254 for(size_t r = 0; r < bmap->getNumMaps(); r++) {
255 Teuchos::RCP<MultiVector> part = nsp->getMultiVector(r);
256 Teuchos::RCP<BlockedMultiVector> bpart = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(part);
257 if(bpart.is_null() == true) {
258 for (int i=0; i<numPDEs; ++i) {
259 ArrayRCP<Scalar> nsValues = part->getDataNonConst(i);
260 int numBlocks = nsValues.size() / numPDEs;
261 for (int j=0; j< numBlocks; ++j) {
262 nsValues[j*numPDEs + i] = 1.0;
263 }
264 }
265 if ( (int) nullspaceDim > numPDEs ) {
266 /* xy rotation */
267 ArrayRCP<Scalar> nsValues = part->getDataNonConst(numPDEs);
268 int numBlocks = nsValues.size() / numPDEs;
269 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks != (int) xvals.size(), Exceptions::RuntimeError, "MueLu::NullspaceFactory::fillNullspaceVector(): number of coordinates does not match ndofs/numPDEs.");
270 for (int j=0; j< numBlocks; ++j) {
271 nsValues[j*numPDEs + 0] = -(yvals[j]-cy);
272 nsValues[j*numPDEs + 1] = (xvals[j]-cx);
273 }
274 }
275 if ( (int) nullspaceDim == numPDEs+3 ) {
276 /* xz rotation */
277 ArrayRCP<Scalar> nsValues = part->getDataNonConst(numPDEs+1);
278 int numBlocks = nsValues.size() / numPDEs;
279 for (int j=0; j< numBlocks; ++j) {
280 nsValues[j*numPDEs + 1] = -(zvals[j]-cz);
281 nsValues[j*numPDEs + 2] = (yvals[j]-cy);
282 }
283 /* yz rotation */
284 nsValues = part->getDataNonConst(numPDEs+2);
285 numBlocks = nsValues.size() / numPDEs;
286 for (int j=0; j< numBlocks; ++j) {
287 nsValues[j*numPDEs + 0] = -(zvals[j]-cz);
288 nsValues[j*numPDEs + 2] = (xvals[j]-cx);
289 }
290 }
291 /*
292 // Scale columns to match what Galeri does. Not sure that this is necessary as the qr factorizatoin
293 // of the tentative prolongator also takes care of scaling issues. I'm leaving the code here
294 // just in case.
295 if ( (int) nullspaceDim > numPDEs ) {
296 Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> norms2(nullspaceDim);
297 nullspace->norm2(norms2);
298 Teuchos::Array<Scalar> norms2scalar(nullspaceDim);
299 for (int i = 0; i < nullspaceDim; i++)
300 norms2scalar[i] = norms2[0] / norms2[i];
301 nullspace->scale(norms2scalar);
302 }
303 */
304 } else {
305 // call this routine recursively
306 fillNullspaceVector(bpart,numPDEs,xvals, yvals, zvals, nullspaceDim, cx, cy, cz);
307 }
308 }
309 }
310
311} //namespace MueLu
312
313#endif // MUELU_NULLSPACEFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
Exception indicating invalid cast attempted.
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.
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()
void Build(Level &currentLevel) const
Build an object with this factory.
void fillNullspaceVector(const RCP< BlockedMultiVector > &nsp, LocalOrdinal numPDEs, ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > xvals, ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > yvals, ArrayRCP< const typename Teuchos::ScalarTraits< Scalar >::coordinateType > zvals, LocalOrdinal nullspaceDim, Scalar cx, Scalar cy, Scalar cz) const
Helper function to recursively fill BlockedMultiVector with default null space vectors.
RCP< const ParameterList > GetValidParameterList() const
Define valid parameters for internal factory parameters.
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Namespace for MueLu classes and methods.
@ Runtime1
Description of what is happening (more verbose)